Genome-wide identification and characterization of long non-coding RNAs expressed during sheep fetal and postnatal hair follicle development

Long non-coding RNAs (lncRNAs), >200 nt in length, are transcribed from mammalian genomes. They play important regulatory roles in various biological processes; However, the function and expression profile of lncRNAs involved in the development of hair follicles in the fetus, have been relatively under-explored area. To investigate the specific role of lncRNAs and mRNAs that regulate hair follicle development, we herein performed a comprehensive study on the lncRNA and mRNA expression profiles of sheep at multiple embryonic days (E65, E85, E105, and E135) and six lambs aged one week (D7) and one month (D30) using RNA-seq technology. The number of genes (471 lncRNAs and 12,812 mRNAs) differentially expressed and potential targets of differentially expressed lncRNAs were predicted. Differentially expressed lncRNAs were grouped into 10 clusters based on their expression pattern by K-means clustering. Moreover, Gene Ontology (GO) annotation and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses showed that some differentially expressed mRNAs, such as DKK1, DSG4, FOXE1, Hoxc13, SFRP1, SFRP2, and Wnt10A overlapped with lncRNAs targets, and enriched in important hair follicle developmental pathways, including Wnt, TNF, and MAPK signaling pathways. In addition, 9 differentially expressed lncRNAs and 4 differentially expressed mRNAs were validated using quantitative real-time PCR (qRT-PCR). This study helps enrich the Ovis lncRNA databases and provides a comprehensive lncRNA transcriptome profile of fetal and postnatal skin of sheep. Additionally, it provides a foundation for further experiments on the role of lncRNAs in the regulation of hair growth in sheep.


Results
Overview of sequencing data. In order to identify lncRNAs and mRNAs expressed during the development of HFs in sheep, we constructed 18 cDNA libraries. E65-1, E65-2, and E65-3 as the replicate libraries for the E65 experiment group; E85-1, E85-2, and E85-3 as the replicate libraries for the E85 experiment group; E105-1, E105-2, and E105-3 as the replicate libraries for the E105 experiment group; E135-1, E135-2, and E135-3 as the replicate libraries for the E135 experiment group; D7-1, D7-2, and D7-3 as the replicate libraries for the D7 experiment group; D30-1, D30-2, and D30-3 as the replicate libraries for the D30 experiment group using Subo Merino sheep back-skin samples. Three biological replicates were used for each stage. A total of 1,847,973,012 raw reads were obtained in all 18 libraries. After removing the low-quality reads and adapter fragments, we obtained 1,726,188,919 clean reads, which corresponded to 93.41% of the raw reads from all the libraries. The rRNA content was between 0.03-0.39%; Q20% in all the libraries were above 95%, indicating that clean reads were in good quality. 81-88% of the clean reads in all the libraries were successfully mapped to the reference genome; the average mapping percentage was 81.5%. Out of successfully mapped reads, 82.20% uniquely mapped reads were used in transcript construction. In addition, the sequences of novel transcripts from sheep skin tissue were exist on all chromosomes, except for the mitochondrial genome; the sequences of new transcripts were prevalent on sheep chromosome 1, 2, and 3, which are also the three longest chromosomes.
Identified lncRNAs. After filtering and coding potential evaluation steps, lncRNA transcripts of libraries resulting in 10,193 (1,540 known and 8,653 novel) lncRNAs were considered expressed in skin tissues, respectively, and used in further analyses. In addition, a total of 27,054 mRNAs were identified. To gain insight into the functional role of lncRNAs that are expressed during sheep HF development, we analyzed their genomic context based on their orientation to local protein-coding genes. We identified 95 exonic sense, 1213 intronic sense, 1467 exonic antisense, 988 intronic antisense, 557 bidirectional, and 5873 intergenic lncRNAs (Fig. 1). The exon number, transcript length, and expression level of lncRNAs and mRNAs were calculated and graphed, as shown in Fig. 2a-c. The transcript length (Fig. 2a) and expression level (Fig. 2b) of lncRNAs were significantly shorter than those of mRNAs. 79% of lncRNAs comprised two exons, whereas mRNAs contained a broad range of exon numbers from two to nineteen (Fig. 2c). All in all, lncRNAs were found to be shorter in length, lower expression level, and fewer in exon number than protein-coding genes in skin (P < 0.05)-findings that are in agreement with those of previous studies 9,39 . Differentially expressed lncRNAs. A comparison of lncRNA and mRNAs expression levels were estimated by FPKM (fragments per kilo-base of exon per million fragments mapped) using Cuffdiff. RNA sequencing analyses detected 471 lncRNAs (174 were up-regulated and 297 were down-regulated) (Fig. 3) and 12,812 mRNAs (7,269 were up-regulated and 5,543 were down-regulated) that were differentially expressed in skin during HF development. Additionally, 439 differentially expressed novel lncRNAs were obtained. The number of up-regulated genes of lncRNA was greater than that of down-regulated genes at the E105-VS-E85 comparison group, while the number of down-regulated genes was greater than that of the up-regulated genes in skin at other comparison groups. The expression patterns of differentially expressed lncRNAs were measured by systematic www.nature.com/scientificreports www.nature.com/scientificreports/ cluster analysis, to explore the similarities and to compare the relationships between the different libraries. The replicates for each developmental stage clustered together, and E65, E85, and E105 formed one group and E135, D7, and D30 formed another group Data are expressed as FPKM. Red: relatively high expression; Green: relatively low expression. The bar code on the right represents the color scale of the log2 values (Fig. 4). To further analyze the interactions among the differentially expressed lncRNAs, we constructed a Venn diagram using the 35,55,46,43,267, and 27 lncRNAs that were differentially expressed in E85-VS-E65, E105-VS-E85, E135-VS-E105, D7-VS-E135, D30-VS-E65, and D30-VS-D7, respectively (Fig. 5). We did not detect any co-expressed differentially expressed lncRNA in all six developmental stages, but we identified 7,7,8,8,170, and 2, respectively, stage-specific differentially expressed lncRNAs. We then compared the expression level of these differentially expressed lncRNAs between the samples for every two HF stages. Even though, the level of expression of most of these lncRNAs was equivalent, there were differences in the level of expression between the two stages for some lncRNAs.

