Molecular Vision 2020; 26:117-134
Received 17 September 2019 | Accepted 03 March 2020 | Published 05 March 2020
Chao Geng, Yahong Li, Fang Guo, Jindan Wang, Yiyun Yue, Kewen Zhou, Ruihua Wei, Yan Zhang
Tianjin Key Laboratory of Retinal Functions and Diseases, Eye Institute and School of Optometry, Tianjin Medical University Eye Hospital, Tianjin, China
Correspondence to: Yan Zhang, Tianjin Key Laboratory of Retinal Functions and Diseases, Eye Institute and School of Optometry, Tianjin Medical University Eye Hospital, Tianjin, 300384, China; Phone: +86 2286428863; FAX: +86 2286428777; email: email@example.com
Purpose: To detect the differential expression of long non-coding RNAs (lncRNAs) in the ocular posterior poles of two guinea pig myopia models and explore the pathogenic role of lncRNAs in myopia.
Methods: Form-deprived myopia (FDM) and lens-induced myopia (LIM) models were induced in guinea pig right eyes by wearing a translucent latex balloon head mask and a −10.00 diopter (D) lens, respectively. Ocular biometric parameters were measured biweekly. At 6 weeks after the induction of myopia, the guinea pig eyeballs were processed for hematoxylin and eosin staining to examine the ocular morphology. The ocular posterior poles from the normal control, FDM, and LIM groups were collected to analyze the differential expression of lncRNAs between the groups with high-throughput RNA sequencing (RNA-seq). Further, the lncRNA-mRNA colocation network was delineated to predict the functions of the differentially expressed lncRNAs. Last, Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed on the colocated mRNAs of the differentially expressed lncRNAs. Additionally, the expression of the most differentially expressed lncRNAs in the myopia-induced eyes and the contralateral eyes was validated with quantitative real-time PCR (qPCR).
Results: Compared with the normal controls and the contralateral eyes, the myopia-induced eyes in the FDM and LIM groups exhibited decreased scleral and choroidal thicknesses, reduced refraction, and increased ocular axial length but without changes in the corneal curvature radius at 6 weeks after myopia was induced. RNA-seq showed that 372 and 247 lncRNAs were differentially expressed in the FDM and LIM groups, respectively, in comparison to the normal counterparts. There were 380 differentially expressed lncRNAs in the LIM group compared to the FDM group. The GO and KEGG analyses showed that the colocated mRNAs of the differentially expressed lncRNAs were enriched in cellular components such as the extracellular matrix (ECM) structural constituent; in molecular functions such as kinase activity, metabolism, and growth; as well as in pathways including ECM-receptor interaction, glycosaminoglycan degradation, and mucin type O-Glycan biosynthesis. The expression patterns of the selected lncRNAs were verified with qPCR.
Conclusions: High-throughput RNA-seq revealed previously undescribed lncRNA expression profiling in guinea pig FDM and LIM models. These results may shed light on the molecular pathogenesis of myopia and provide clues for interventional targets for this highly prevalent visual disorder.
Myopia is one of the most common types of refractive errors. In recent years, the incidence of myopia has increased dramatically, particularly in Asia . A study predicted that nearly half of the global population would develop myopia by 2050, and individuals with high myopia, whose spherical equivalent refractive error is more than −6.0 diopters (D), would account for 10% of the global population [2,3]. Therefore, myopia has become a major public health issue worldwide. Moreover, high myopia can incur a series of ocular pathologies, including cataract, glaucoma, retinal detachment, myopic choroidal neovascularization, and posterior scleral staphyloma, thus leading to visual dysfunctions and irreversible vision loss [4,5].
The exact pathogenic mechanism underlying myopia remains unclear. Recent studies revealed that in addition to long-known genetic factors, environmental factors and interactions between the two play important roles in the initiation and progression of myopia . Epigenetics refers to the generation of hereditable phenotypic changes under the influence of environmental factors without changes in the genetic code. Epigenetics may alter gene expression through DNA methylation, posttranslational modifications of histone proteins, chromosomal remodeling, and regulation of non-coding RNAs, such as microRNA (miRNA) and long non-coding RNAs (lncRNAs) [6-8]. The role of miRNA in the pathogenesis of myopia has been examined with high-throughput approaches. For instance, Metlapally et al.  used microarrays to analyze genome-wide miRNA and mRNA expression changes in scleral tissues of a murine model of form-deprived myopia (FDM), suggesting the involvement of miRNAs in scleral remodeling and eyeball growth during the progression of myopia. In addition, a miRNA microarray of a mouse model of lens-induced myopia (LIM) suggested that miRNAs with overlapping expression changes in different eye tissues during the progression of myopia may serve as biomarkers for predicting myopia . However, the role of lncRNAs in the pathogenesis and development of myopia remains unknown.
RNA sequencing (RNA-seq) is a high-throughput sequencing technology of whole transcriptome (including mRNAs and non-coding RNAs) that can be used to determine gene transcriptional structures and RNA expression levels, as well as to study complex transcripts and discover novel transcripts . With the advances in RNA-seq technology, an expanding body of lncRNAs has been discovered. Studies have shown that lncRNAs may regulate RNA alternative splicing, modulate protein activities, participate in transcriptional regulation, and affect post-transcriptional modifications, ultimately leading to the occurrence and development of diseases [12,13]. RNA-seq has been applied to explore the molecular mechanism of myopia. For example, Srinivasalu et al.  used this technology to analyze gene expression changes in different scleral regions of FDM guinea pigs. The results revealed differential gene expression in the sclera around the optic nerve, which may be related to axis elongation in myopia. In the present study, the guinea pig FDM and LIM models were established and characterized. RNA-seq technology was employed to compare lncRNA expression profiling in the ocular posterior pole, mainly composed of the neural parenchyma and vasculature in the retina, choroid, and sclera, among the normal control, FDM, and LIM groups. Furthermore, Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of mRNAs adjacent to the differentially expressed lncRNAs were performed to predict the regulatory functions of the lncRNAs. This is the first report on the differential expression of lncRNAs in rodent models of myopia. The results provide clues to the possible roles of lncRNAs in the pathogenesis of myopia and interventional targets to delay the progression of this rapidly developing epidemic.
One hundred male 3-week-old colored guinea pigs (body weight 150–180 g) were purchased from Aochen experimental animal farm (Danyang, China). The animals were subjected to ocular and optometric examinations upon arrival, and those with eye diseases, myopia, and anisometropia more than 1.5 D were excluded. The rest were maintained in a temperature- and humidity-conditioned facility (25 ± 1.0 °C with relative humidity of 40–70%) at Tianjin Medical University Eye Institute. They were fed food, water, and vegetables ad libitum. The illumination was under 12-h:12-h light-dark cycles. The study adhered to the ARVO Statement for Use of Animals in Research, and all experimental procedures were approved by the Laboratory Animal Care and Use Committee of Tianjin Medical University (permission number: SYXK2009–0001).
The guinea pigs were divided into three groups: the normal control group, the FDM group, and the LIM group (n = 17/group). All eyes in the normal control group were left untreated. The guinea pigs in the FDM group wore a head mask made of translucent latex balloon as previously described . The right eyes were fully covered but without extra pressure on the corneas and eyelids. The left eyes, noses, mouths, and both ears were uncovered. During the entire study, care was taken to ensure that the right eyes were continuously and completely covered but free to blink. Myopia was induced in the guinea pigs in the LIM group by the wearing of a −10.00 D lens according to the literature . Specifically, another type of head mask, also made of translucent latex balloon, was worn by each animal to expose the eyes, ears, and noses. Nylon buckles were sewn on the head mask around the right eye. The lens was fixed on the head mask in front of the right eye by the nylon buckles. Special care was taken not to compress the cornea and eyelid of the eye, and to ensure that the optical center of the lens matched the center of the pupil. The left eyes of the animals in the FDM and LIM groups did not receive myopic induction and were used as contralateral controls.
Following the induction of myopia, biometric parameters in the eye, including the refraction, axial length, and radius of corneal curvature, were measured every 2 weeks. Briefly, refraction was measured with a streak retinoscope (YZ24; Six six vision technology Co., Ltd., Suzhou, China) in darkness after dilatation of the pupils with 0.5% compound tropicamide. The working distance of the measurement was 50 cm. The mean refractions of the horizontal and vertical meridians were deemed the final equivalent spherical degrees. The measurement was repeated three times for each eye, including the myopia-induced right eye and the contralateral left eye unless otherwise stated, and the average was recorded.
Ophthalmic A-scan ultrasonography (11 MHz, AVISO Echograph class I-Type Bat; Quantel Medical, Clermont-Ferrand, France) was performed to measure the axial length following topical anesthesia of the ocular surface with 0.4% bupivacaine hydrochloride eye drops. The ultrasonic probe was aimed at the center of the pupil and touched the cornea apex vertically to record stable waveforms. The axial length was recorded as the distance from the anterior apex of the cornea to the anterior surface of the retina. Each eye was measured five times, and the average was recorded.
Based on the literature , the corneal curvature radius was measured with a keratometer (OM-4; Topcon, Tokyo, Japan). A +8.00 D lens was placed in the front of the keratometer to amplify the guinea pig’s cornea. The measuring window was in parallel to the eye to be measured. On the keratometer handle, the horizontal and vertical adjustment knobs were appropriately adjusted so that the three light circles were clearly imaged on the cornea and coincided. The horizontal and vertical readings were averaged and then multiplied by the correction coefficient of 0.451 to obtain the corneal curvature radius of the measured eye. Each eye was measured three times, and the average was recorded.
At week 6 following the induction of myopia, all animals were euthanized with an intraperitoneal injection of an overdosed pentobarbital sodium (100 mg/kg). The eyes of the animals were rapidly collected and dissected from the muscle, fascia, and other tissues attached to the outer wall of the eyeballs. Some eyeballs were processed for morphology analysis. For RNA-seq and subsequent quantitative real-time PCR (qPCR), the dissected eyeballs were cut at limbi with a sharp eyebrow trimmer that had been treated with RNaseZap (Thermo Fisher Scientific, Waltham, MA). The resulting ocular posterior pole tissues were snap frozen in liquid nitrogen and stored at –80 °C until further usage.
The guinea pig eyeballs from the normal control, FDM, and LIM groups (n = 6/group) were fixed in 4% paraformaldehyde, paraffin embedded, and serially sectioned (3 µm/section). Ten sections at the comparable position of each eyeball were selected and stained with hematoxylin and eosin (H&E). Pictures were captured under the bright-field of a BX51 microscope (Olympus Optical Co. Ltd., Tokyo, Japan) using identical optical parameters. The thicknesses of the choroid and the sclera in each picture were measured using the cellSens Standard electronic system (Olympus Optical Co. Ltd.). Briefly, three different positions (left third, middle third, and right third) on the choroid or sclera were selected in each picture, and the thickness at each position was measured and labeled by the software. The averaged thickness of the three positions was recorded and represented the thickness of the choroid or sclera on that picture.
Total RNA of the stored posterior poles of the eyes, mainly composed of the parenchyma and vasculature of the retina, choroid, and sclera, in the normal control, FDM, and LIM groups (n = 11/group) was extracted using a GeneJET RNA Extraction Kit (Thermo Fisher Scientific) according to the manufacturer’s protocol. RNA degradation and contamination were monitored on 1% agarose gels. A NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA), a Qubit RNA Assay Kit (Thermo Fisher Scientific, Carlsbad, CA), and an RNA Nano 6000 Assay Kit (Agilent Technologies, Santa Clara, CA) were used to measure the concentration and check the purity and the integrity of the extracted RNA, respectively.
An Epicenter Ribo-Zero rRNA Removal Kit (Epicenter, Charlotte, NC) was used to remove the rRNA from the total RNA extracted from the right eye posterior poles of the normal control, FDM, and LIM groups (n = 3/group). The rRNA-depleted RNA was precipitated with ethanol and used to generate sequencing libraries with a NEBNext Ultra Directional RNA Library Prep Kit (New England Biolabs, Ipswich, MA) following the manufacturer’s procedures. Finally, the library quality was assessed on an Agilent Bioanalyzer 2100 system, and the libraries were sequenced on a HiSeq 4000 platform (Illumina, San Diego, CA).
Clean data (clean reads) were obtained by deleting the reads containing adapters or multiple unidentifiable nucleotides, and low-quality reads from the raw reads. The clean reads were then aligned to the reference genome using HISAT2 v2.0.4 software. Picard tools v1.96 and Sequence Alignment/Map (SAM) tools v0.1.18 were used to sort and mark the duplicated reads and record the Binary Alignment/Map (BAM) alignments in each sample, thus arranging the chromosome coordinates and removing the duplicate reads. GATK2 software was used to perform single nucleotide polymorphism (SNP) calling and filter the original results.
The mapped reads for each sample were assembled with StringTie software (v1.3.1), and the transcripts were split using the HISAT2 alignment results. Then, the splicing variants of each transcript were merged with Cuffmerge software, and the transcripts without a definite direction of a positive or negative chain were excluded. Then the lncRNAs were screened in five steps: First, the transcripts should have more than two exons. Second, the transcripts should be longer than 200 nt. Third, the transcripts overlapping with the annotated exon regions in database were excluded with Cuffcompare software, and the lncRNAs overlapping with the exon regions of the spliced transcripts were included as the annotated lncRNAs in the database for subsequent analyses. Fourth, the expression amplitude of each transcript was calculated as fragments per kilobase million (FPKM) with Cuffquant software, and the transcripts with FPKM not less than 0.5 were selected. Fifth, the transcripts predicted to be without protein-coding capability by the four tools, including Coding-Non-Coding-Index, Coding Potential Calculator, Pfam-sca, and phylogenetic codon substitution frequency, were eventually deemed lncRNAs. Cuffdiff (v2.1.1) software was used to calculate the FPKM of the lncRNAs in each sample. Cuffdiff used a negative binomial distribution model to generate the statistical criteria for differential expression: lncRNAs with a Q value of less than 0.05 were deemed differentially expressed. All the differentially expressed lncRNAs were then sorted in descending order of the fold changes between the two experimental groups.
According to the cis-regulatory effects of lncRNAs on the expression of adjacent genes, the main function of the differentially expressed lncRNAs was predicted by the colocated target genes, which refer to the coding genes within 10 to 100 kilobases (kb) upstream and downstream of the lncRNA of interest. The GOseq R package and the KEGG Orthology Based Annotation System (KOBAS) were used to examine statistical enrichment of the lncRNA target genes in the GO and KEGG pathway analyses, respectively.
According to the RNA-seq results, the five most differentially expressed lncRNAs between the experimental groups were selected based on the statistical significance and fold changes. The expression patterns of these lncRNAs in the normal control group and the myopia-induced eyes and the contralateral untreated eyes in the FDM and LIM groups (n = 8/group) were verified with qPCR. Briefly, 1 μg of the extracted total RNAs were reverse transcribed into cDNAs using random hexamers and reagents in a RevertAid RT Reverse Transcription Kit (Thermo Fisher Scientific, Waltham, MA) according to the manufacturer’s protocol. qPCR was conducted in a 96-well format with the HT7900 Real-Time PCR system (Applied Biosystem, Foster City, CA). The cDNA content of the target genes was normalized to the internal standard glyceraldehyde 3-phosphate dehydrogenase (Gapdh) gene. The reaction mixture comprised the 4 μl cDNA template, gene-specific qPCR primers (Table 1), and 12.5 μl EvaGreen 2X Master Mix (Applied Biologic Materials Inc., Richmond, Canada). One microliter of cDNA from each sample was pooled, serially diluted, and used as a template to generate a standard curve between the Ct values of each target gene or the Gapdh gene and the logarithms of the cDNA template concentrations. The standard curves demonstrated similar priming efficiencies between each target gene and the Gapdh gene. The standard curves also served as positive controls for qPCR, and the reactions using water as the template served as negative controls. The qPCR program consisted of preincubation at 50 °C for 2 min, denaturation at 95 °C for 10 min, 40 cycles of denaturation at 95 °C for 15 s, and extension at 60 °C for 1 min. A dissociation stage was added to check the amplicon specificity. The relative expression levels of the lncRNAs were analyzed using the comparative threshold cycle (2−∆∆Ct) method.
Statistical Program for Social Sciences 20.0 software (IBM SPSS Inc., New York, NY) was used for statistical analysis of the biometric parameters, scleral and choroidal thicknesses, and qPCR results. The data of the three experimental groups were examined with the Shapiro-Wilk test to confirm normal distribution. All data are expressed as mean ± standard error of the mean (SEM). The homogeneity of variance was confirmed with the Levene test. The ocular biometric parameters, sclera and choroid thicknesses, and the abundance of the selected differentially expressed lncRNAs in the myopia-induced right eyes and the untreated contralateral eyes were compared using the paired t test. The overall comparisons of these indices in the right eyes among the experimental groups were analyzed with two-way analysis of variance (ANOVA) or one-way ANOVA, and the pairwise comparisons were performed with the Tukey post hoc test. A p value of less than 0.05 was considered statistically significant.
There was no statistically significant difference in the refraction, axial length, and corneal curvature radius between the left and right eyes in the three groups before the induction of myopia (n = 17/group, paired T test, all p>0.05, Figure 1A–F). At 2, 4, and 6 weeks following the induction of myopia, the diopters of the induced right eyes were statistically significantly decreased, and the axial length was statistically significantly increased compared to the untreated left eyes in the FDM and LIM groups (n = 17/group, paired T test, all p<0.001, Figure 1A,B,D,E). However, the corneal curvature radius between the right eye and the left eye in the two myopia groups did not exhibit statistically significant differences throughout the experimental period (n = 17/group, paired T test, all p>0.05, Figure 1C,F). Compared with the normal control group, the diopters of the induced right eyes in the FDM and LIM groups decreased statistically significantly (n = 17/group, two-way ANOVA and Tukey post hoc, all p<0.001, FDM versus NOR, LIM versus NOR, Figure 1G), whereas the axial length increased statistically significantly (n = 17/group, two-way ANOVA and Tukey post hoc, all p<0.001, FDM versus NOR, LIM versus NOR, Figure 1H) over time. Furthermore, the extent of the myopia in the right eyes of the FDM group was statistically significantly greater than that in the corresponding eyes of the LIM group at all time points examined (n = 17/group, two-way ANOVA and Tukey post hoc, for refraction and axial length, all p<0.05, FDM versus LIM, Figure 1G,H). In addition, no statistically significant difference was found in the corneal curvature radius of the right eyes among the three experimental groups (n = 17/group, two-way ANOVA and Tukey post hoc, all p>0.05, Figure 1I), indicating that the form-deprived and lens-induced myopia interventions, at least in this experimental paradigm, did not affect the corneal curvature radius.
The H&E staining results revealed the normal thickness of the posterior sclera in the normal control group. In the normal control group, the scleral collagen fibers were densely arranged and without a break (Figure 2A). However, in the right eyes of the FDM and LIM groups, the sclera thickness was thinner, and the collagen fibers were sparse and frequently broken (Figure 2B,C). Quantitatively, the thickness of the posterior sclera of the right eyes in the FDM and LIM groups was 51.91 ± 5.910 μm and 55.06 ± 9.540 μm, respectively. There was no statistically significant difference between the two groups (n = 6/group, one-way ANOVA and Tukey post hoc, p>0.05, FDM versus LIM, Figure 2D). In contrast, the posterior sclera thickness in the normal control group was 78.09 ± 9.660 μm, statistically significantly thicker than that of the right eyes in the FDM and LIM groups (n = 6/group, one-way ANOVA and Tukey post hoc, all p<0.001, FDM versus NOR, LIM versus NOR, Figure 2D). Similarly, the thickness of the choroid in the normal control, FDM, and LIM right eyes was 63.42 ± 10.96 μm, 36.50 ± 6.580 μm, and 43.76 ± 12.330 μm, respectively, with the thicknesses of the FDM and LIM right eyes statistically significantly less than that of the normal controls (n = 6/group, one-way ANOVA and Tukey post hoc, both p<0.001, FDM versus NOR, LIM versus NOR, Figure 2D). However, unlike the sclera thickness, the choroid thickness of the FDM right eyes was statistically significantly decreased compared to that of the LIM right eyes (n = 6/group, one-way ANOVA and Tukey post hoc, p<0.001, FDM versus LIM, Figure 2D). The scleral and choroidal thicknesses of the induced right eyes were all statistically significantly decreased in comparison to those of the untreated left eyes in both myopia groups (n = 6/group, paired T test, all p<0.001, Figure 2B,C,E–G).
After RNA-seq, the raw reads of the normal control, FDM, and LIM groups were 102,431,862, 109,706,556, and 101,301,838, respectively. Clean reads and clean bases were obtained after low-quality reads were filtered from the raw reads. Detailed data quality results are shown in Table 2.
The top 30 chromosomes harboring the largest number of SNPs are shown in Figure 3A. The violin plot shows the distribution trends of the SNPs in the top 30 chromosomes (Figure 3B). The SNPs in the normal control group are mainly distributed on the chromosomes with 6,000 to 11,000 SNPs. The SNPs in the FDM group are mainly enriched on the chromosomes with 5,000 SNPs. Most of the SNPs in the LIM group are found on the chromosomes with 6,000 to 13,000 SNPs (Figure 3B).
Based on the spliced transcript results, as well as the structural and functional characteristics of lncRNAs, a rigorous screening process was set up. After five screening steps (Figure 3C), the selected lncRNAs were used for the subsequent analyses. The proportion of different types of lncRNAs, such as long intergenic non-coding RNA (lincRNAs), antisense RNAs, and intronic lncRNAs, is illustrated (Figure 3D).
Hierarchical clustering analysis showed that a total of 703 lncRNAs were differentially expressed (negative binomial distribution model, Q < 0.05), and the lncRNAs were clustered in three groups (Figure 4). Three hundred and seventy-two lncRNAs were differentially expressed between the normal control and FDM groups, of which 228 were upregulated and 144 downregulated (Figure 5A). The top 20 differentially expressed lncRNAs selected according to statistical significance and fold changes, such as LNC_000906, rna10561, rna39119, and LNC_003158, between the normal control and FDM groups are listed in Table 3. There were 247 differentially expressed lncRNAs between the normal control group and the LIM group, including 119 upregulated and 128 downregulated lncRNAs (Figure 5B). The top 20 differentially expressed lncRNAs in the LIM group, compared to the normal control group, are listed in Table 4. A total of 380 lncRNAs (146 upregulated and 234 downregulated) showed differential expression between the LIM group and the FDM group (Figure 5C). In Table 5, the top 20 differentially expressed lncRNAs are shown, including LNC_003158, LNC_000049, rna39119, and LNC_002414. In addition, 287 lncRNAs were differentially expressed between only the FDM and normal control groups, 162 between only the LIM and normal control groups; whereas 85 were differentially expressed in both myopia groups compared to the normal control group (Figure 5D). Among these 85 lncRNAs, 24 were annotated lncRNAs, and 61 were novel. Furthermore, 80% of them (68 of the 85 lncRNAs) exhibited similar trends in both myopia groups (35 upregulated and 33 downregulated) in comparison to the normal control group (Figure 5E,F).
LncRNAs can regulate the expression of neighboring protein-coding genes in a cis-manner; therefore, the protein-coding genes colocated with the differentially expressed lncRNAs were designated as the target genes of the lncRNAs. The target genes were analyzed to predict the possible functions of the lncRNAs under the pathological conditions of FDM and LIM. The target genes of the most upregulated and downregulated lncRNAs in the FDM and LIM groups are illustrated in Figure 6.
GO enrichment analysis was conducted to analyze the target gene functions in three categories, including cell components, biologic process, and molecular function. Compared with the normal control group, the target genes in the FDM group were enriched in 1,972 terms, which mainly involved the cellular components, such as extracellular matrix (ECM), collagen, the cytoplasmic membrane, and ribosomes. Moreover, the target genes in this myopia group participated in biologic processes such as cobalamin, peptide, vitamin transport, cell division, biologic metabolism, and biologic regulation. The target genes were also involved in molecular functions, including kinase activity, protein binding, and ion binding (Figure 7A). However, the target genes in the LIM group, in comparison to the normal control group, were enriched to 1,531 terms, and contributed to cell components including ribosomes, cytoplasm, protein complex, and plasma membrane. The target genes in the LIM group participated in biologic processes such as growth, development, signal regulation, immunity, biosynthesis, and metabolism. At the molecular level, the target genes functioned in insulin-like growth factor binding, kinase activity, and nucleoside binding (Figure 7B). As for the LIM group relative to the FDM group, the target genes were mainly involved in primary metabolic processes and contributed to the molecular function of purine ribonucleoside and its triphosphate binding (Figure 7C).
KEGG pathway analysis revealed the top ten signaling pathways enriched with the target genes of the differentially expressed lncRNAs (Figure 8). The FDM group participated in ECM-receptor interaction, glycosaminoglycan degradation, the PI3K-Akt signaling pathway, histidine metabolism, glutamatergic synapse, the Notch signaling pathway, and so forth. Among these pathways, ECM-receptor interaction and glycosaminoglycan degradation were the most significantly enriched signaling pathways (Figure 8A). The lncRNA target genes in the LIM group were mainly enriched in antigen processing and presentation, phagosome, mucin type O-Glycan biosynthesis, the AMPK signaling pathway, etc. (Figure 8B). Taken together, the target genes in the two myopia groups were statistically significantly enriched in ECM-receptor interaction and mucin type O-Glycan biosynthesis pathways. When the LIM group was compared to the FDM group, the differentially expressed lncRNA target genes were enriched in the endocytosis, ECM-receptor interaction, and glycerolipid metabolism pathways (Figure 8C).
Five differentially expressed lncRNAs with the greatest fold changes between the experimental groups were selected, and the expression patterns were examined with qPCR. The results demonstrated that the expression of lnc000906 and lnc000049 was significantly upregulated in the right eyes of the FDM and LIM groups, respectively, compared to those of the normal control group n = 8/group, one-way ANOVA and Tukey post hoc, p<0.05, for FDM-R versus NOR; p<0.01, for LIM-R versus NOR; Figure 9A,B). In addition, the abundance of these two lncRNAs in the right eyes was statistically significantly higher than in the contralateral eyes (n = 8/group, paired t test, both p<0.05, for FDM-R versus FDM-L, LIM-R versus LIM-L; Figure 9A,B). Conversely, the expression of lnc0031538 and lnc002905 was downregulated in the FDM and LIM right eyes, respectively, compared to the normal and contralateral counterparts (n = 8/group, all p<0.05, one-way ANOVA and Tukey post hoc for FDM-R versus NOR, and LIM-R versus NOR; paired t test for FDM-R versus FDM-L, and LIM-R versus LIM-L; Figure 9C,D). Moreover, lnc0031538 and rna39119 exhibited differential expression between the right eyes of the two myopia groups (n = 8/group, one-way ANOVA and Tukey post hoc, p<0.01, for lnc0031538, FDM-R versus LIM-R; p<0.05, for rna39119, FDM-R versus LIM-R; Figure 9C,E). The qPCR results verified those of RNA-seq (Tables 3, 4, and 5), indicating the reliability of this high-throughput gene expression analysis technology.
Numerous epidemiological studies have shown that environmental factors, such as close reading, high intensity study, and excessive use of electronic devices, have a close relationship with the occurrence and development of myopia . However, the mechanism underpinning the environmental factors in the pathogenesis of myopia remains elusive. A few clues have implicated that the environmental factors might influence the initiation and progression of myopia through epigenetic regulatory mechanisms. For instance, Vishweswaraiah et al.  found a significantly increased number of hypermethylated genes in the peripheral blood of patients with high myopia. Furthermore, lncRNAs, as an important regulatory mechanism of epigenetics, have not been reported to be associated with myopia pathogenesis. Therefore, in this study, we employed high-throughput RNA-seq technology and for the first time, analyzed the differential expression of lncRNAs in the two well-recognized guinea pig models of myopia, the FDM and LIM models, in the hope that this study might shed light on the pathogenesis of myopia and identify potential interventional targets for this highly prevalent visual disorder.
Studies have shown that form deprivation and negative lens can disrupt the formation of clear images on the retina, thus promoting excessive ocular axis growth and scleral and choroidal thinning [19-22]. Consistent with these studies, the present results showed that under form deprivation and lens induction, the refraction was dramatically decreased, the ocular axial length significantly extended, and the sclera and choroid thicknesses significantly reduced compared with those of the normal controls. However, the molecular mechanisms underlying FDM and LIM may be different. Form deprivation is an open-loop system, and myopia continues to progress during deprivation. However, lens-induced optical defocusing is a closed-loop system. Once the stimulus is eliminated, the axial length counteracts the hyperopia defocus induced by the negative lens and ceases the excessive growth of the eyeball . One aim of the present study was to perform RNA-seq analysis on the differential expression of lncRNAs in the two myopia models and reveal their underlying molecular mechanisms.
The present study showed that lncRNAs were aberrantly expressed in the FDM and LIM groups. A total of 703 differentially expressed lncRNAs were further analyzed. It has been reported that lncRNAs can regulate the expression of adjacent genes at the transcriptional level in a cis-acting manner . Therefore, the functions of the differentially expressed lncRNAs may be related to the adjacent protein-coding genes (designated as the target genes of the lncRNAs), and the functional analyses of the target mRNAs by GO and signaling pathways can predict lncRNA functions.
In the FDM group, the target genes of the significantly upregulated lncRNAs included Sgsh, Col4a1, and Col4a2, which may be involved in the development of myopia. The mutation in Sgsh has been reported to be associated with mucopolysaccharidosis type IIIA disease, which is characterized by the accumulation of glycosaminoglycans (GAGs) in human organs . GAGs are unbranched polysaccharides distributed in scleral ECM. GAGs may interact with collagen fiber in ECM, fill the gaps in the fiber network, and bind to water molecules, thus maintaining scleral mechanical properties [26,27]. Previous studies have revealed reduced sulfated GAGs in sclera of myopic eyes . Consistently, the pathway analysis of this study showed that the lncRNA target genes in the FDM group were significantly enriched in the glycosaminoglycan degradation pathway. Therefore, we speculate that abnormal expression of lncRNAs might accelerate myopia progression through altering scleral ECM content. In addition, the α1 and α2 chains of type IV collagen, encoded by the Col4a1 and Col4a2 genes, respectively, are the two isoforms of type IV collagen and the main components of tissue basement membrane. α1(IV) and α2(IV) coexist in the basement membranes of multiple ocular tissues, including the RPE, choroidal vessels, and episcleral veins . According to previous studies, reduction in choroidal blood flow causes hypoxia in the sclera, which results in scleral ECM remodeling and subsequent formation of axial myopia . Therefore, we hypothesize that lncRNA upregulation may impair the basement membranes of choroid vessels by modulating the expression of Col4a1 and Col4a2, which may further alter the morphology and stability of choroidal vessels and lead to the initiation and progression of myopia. However, the target gene LOC106027164 of an upregulated lncRNA in the LIM group was annotated as a collagen α1(I) chain-like gene in the guinea pig. Type I collagen, generally composed of two chains of collagen type I α1 (COL1A1) and collagen type I α2 (COL1A2), is the main component of scleral ECM, and the synthesis and degradation of type I collagen are strongly associated with scleral remodeling under myopia . Moreover, the scleral Col1a1 gene promoter and exon 1 cytosine-phosphate-guanine (CpG) islands have been shown to be hypermethylated in a mouse model of diffuser lens-induced myopia . Therefore, the results of this work and that of others suggest that under LIM, type I collagen may subserve a molecular target regulated by epigenetic mechanisms, such as DNA methylation and lncRNA.
Moreover, the GO analysis in this study revealed that the target genes of the differentially expressed lncRNAs in the two myopic groups were involved in collagen turnover, ECM structural component degradation, kinase activity, growth, and development. The signal pathway enrichment analysis also showed that in the FDM and LIM groups, the lncRNA target genes were statistically significantly enriched in the ECM-receptor interaction signaling pathway. These results are consistent with previous reports, showing that the expression of sclera-related genes, such as matrix metallopeptidase-2 , gelatinase A, and tissue inhibitor of metalloproteinase 2 , is changed under myopic conditions.
Taken together, the bioinformatic analyses suggest that the ECM in the sclera, choroid, RPE, and choroidal vessels, in particular collagen, a major component of ECM, seems to be a similarity in the pathogenic mechanisms between FDM and LIM. However, compared to the normal controls, there were 372 and 247 differentially expressed lncRNAs in the FDM and LIM groups, respectively. Yet only 68 of these differentially expressed lncRNAs were shared by and exhibited similar trends in the two myopia groups. Further, the direct comparison of the FDM and LIM groups resulted in identification of 380 differentially expressed lncRNAs. These results implicate that different molecular mechanisms may underlie FDM and LIM despite the partial similarities the models share.
However, one limitation in this study is that the ocular posterior pole was not further dissected into the sclera, choroid, and retina, and lncRNA expression profiling was detected in the whole posterior pole. As shown in the morphology analysis, the thickness of the posterior pole, including the thicknesses of the sclera and the choroid, was dramatically reduced 6 weeks after the induction of myopia. Additionally, the amplitude of lncRNA expression was low compared to the RNA species we have been familiar with, such as mRNA and microRNA. Therefore, the ensemble of the ocular posterior pole was employed in an attempt to detect thorough expression profiling of low-amplitude lncRNAs. It is possible that this approach may mask the differential gene expression in individual tissues of the posterior pole. Nevertheless, the lncRNA expression profiling of the ocular posterior pole was a pool of clues from which major pathogenic or interventional targets can be selected. In future studies, it would be interesting to confirm the expression patterns of the selected targets in different tissues of the ocular posterior pole using microdissection and qPCR, as well as RNAscope in situ hybridization.
In conclusion, high through-put RNA-seq was performed in the ocular posterior poles of two guinea pig models of myopia, FDM and LIM. The lncRNAs were aberrantly expressed in the FDM and LIM groups compared to the healthy and contralateral controls. The functions of the lncRNA target genes may be related to myopia. The results, for the first time, suggest a link between lncRNAs and pathogenesis of myopia. These differentially expressed lncRNAs may provide a clue to the epigenetic pathogenesis of myopia and to interventional targets for this epidemic visual disorder.
This work was supported by The Key Project of Tianjin Municipal Science and Technology Commission (17JCZDJC35600); The High-level Innovative Talent Program for Distinguished Scholar (YDYYRCXMB201802); The project from National Natural Science Foundation of China (81970827); National and Municipal College Students' Innovation and Entrepreneurship Training Program (201910062021).