RNA sequencing analysis of long non-coding RNA expression in ocular posterior poles of guinea pig myopia models

Animals

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).

Guinea pig FDM and LIM models

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 [15]. 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 [16]. 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.

Ocular biometric parameters

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 [17], 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.

Sample collection

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.

Hematoxylin and eosin staining

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.

RNA isolation

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.

Library preparation and sequencing

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).

Data analysis for RNA sequencing

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.

Functional group analysis

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.

Quantitative real-time PCR

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 analysis

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.

Ocular biometric parameters

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.

H&E staining

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).

RNA-seq data quality

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.

Analysis of SNPs

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).

Screening of lncRNAs

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).

Expression profiling of lncRNAs in experimental groups

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).

LncRNA target gene analysis

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

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 enrichment analysis

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).

Verification of lncRNA expression

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.