Target gene prediction of lncRNAs.
To understand the function of the differentially expressed lncRNAs in different HF developmental stages, we predicted the target genes of cis-and trans-acting lncRNAs. A total of 3,685 potential target genes were identified for 5,852 cis-acting lncRNAs, while 1,722 potential target genes were identified for 4,117 trans-acting lncRNAs. These 3,685 and 1,722 candidate target genes were subjected to GO and pathway analyses. GO enrichment and KEGG pathway analysis of lncRNA targets and mRNAs. Functional analysis showed that cis and trans target genes of lncRNAs were enriched in 5,509 GO terms. GO terms of lncRNA potential targets were highly enriched for several processes, such as organic substance metabolic, regulation of stimulus, and cellular metabolic process. Importantly, 13, 5, 12, and 7 predicted lncRNA targets overlapped with the differentially expressed mRNAs, and predicted to function in HF development and morphogenesis, epidermis development and morphogenesis, skin development and morphogenesis, and cell differentiation and migration, respectively (see Table 1).
Based on the KEGG pathway database, pathway analysis was performed to predict the significantly enriched pathways involving putative target genes of differentially expressed lncRNAs. Pathway analysis showed that 268 pathways involving the target genes of differentially expressed lncRNAs were enriched; some of these pathways are known to play an important role in the regulation of HF development, including Wnt, TNF, MAPK, and Hedgehog which shared 7, 1, 2, and 1 genes with lncRNA targets, respectively (see Table 1). To characterize the GO terms and pathways, we listed the top 30 significantly enriched GO terms and pathways (P < 0.05) in six comparative groups.
It is worth mentioning that SFRP2 40 , DKK1 41 , and Hoxc13 42 obtained from the GO enrichment and KEGG pathway analysis of lncRNA targets and mRNAs have been proved to be key genes related to HF development, and in this experiment, these three genes are enriched in GO items and Wnt signaling pathway related to HF development. SFRP2, DKK1, and Hoxc13 are the predicted target genes of the differentially expressed lncRNAs such as TCONS_00202353, TCONS_00453233, and TCONS_00428783, which may regulate the development of HF (Fig. 6).

