Analysis of dynamic and widespread lncRNA and miRNA expression in fetal sheep skeletal muscle

Ethics statement

All experimental protocols and procedures were approved by the Institutional Animal Care and Use Committee of Lanzhou Institute of Husbandry and Pharmaceutical Science of the Chinese Academy of Agricultural Sciences (approval no. NKMYD201805, dated 18 October 2018).

Animal and tissue samples

A total of 15 Gansu Alpine fine wool sheep were used in this study. All Gansu Alpine fine wool sheep were raised in the experimental facilities of Gansu Provincial Sheep Breeding Technology Extension Station (Huangcheng, Gansu, China) under the same conditions with free access to food and water in natural lighting. Caesarean sections were performed on three pregnant ewes from each developmental stage to collect female foetuses at 60, 90, and 120 days of gestation, and six female lambs were collected 0 and 360 days after birth. All animals were slaughtered after being anesthetized with xylazine chlorhydrate. We made all efforts to minimize animal suffering, and the slaughter procedures were carried out in accordance with animal welfare procedures. After slaughter, we collected three gestational (E60, E90, and E120) and two postnatal stage (D0 and D360) LD muscle samples. The tissue samples were snap frozen in liquid nitrogen and stored at −80 °C until analysis.

Muscle staining

We prepared the LD samples for histological sectioning using Carter & Clarke’s (1957) method. The LD samples were first placed in tubes containing 4% paraformaldehyde solution. After fixation, the samples were embedded, cut into slices, baked, H&E strained, and mounted (Auber, 1952). We examined the sections using a digital trinocular camera microscope (BA400Digital, McAudi Industrial Co., Ltd., Xiamen, China), and used the image analysis software Motic Images Advanced 3.2 to take and import an image. We then selected the objective lens magnification (40 ×) and the unit of measurement (um), used the polyline tool to measure the data, and exported the measured raw data in the .xls format for sorting and analysis.

RNA isolation, library construction, sequencing, and data analysis

The total RNA from 15 muscle samples was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA, USA), according to the manufacturer’s instructions. We monitored RNA degradation and contamination using 1% agarose gels. RNA purity and concentration were detected using NanoDrop 2000 (Thermo Fisher Scientific Inc., Waltham, MA, USA), and were further measured using an Agilent 2100 bioanalyzer. Samples with an RNA integrity number (RIN) value greater than 8.0 were used for sequencing. RNA was digested by TruSeq stranded total RNA and a Ribo-Zero Gold Kit (Illumina, San Diego, CA, USA). After the total RNA was extracted, we removed the ribosome RNAs (rRNAs) to retain mRNAs and ncRNAs. The enriched mRNAs and ncRNAs were then broken down into short fragments using fragmentation buffer and turned into cDNA via reverse transcription and random primers. The second-strand cDNA was synthesized using DNA polymerase I, RNase H, dNTP (dUTP instead of dTTP), and buffer. Next, the cDNA fragments were purified using the Qiaquick PCR extraction kit, end-repaired, combined with poly A, and ligated to Illumina sequencing adapters. We used uracil-N-glycosylase (UNG) to digest second-strand cDNA that had been size-selected by agarose gel electrophoresis, amplified by PCR, and sequenced using an Illumina HiSeq TM 2500 from OE Biotechnology Corporation (Shanghai, China) (Zhou et al., 2016; Zhu et al., 2017). We pooled equal amounts of total RNA from muscle samples across different stages (i.e., E60, E90, E120, D0, and D360; n = 3) into one sample. To get high quality clean reads from the sequencing machines, we filtered out those that: (1) contained adapters, (2) had more than 10% unknown nucleotides (N), or (3) contained more than 50% low quality (Q-value ≤ 20) bases. We used the short reads alignment tool Bowtie 2 (Langmead & Salzberg, 2012) to map reads to the rRNA database, and then removed the rRNA mapped reads. The remaining reads were used in the transcriptome assembly and analysis. We then mapped the rRNA-removed reads from each sample to the reference genome using TopHat2 (Kim et al., 2013) and the following alignment parameters: (1) a maximum read mismatch of 2, (2) a 50 bp distance between mate-pair reads, and (3) a ± 80 bp error of distance between mate-pair reads. After aligning them with the reference genome, we re-aligned the unmapped reads (or very poorly mapped reads) using Bowtie2, and split the enriched unmapped reads into smaller groups to find potential splice sites. The section and section positions of these short segments were also predicted. We built a set of splice sites using initial unmapped reads from TopHat2 without relying on known gene annotation (Trapnell et al., 2010). The sequence alignments are not only useful for identifying expressed genes and their quantitative expressions, but also for identifying alternative splicing and new transcripts.