K-means cluster analysis of lncRNAs.
To identify the expression pattern of differentially expressed lncR-NAs during six HF developmental stages, we performed k-means clustering of all differential lncRNAs. LncRNAs were classified into 10 clusters based on their group mean and changes with respect to other compartments. Each box is a differential lncRNA cluster, each line in the box represents a expression trend of differential lncRNAs, a black bold lines is the average trend of all differential lncRNAs in the box (Fig. 7). These clusters can be further assigned to four main groups of different dynamic patterns. The first pattern group, comprising clusters 2 and 4 included 29 and 64 differentially expressed lncRNAs, represents lncRNAs whose expression level gradually www.nature.com/scientificreports www.nature.com/scientificreports/ increased from that of lncRNAs in Group 2. The second pattern group, comprising of clusters 1, 3, 6, and 7 included 18, 31, 58, and 112 differentially expressed lncRNAs respectively, represents lncRNAs that showed an overall decreasing trend. The third pattern group, consisting of clusters 5, 8, and 10 included 12, 12, and 11 differentially expressed lncRNAs respectively, represents lncRNAs that were up-and down-regulated continuously in all stages. The expression level of 10 differentially expressed lncRNAs in cluster 9, which is the last pattern group, remained stable at the entire stages.
Validation of the sequencing data by qRT-PCR. To validate the sequencing data, we investigated the expression pattern of 9 randomly selected lncRNAs and 4 mRNAs, which were identified to be differentially www.nature.com/scientificreports www.nature.com/scientificreports/ expressed during two HF developmental stages (E85 and E105; because these two groups are the important stages for the primary HF development and secondary HF formation) 43 using qRT-PCR. As shown in Fig. 8, all of the selected lncRNAs and mRNAs were detected in skin tissues, and of these, 12 showed differential expression at two different time points. In addition, TCONS_00280360 did not show differential expression, and this result was inconsistent with that of RNA-seq. Collectively, the sequencing results were reliable.

Discussion
Several lncRNAs relevant to skin biology have been reported, such as ANCR and SPRY4-IT1 44 , the detailed functional roles of lncRNAs in hair and wool follicle development are based on limited reports and, thus, largely unknown. A recent study detected 15 differentially expressed lncRNAs during secondary wool follicle induction in super-fine wool sheep, which possesses no primary wool follicles due to selective breeding 12 . In recent years, increasing evidence has demonstrated that lncRNAs play important roles in the development of various tissues, such as brain 45 , liver 46,47 , heart 48 , testis [49][50][51] , ovary 52 , kidney 53 , uterus 54 , muscle [55][56][57] , mammary gland 58 , and inner ear 59 . However, very few studies have been conducted on the potential role of lncRNAs in mammalian skin HF development 3,36 . Regarding sheep  www.nature.com/scientificreports www.nature.com/scientificreports/ skin tissues, investigations into the molecular mechanisms of HF development have focused more on miRNAs 60-62 than on lncRNAs, and therefore, the role of lncRNAs in sheep HF development has not been fully illustrated.
In the present study, we systematically identified the lncRNAs and mRNAs involved in sheep (Ovis aries) skin during different HF development across four fetal and two postnatal stages, using RNA-sequencing technology. We subsequently characterized putative lncRNAs to elucidate their diverse features in order to provide an insight into the uncharacterized sheep genome regions and their relationship with sheep HF development. To the best of our knowledge, this is the first study to systematically identify lncRNAs from RNA-seq data of sheep HF development during fetal and postnatal stages. In addition, we identified a total of 10,193 and 27,054 expressed lncRNAs and mRNAs in different sheep HF developmental stages; in which, many lncRNAs that were not reported in previous data were identified.
Unlike protein-coding genes or miRNAs, the function of lncRNAs could not be inferred from previous sequencing or structural studies performed by us and others 63,64 . Therefore, we analyzed the genomic context of identified lncRNAs and classified them based on their genomic relationship with protein-coding genes as exonic sense, intronic sense, exonic antisense, intronic antisense, bidirectional, and intergenic. We found that all six categories of lncRNAs could be examined. LcRNAs identified in sheep skin shared many of the features of other mammalian counterparts: shorter length, lower expression level, and fewer exon number than those of protein-coding genes 37,62 . These clues further indicate that mammalian lncRNAs share some specific characteristics, and provide a reference for studying lncRNAs in other species.
A total of 471 differentially expressed lncRNA genes were identified in fetal and postnatal skin tissues. We also performed clustering analysis on 18 samples showing differential expression of lncRNAs at different HF developmental stages. The results indicated that the early stages of HF development (E65, E85, and E105) are clustered together, while the later stages of HF development (E135, D7, and D30) are clustered together. The above findings suggest that the changes in the expression of these lncRNAs in early stages of HF development are much more dramatic than those in later stages. To further survey the interaction of differentially expressed lncRNAs genes, we constructed a Venn diagram using lncRNAs among all the developmental stages. As a result, 7,7,8,8,170, and 2 stage-specific differentially expressed lncRNA genes were identified out of the 35, 55, 46, 43, 267, and 27 lncRNA genes, respectively. These results suggested that they may have an important role in HF development.
Most evidence suggests that the expression of lncRNA-coding genes can regulate and have a high correlation with the expression of neighboring mRNA-coding genes 63,64 . Based on this, we searched for coding genes 10-kb upstream and downstream of lncRNA-coding genes as target genes to predict the function of the cis-acting lncR-NAs. We also searched for mRNA-coding genes using co-expression prediction to predict function of trans-acting lncRNAs. Consequently, we found that many target genes of differentially expressed lncRNAs were also differentially expressed. This suggested that lncRNAs may function through cis/trans-regulation of target genes, and this may play a critical role in the development of sheep skin HFs.
We identified 471 differentially expressed lncRNAs genes and predicted their potential targets which partially overlapped with differentially expressed mRNAs and were enriched in HF developmental GO terms. Pathway analysis showed that 268 pathways involving the target genes of differentially expressed lncRNAs were enriched; some of these pathways are known to play an important role in the regulation of HF development, such as WNT 65 , TNF 66 , and MAPK signaling pathways 67 (Table 1). Interestingly, we found potential target genes, including DKK1, Hoxc13, DSG4, Wnt10A, FOXE1, SFRP1, and SFRP2 were overlapped with mRNAs, and enriched in the HF developmental pathways mentioned above, and were found to play important roles in HF development and morphogenesis [68][69][70][71][72][73][74] . Recently, many studies claimed that lncRNAs could promote or inhibit the expression of target genes by elevating or reducing the activity of pathways [75][76][77] . It was inferred that the above-mentioned lncRNAs, and their target genes could regulate the HF developmental process via regulating various pathways. www.nature.com/scientificreports www.nature.com/scientificreports/ K-means cluster analysis indicated that differentially expressed lncRNAs in the first pattern group had a positive correlation, and play important roles in the development of HF, while the differentially expressed lncRNAs in the second pattern group had a negative correlation, and lack of involvement in HF development. In addition, differentially expressed lncRNAs in the third pattern group had a complex regulatory role on the development of HF. Furthermore, differentially expressed lncRNAs in the last pattern group might be involved in the HF development at the postnatal stage. To validate the sequencing data, we randomly selected 9 different lncRNAs and 4 mRNAs for qRT-PCR analyses. Except TCONS_00280360, the qRT-PCR results of other selected lncRNAs and mRNAs are consistent with the sequencing data. RNA-seq is used for large-scale screening, which reflects the overall trend of gene expression in the sample, but does not guarantee that the trend of each gene is consistent with qRT-PCR. Generally speaking, the qRT-PCR result is more accurate. RNA-Seq and qRT-PCR have different experimental conditions and may get different results. Collectively, the sequencing results were reliable.
In conclusion, we systematically identified the lncRNAs and mRNAs involved in sheep (Ovis aries) skin during different HF development across four fetal and two postnatal stages using RNA-sequencing technology. A total of 471 differentially expressed lncRNAs and 12,812 differentially expressed mRNAs were identified from all the lncRNA libraries. GO and KEGG function enrichment analysis indicated that some target genes of these lncR-NAs as DKK1, Hoxc13, DSG4, Wnt10A, FOXE1, SFRP1, and SFRP2 were overlapped with mRNAs, and involved in HF developmental GO terms and pathways. The dynamic expression profile and analysis of certain features of these genes indicated that lncRNAs may play important roles in HF development of sheep. These results will serve as a very useful resource for genetic information that may improve our understanding of the mechanism of action lncRNAs involved in skin HF development.