Transcript reconstruction was carried out using Cufflinks (Trapnell et al., 2012) which, together with TopHat2, allows biologists to identify new genes and new splice variants of known genes. We preferred to use the program’s reference annotation-based transcripts (RABT). Cufflinks constructed faux reads according to the reference gene to compensate for low coverage sequencing. During the last step of assembly, all reassembled fragments were aligned with reference genes and any similar fragments were removed. We used Cuff merge to combine transcripts from different replicas of a group into a comprehensive set of transcripts, and then to merge the transcripts from multiple groups into a comprehensive set of transcripts for further downstream differential expression analysis.

We quantified transcript abundances using RSEM (Li & Dewey, 2011). Transcript expression levels were normalized using the Fragments Per Kilobase of transcript per Million mapped reads (FPKM) method, and the differentially expressed transcripts of coding RNAs and lncRNAs were individually analysed. To identify differentially expressed transcripts across samples or groups, we used the edge R package (http://www.r-project.org/). We identified transcripts with a fold change ≥ 2 and a false discovery rate (FDR) <0.05 as significant DEGs. DEGs were then subjected to Gene Ontology (GO) function and KEGG pathway enrichment analysis. We performed gene expression pattern analysis to cluster genes of similar expression patterns from multiple samples (at least three from a specific time point, space, or treatment dose size order). To examine the DEG expression patterns, we normalized the expression data from each sample (in the order of treatment) to 0, log2(v1/v0), log2(v2/v0), and then clustered them using Short Time-series Expression Miner software (STEM) Ernst & BarJoseph, 2006. The clustered profiles with p-values ≤ 0.05 were considered significant, and the DEGs from each profile underwent GO and KEGG pathway enrichment analysis. Using the hypothesis test for p-value calculation and FDR correction (Saldanha, 2004), we defined GO terms or pathways with Q values ≤ 0.05 as significant enriched GO terms or pathways.

Real-time quantitative RT-PCR

The 15 muscle samples were stored at −80 °C prior to RNA extraction. We extracted the total RNA from the 15 muscle samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA purity and concentration were determined using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific), and RNA integrity was evaluated using agarose gel electrophoresis staining with ethidium bromide.