Materials and Methods
Ethics statement. The methods were performed in accordance with the guidelines of the good experimental practices adopted by the Institute of Animal Husbandry. All experimental protocols were approved by the Institute of Animal Husbandry of the Xinjiang Academy of Animal Science. www.nature.com/scientificreports www.nature.com/scientificreports/ Animal selection and skin tissue preparation. The Subo Merino sheep, a breed of sheep in China, is famous for its excellent quality and yield of wool, and high survival rate. The tested individuals were selected from the Subo Merino sheep herd located at Xinjiang Kechuang animal husbandry breeding center, Urumqi county (Latitude 43°01′08″-44°06′11″N, Longitude 86°37′56″-88°58′22″E), Xinjiang Province, China. 40 Healthy Subo Merino ewes (2~3 years old), which were selected from the same flock and reared using the same management practices, were artificially inseminated with fresh sperm from a single ram. The day of insemination was designated as embryonic day 0 (E0). Three embryos and three lambs (three samples at each stage) were collected. Skin samples (~2 cm 2 for each individual) were collected from the right mid-side of embryos at four different embryonic days (E65, E85, E105, and E135) and from lambs at two postnatal days (D7 and D30). The samples were cut into small pieces, quickly placed into liquid nitrogen, and subsequently stored at −80 °C for RNA-seq and qRT-PCR analyses.
Total RNA isolation, library construction, and sequencing. Total RNA was isolated from the tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA quality was verified using a 2100 Bioanalyzer RNA Nano Chip (Agilent technologies, Santa Clara, CA, USA), and the total RNA content was measured using Quality control and transcriptome assembly. Sequencing-received raw image data were transformed by base calling into sequence data, which was called raw data. The raw data were subjected to quality control using FastQC (v0.11.4) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Raw reads of all 18 samples were pre-processed using in-house Perl scripts developed by the Novogene Bioinformatics Institute (Beijing, China). In this step, clean data were obtained after filtering out reads with adapters and poly-N > 10%, and low quality reads (>50% of bases whose Phred quality scores were <5%) from the raw data. At the same time, the Q20 (the 1% base error rate) and GC (the percent of guanine and cytosine) content of the clean data were calculated.
All subsequent analyses were based on the high-quality data. The clean reads were then mapped to the sheep reference genome (ensembl Oar v3.1) using TopHat (v2.0.9) 10 . Reference genome and gene model annotation files were downloaded from Ensembl browser (ftp://ftp.ensembl.org/pub/release-85/fasta/Ovis_aries/dna/Ovis_aries. Oar_v3.1.dna.toplevel.fa). The mapped reads from each sample were assembled by both Scripture (beta2) 78 and Cufflinks (v2.2.1) 9 in a reference-based approach. Scripture was run with default parameters. Cufflinks was run with 'min-frags-per-transfrag = 0' and '-library-type fr-firststrand' , and other parameters were set as default. Read counts for different gene and other known transcripts (misc_RNA, pseudogene, rRNA, tRNA, and others) were extracted by HTSeq (v0.6.1) 79 with default parameters and allowing one mismatch. To monitor mapping events on both strands, both sense and antisense sequences were included 80 .

Identification of lncRNAs.
Putative lncRNAs were defined as novel transcripts if they satisfied the following criteria: (1) Transcripts with single-exon and less than 200 bp in length were removed first; (2) All transcripts encoding complete open reading frames (ORFs) of more than 300 bp were excluded from our dataset; (3) Transcripts that did not belong to any other type of RNAs (rRNA, tRNA, snRNA, snoRNA, pre-miRNA, and those encoded by known pseudogenes) were excluded; (4) Among the different classes of class_code (http:// cufflinks.cbcb.umd.edu/manual.html#class_codes), only transcripts annotated as "u", "i", and "x" were retained, which represent potential novel intergenic, intronic and anti-sense lncRNAs, respectively; (5) In order to identify  81 , a database with a focus on the domestic-animal lncR-NAs, with the settings of identity = 100%, mismatch = 0, E-value; (6) Three coding potential analyses software, Coding-Non-Coding-Index (CNCI) (v2) 82 (score < 0), Coding Potential Calculator (CPC) (0.9-r2) 83 (score < 0), and Pfam-scan (v30.0) 84 (E-value < 0.001) were used to discover the candidate lncRNAs. Transcripts predicted with coding potential by any of the three tools earlier described were filtered out, and those without coding potential were retained. Then, we selected those shared by three tools as the final candidate lncRNAs and use for the further analysis. For mRNA identification, uniquely mapped and properly paired reads were used in transcript construction with Cufflinks (v2.2.1) 85 for mRNA identification. Constructed transcripts were compared with Ensembl Ovis aries gene annotation to identify expressed mRNAs using Cuffcompare 85 . The detailed workflow of novel lncRNA prediction is presented in Fig. 9.   Figure 7. K-means clustering analysis of differentially expressed lncRNAs among six comparison groups. Each box is a differential lncRNA cluster, each line in the box represents a expression trend of differential lncRNAs, a black bold lines is the average trend of all differential lncRNAs in the box.   www.nature.com/scientificreports www.nature.com/scientificreports/ Differential expression analysis and hierarchical clustering of lncRNAs. Cuffdiff (v2.1.1) was used to calculate fragments per kb for a million reads (FPKM) of both lncRNAs and mRNAs in each sample 10 , based on the TopHat BAM files and the reference GTF file. These counts were used in analyses of differential gene expression using the edgeR v3. 16.0 software 86 . LncRNAs and coding genes were differentially expressed between any two comparison groups of samples which representing a specific developmental stages were identified using edgeR 86 . For biological replicates, lncRNAs and mRNAs with a P-adjust ≤ 0.05 and an absolute value of the log2(fold change) ≥2 were assigned as being differentially expressed between different comparison groups 13 . The expression patterns of differentially expressed lncRNAs were measured by systematic cluster analysis, to explore the similarities and to compare the relationships between the different libraries. Scatter plots were used to demonstrate lncRNAs showing differential expression between every two HF developmental stages. Furthermore, the differentially expressed lncRNAs were subjected to K-means clustering using the Euclidean distance method associated with complete linkage on the BMKCloud platform (https://www.biocloud.net/). lncRNA genomic context analysis. We determined the genomic context of lncRNAs in relation to protein-coding genes according to a protocol that was updated for this study 87 . Briefly, intronic sense and antisense lncRNAs were defined, and the sequence of each corresponding transcript was mapped to the positive and opposite strand, respectively, of an intron of a protein-coding gene. Exonic sense and antisense lncRNAs were defined, and the sequence of each corresponding transcript was mapped to the positive and opposite strand, respectively, of a RefSeq-annotated exon [includes 5′-untranslated region (UTR), coding exon, and 3′-UTR]. Bidirectional lncRNAs were defined, and the sequence of each corresponding transcript was mapped head-head to protein-coding gene at a distance of <1000 bp. Intergenic lncRNAs were defined, and the sequence of each corresponding transcript was mapped within an intergenic region, without the presence of any overlapping or bidirectionally coded sequences for transcripts nearby.
Target gene prediction. Since lncRNAs can cis target the neighboring genes 88,89 , we searched coding genes 10-Kb upstream and downstream of all identified lncRNAs and analyzed their functional roles. Firstly, we computed Pearson's correlation coefficients between each pair of lncRNA-protein coding genes. The protein coding genes having significant correlations with lncRNAs at p.BH < 0.05 were considered potential cis target genes for those lncRNAs. The trans-acting lncRNAs influence the expression of other genes. For the trans target gene prediction, the blast ratio (e < 1E-5) between lncRNAs and protein coding genes was calculated. RNAplex software was then used to select trans-acting target genes 90 . GO and KEGG pathway enrichment analysis. To improve our understanding of the functions of the lncRNAs and mRNAs, mRNAs and potential cis and trans target genes of lncRNAs were mapped to GO terms in the Gene Ontology (GO) database, and the gene numbers for each GO term were calculated using the GOseq R package (v1.18.0) 91 . In addition, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of target genes and mRNAs was performed using the KOBAS software (v2.0) 92 using a hypergeometric test. GO terms and KEGG pathways with P < 0.05 were considered significantly enriched. The calculating formula is: where ʻNʼ is the number of all lncRNAs with GO or KEGG annotation; ʻnʼ is the number of lncRNAs in ʻNʼ; ʻMʼ is the number of genes that are annotated to the certain GO terms or specific pathways; and ʻmʼ is the number of DEGs that are annotated to the certain GO terms or specific pathways in ʻMʼ.
qRT-PCR. Several differentially expressed mRNAs and lncRNAs involved in HF development were selected and confirmed by qRT-PCR with Ovis-actin used as an internal reference. The primers for qRT-PCR are listed in Supporting Information 1. Total cDNA was synthesized using ReverTra Ace qPCR kit (TOYOBO). qPCR was performed on the 7900 HT Sequence Detection System (ABI, USA) using ABI Power SYBR Green PCR Master Mix (ABI, USA). Quantitative PCR was performed under the following conditions: denaturation at 98 °C for 5 min, followed by 40 cycles of amplification and quantification (95 °C for 15 s, and 60 °C for 1 min). Melting curves were analyzed after amplification. Except for the cDNA conversions, the other reactions were run in triplicate for each sample. The delta-delta-Ct method was used to analyze the expression level 93 .