Quantification was performed using a two-step reaction process: reverse transcription (RT) and PCR. Each RT reaction consisted of two steps. The first step was to combine 0.5 µg of RNA, 2 µl of 4 ×g DNA wiper Mix, and 8 µl of nuclease-free H2O. Reactions were performed in a GeneAmp^® PCR System 9700 (Applied Biosystems, Foster City, CA, USA) for 2 min at 42 °C. The second step was to add 2µl of 5 × HiScript II Q RT SuperMix IIa. Reactions were performed in a GeneAmp^® PCR System 9700 (Applied Biosystems) for 15 min at 50 °C, then for 5 s at 85 °C. The 10 µl RT reaction mix was then diluted ×10 in nuclease-free water and maintained at −20 °C. Real-time PCR was performed using a LightCycler^® 480 II Real-time PCR Instrument (Roche, Basel, Switzerland) and 10 µl of a PCR reaction mixture that included 1 µl of cDNA, 5 µl of 2 ×ChamQ SYBR qPCR Master Mix, 0.2 µl of forward primer, 0.2 µl of reverse primer, and 3.6 µl of nuclease-free water. Reactions were incubated in a 384-well optical plate (Roche) at 95 °C for 30 s, followed by 40 cycles of 10 s at 95 °C, and 30 s at 60 °C. Each sample was run in triplicate for analysis. After the PCR cycles, we performed melting curve analysis to validate the specific generation of the expected PCR product. The primer sequences were designed and synthesized by Generay Biotech (Shanghai, China) using the lncRNA and miRNA sequences obtained from the NCBI database (Tables S15 and S16 ; (Bai et al., 2017). We estimated the PCR efficiency of each gene using standard curve calculation of the four cDNA serial dilution points. Cycle threshold (Ct) values were transformed to quantities using the comparative Ct method described by Chen et al. (2017). We carried out lncRNA and mRNA data normalization using the GAPDH reference gene, and miRNA data normalization using the U6 reference gene. The lncRNA and miRNA expression levels were normalized (using GAPDH and U6) and estimated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Statistical analyses

After comparing the clean reads to the transcription template using Bowtie2 software (Langmead & Salzberg, 2012), we quantified the transcripts using eXpress software (Roberts & Pachter, 2013) to obtain FPKM values and mRNA and lncRNA counts. For samples without biological replicates, we calculated the p-value using the Audic_Claverie formula (Tiňo, 2009). miRNAs with p-values <  0.05 and TMP difference multiples >  2 were screened.

Developmental changes in LD muscle across different stages

In this study, we stained the 12 LD in histological sections to analyze the area and circumference of muscle fibers at the various development stages from embryonic to postnatal. Our analysis showed that the area and circumferences of the muscle fibers did not increase over time, but it did decrease during the E90 and E120 periods (Fig. 1A). However, the tightness and evenness of muscle fibers significantly improved as the fetus continued to develop (Fig. 1B). These results suggest that embryonic skeletal muscle growth can be characterized by time.

LncRNA and small RNA sequencing

To comprehensively analyze sheep lncRNAs and miRNAs across different developmental stages, we constructed five cDNA libraries and five small RNA libraries (E60, E90, E120, D0, and D360) from three LD samples taken at 60, 90, and 120 days of gestation (E60, E90, and E120), and two from the postnatal (D0, D360) developmental stage. During lncRNA sequencing, we generated a total of 504,449,178 raw reads using all five libraries. After discarding adaptor sequences and low-quality reads, we obtained 477,031,266 clean reads (Table S1). Additionally, about 91.87% (91.70% - 92.48%) of the clean reads from each library were mapped to the sheep reference genome (Table S1), confirming the reliability of the sequencing data. During small RNA sequencing, we generated a total of 84,426,880 raw reads, and we obtained 10,623,459 - 26,521,572 clean reads (97%) ranging in size from 18 to 30 nt (Table S1). About 8,513,544 - 24,586,094 reads in the LW samples were perfectly mapped to the reference genome (GCF_000298735.2_Oar_v4.0), amounting to 72.2% - 92.7% of the clean reads (Table S2). Our sequencing results showed that most of the small RNA sequences were between 21 and 23 nt in length, which was consistent with the length distribution of the Dicer products and indicated that the sequencing results were good quality and could be used for follow-up analysis.

Differential expression analysis of sheep mRNAs and lncRNAs

To investigate the key mRNAs and lncRNAs involved in regulating sheep skeletal muscle development, we used RNA-seq datasets from five time points (three gestation stages and two postnatal stages) to characterize their time-specific expression patterns. When comparing the gene expression levels across the five developmental stages, we found 693 (352 upregulated) DEGs between E60 and E90 (Table S3), 799 (278 upregulated) DEGs between E90 and E120 (Table S4), 929 (672 upregulated) DEGs between E120 and D0 (Table S5), and 815 (257 upregulated) DEGs between D0 and D360 (Figs. 2A & 2B, Table S6). When analyzing these DEGs, we found that 478, 535, 633, and 580 genes were uniquely expressed in one of the two samples in E60 vs. E90, E90 vs. E120, E120 vs. D0, and D0 vs. D360, respectively (Fig. 2A). We detected two of these DEGs, FBN2 (gene ID: 101104991) and HR (gene ID: 443241), across all four comparisons. We also analyzed the differentially expressed lncRNAs (DE-lncRNAs) between E60 vs. E90, E90 vs. E120, E120 vs. D0, and D0 vs. D360, and detected 206, 239, 148, and 101 DE-lncRNAs, respectively (Fig. 2C & 2D, Tables S7–S10). We detected six of these DE-lncRNAs, TCONS_00009825, TCONS_00028121, TCONS_00122104, TCONS_00146167, TCONS_00189359, and TCONS_00195408, across all four comparisons. In the embryonic stage, the following gene expressions were significant: TCONS_00054028, TCONS_00023110, and TCONS_00150168 (Fig. 2E). Similarly, we found post-birth differential expression of other related genes, including TCONS_00016124, TCONS_00127579, TCONS_00009472, and TCONS_00014919 (Fig. 2F). Five differentially expressed lncRNAs were selected for qRT-PCR analysis (Table 1), and their differential expression was consistent with the RNA-Seq results. In summary, we found temporal-specific expression of the detected lincRNAs during the gestational and postnatal stages.

GO and KEGG pathway analysis

We used KEGG pathway analysis of DE-lncRNA target genes to identify the pathways that were enriched in DE-lncRNA target genes. When comparing E60 vs. E90, E90 vs. E120, E120 vs. D0, and D0 vs. D360, we found that in the most significantly enriched pathways, DE-lncRNA target genes participated in signal transduction, the endocrine system, the nervous system, cell growth and death folding, sorting and degradation, the immune system, cellular community eukaryotes, translation, amino acid metabolism, and the carbohydrate metabolism pathways (Fig. 3). The top 20 significantly enriched KEGG analyses for each comparison’s DE-lncRNAs are shown in Fig. S1. GO enrichment analysis of DE-lncRNA targeted genes also delivered a large number of significant annotations, but in order to examine temporal changes in skeletal muscle development, we provided the detailed results from four adjacent DE-lncRNA comparisons. In the E60 vs. E90 comparison, T=the top 30 GO terms that were significantly related to genes targeted by total lncRNAs included male meiotic nuclear division, mitotic spindle assembly checkpoint, kinetochore, cell, protein kinase activity, and microtubule plus −end binding (Fig. S2). In the E90 vs. E120 comparison, the most significantly enriched GO terms for the total lncRNAs were associated with muscle development, and included signal peptide processing, activation of MAPKK activity, L −type voltage −gated calcium channel complex, voltage −gated calcium channel complex, histone acetyltransferase complex, receptor signaling complex scaffold activity, and Rho GTPase binding (Fig. S3). In the E120 vs. D0 comparison, the top 30 processes for down and upregulated lncRNAs were associated with male meiotic nuclear division, cellular response to leukemia inhibitory factor, intercellular bridge, and integral component of endoplasmic reticulum membrane (Fig. S4). In the D0 vs. D360 comparison, we found associations with male meiotic nuclear division, cellular response to leukemia inhibitory factor, microtubule organizing center, and integral component of endoplasmic reticulum membrane (Fig. S5). Across the four comparisons, we observed significant changes in skeletal muscle development during the prenatal and neonatal stages.

LncRNA-gene interaction network construction

We extracted the candidate sequences by screening out lncRNAs and genes that were not on the same chromosome as the candidate targets. We used the RNA interaction software RIsearch (Wenzel, Akbaşli & Gorodkin, 2012) to predict the binding of candidate lncRNAs and genes at the nucleic acid level, with the number of bases directly interacting with each other according to the two nucleic acid molecules was no less than 10, and the free energy of base binding was no more than -50. The possible regulatory interaction networks between lncRNAs and their target genes (mRNAs) were constructed. In the E60 vs. E90 comparison, the lncRNA-gene interaction network was comprised of 15 lncRNAs and 66 protein-coding genes with close networks; among these, TCONS_00196403 and TCONS_00196407 constructed complex network relationships with targeted coding genes (Fig. 4A). In the E90 vs. E120 comparison, the lncRNA-gene interaction network was made up of complex network nodes and lncRNA-gene connections between 34 lncRNAs and 79 protein-coding genes. TCONS_00029967 and TCONS_00088616 had a strong correlation with protein-coding genes (Fig. 4B). In the E120 vs. D0 comparison, the lncRNA-gene interaction network was comprised of complex network nodes and lncRNA-gene connections between 25 lncRNAs and 95 protein-coding genes. TCONS_00181895, TCONS_00149972, and TCONS_00098734 had a strong correlation with protein-coding genes (Fig. 4C). In the D0 vs. D360 comparison, the lncRNA-gene interaction network was comprised of complex network nodes and lncRNA-gene connections between 11 lncRNAs and 86 protein-coding genes. TCONS_00196411, TCONS_00157638, and TCONS_00126396 had a strong correlation with protein-coding genes (Fig. 4D).

Temporal-specific expression of the detected sheep miRNAs

The miRNA expression level comparisons across the five developmental stages revealed that there were 328 (265 upregulated) DEGs between E60 and E90 (Table S11), 302 (243 upregulated) DEGs between E90 and E120 (Table S12), 295 (166 upregulated) DEGs between E120 and D0 (Table S13), and 184 (110 upregulated) DEGs between D0 and D360 (Fig. 5A, Table S14). We found that the miRNA expression levels during the embryonic and postnatal stages differed significantly. During the embryonic stage, the differentially expressed miRNAs were miR-134, miR-2387, miR-105, miR-3957, miR-493, miR-541, miR-767, miR-432, and miR-433 (Fig. 5B). We also conducted a differential analysis of the miRNAs based on their expression levels. The miRNAs with expression levels that changed more than two times and with P values <0.05 were deleted. The results showed that some miRNAs were differentially expressed after birth, specifically miR-22, miR-365, miR-3556, miR-29, miR-193, miR-150, miR-1, and miR-133 (Fig. 5C). Five differentially expressed miRNAs were selected for qRT-PCR analysis (Table 1), and their differential expression was consistent with the RNA-Seq results. We generated a heatmap based on the expression patterns at all time points and across all known expressed miRNAs (Fig. 5D). The significantly enriched processes and observed pathways were consistent with the temporal changes in skeletal muscle development across all sampled stages.

miRNA target gene prediction and functional analysis

Most animal miRNAs were not completely complementary to their target mRNA, mainly in the 3′ non-coding region (3′ UTR) of the target mRNA, and their mechanism of action is translation inhibition. Previous studies also found that animal miRNAs could target the 5′ end of the mRNA as well as the coding region (Aleksandra et al., 2013). In this study, we used the miRanda algorithm (John et al., 2004) to predict miRNA target genes. We detected 216, 143, 198, and 204 differentially expressed miRNAs and 33,106, 26,431, 29,049 and 27,554 target genes in the E60 vs. E90, E90 vs. E120, E120 vs. D0, and D0 vs. D360 comparisons, respectively. The high expression of novel 129_mature in E60 vs. E90 indicates that the expression product of the target AHNAK2 gene is related to AHNAK nucleoprotein 2 and transcript variant X1, and is involved in the nucleus, cytoplasm, cytosol, plasma membrane, Z disc, T-tubule, cytoplasmic vesicle membrane, and sarcolemma pathways. The high expression of novel 410_mature in E90 vs. E120 suggests that the expression product of the target AHNAK2 gene is mainly related to the S-phase response (cyclin-related). The high expression of novel 360_mature in E120 vs. D0 and D0 vs. D360 shows that the expression product of the target PASD1 gene is related to the PAS domain containing 1 and transcript variant X1, and is involved in the transcription coactivator binding, nucleus, nuclear speck, negative regulation of circadian rhythm, negative regulation of transcription, DNA-templated, rhythmic process, and Cry-Per complex pathways. We performed GO enrichment analysis at GO Level 2 for all target genes and the genes of differentially expressed miRNAs. The biological processes involved included cellular process, biological regulation, metabolic process, response to stimulus, and regulation of biological process (Figs. S6 & S7).