Identification and analysis of differentially expressed long non-coding RNAs between multiparous and uniparous goat (Capra hircus) ovaries

Long non-coding RNAs (lncRNAs) play important roles in almost all biological processes. However, there is little information on the effects of lncRNAs on ovulation and lambing rates. In the present study, we used high-throughput RNA sequencing to identify differentially expressed lncRNAs between the ovaries of multiparous (Mul) and uniparous (Uni) Anhui White goats. Among the 107,255,422 clean reads, 183,754 lncRNAs were significantly differentially expressed between the Uni and Mul. Among them, 455 lncRNAs were co-expressed between the two samples, whereas, 157,523 lncRNAs were uniquely expressed in the Uni, and 25,776 uniquely lncRNAs were expressed in the Mul. Through Cis role analysis, 24 lncRNAs were predicted to overlap with cis-regulatory elements, which involved in Progesterone-mediated oocyte maturation, Steroid biosynthesis, Oocyte meiosis, and gonadotropin-releasing hormone (GnRH) signaling pathway. These 4 pathways were related to ovulation, and the KEGG pathway analysis on target genes of the differentially expressed lncRNAs confirmed this results. In addition, 10 lncRNAs harbored precursors of 40 miRNAs, such as TCONS_00320849 related to a mature miRNA sequence, miR-365a, which was reported to be related to proliferation, were annotated in the precursor analysis of miRNAs. The present expand the understanding of lncRNA biology and contribute to the annotation of the goat genome. The study will provide a resource for lncRNA studies of ovulation and lambing.

Introduction Long non-coding RNAs (lncRNAs) are non-coding RNA transcripts of more than 200 nucleotides in length. Initially, lncRNA was regarded as transcriptional "noise" that lacked a complete open reading frame and had no ability to encode proteins [1]. Human homologue 19 (H19), found in mice, was the first lncRNA which was discovered in mammals in 1998 [2], then many more lncRNAs were discovered in the following 20 years. Recent studies have shown that lncRNA participates in many important physiological processes, such as X chromosome inactivation [3], dosage compensation [4], and gene imprinting [5]. Additionally, lncRNAs have a significant impact on the diagnosis, treatment, and development of disease. For example, lncRNA plays an important role in epigenetics [6,7], post-transcriptional regulation [8], maintenance of stem cell pluripotency, and regulation of gene expression during disease processes [9].
Although lncRNAs are widely found in animals, their mechanisms of action are unclear [10,11]. Using high-throughput poly (A)-independent and strand-specific RNA-seq of rats, lncRNA expression was found to be more tissue specific than that of mRNA [12]. Knockdown of several lncRNAs in mature oocytes increases developmental rates in cattle and leads to larger blastocysts [13]. The sixth CCCTC binding-factor (CTCF) binding site (CCCTC) was identified in H19 DMR, however, had a significant methylation difference between the highand low-fertility bulls [14]. Boulanger promised that Forkhead box L2 (Foxl2) loss of function dissociated from loss of lncRNA expression is sufficient to cause an XX female-to-male sex reversal in the goat model and, as in the mouse model, an agenesis of eyelids [15]. Ovarian hormones regulate Foxl2 expression and thereby expand the number of genes controlled by the hypothalamic-pituitary-gonad axis that ultimately dictates reproductive fitness [16]. Although lncRNA studies have been begun in several species, the research conducted into the effects of lncRNA on the reproductive rate was relatively limited. High-throughput sequencing can be used to identify the lncRNAs, and bioinformatics analysis can predict the function of lncRNAs, which would provide a basis for animal reproduction study and possibly improve the reproductive rate of goats.
Compared with other goats, Anhui White goat (AWG) is known for its higher fertility, precocious puberty, and higher leather quality. The kidding rate of AWG is 2.27-2.39 and the ewes can be in estrus all year round, which makes it to be an ideal model for the study of goat breeding traits [17]. To explore the function of lncRNAs in ovulation and lambing, we used Solexa sequencing to identify differentially expressed lncRNAs between the ovaries of multiparous and uniparous AWGs in the study, and predict the target genes of the differentially expressed lncRNAs. In addition, classified annotation, GO analysis and KEGG pathway were used to analyze the function of target genes. Through Cis role analysis, 24 lncRNAs were predicted to overlap with cis-regulatory elements, which involved in Progesterone-mediated oocyte maturation, Steroid biosynthesis, Oocyte meiosis, and gonadotropin-releasing hormone (GnRH) signaling pathway. These 4 pathways were related to ovulation, and the KEGG pathway analysis on target genes of the differentially expressed lncRNAs confirmed this results. These results will help us to further understand the role of lncRNAs in ovulation traits.

Experimental animal preparation and library construction
The experimental goats in this study, Anhui White goats (AWG, an indigenous Chinese breed), were obtained from the College of Animal Science and Technology, Anhui Agricultural University, Hefei, China. Based on the long-term observation, laparoscopic and ultrasound detection, 6 target goats, 3 multiparous goats (Mul) and 3 uniparous goats (Uni), which were undergoing oestrus, were selected for their all ovaries. According to the at least 3-year records of lambing, the litter size of Mul and Uni were more than one and only one, respectively. Meanwhile, the physical condition and age, which is 3 to 4 years old, of the experimental samples were basically consistent. To further eliminate the influence of other factors, the unified management system of the field was adopted in the study farm feeding and stabling. The selected goats were killed and dissected, and both whole ovaries from each goat were collected immediately, snap-frozen in liquid nitrogen, and stored at -80˚C until use. One whole ovary from each goat was used for RNA-Seq, meanwhile, the other one from the same goat was used for real-time PCR. All experimental procedures involving AWGs performed in the present study had been given prior approval by the ethics committee of Anhui Agricultural University under permit no. AHAU20101025.
The total RNA of the ovaries (12 samples) from the 6 goats was extracted using RNAiso Plus (Takara) following the manufacturer's protocol. The quantity and quality of the total RNA were measured using the Agilent 2100 Bioanalyzer system. To minimize differences among the goats, an equal quantity (10 μg) of total RNA isolated from one ovary from the 3 individual multiparous and uniparous goats, respectively, were pooled for library preparation and Illumina sequence (The Beijing Genomics Institute). In addition, the total RNA from the other one of the ovaries from the 3 individual multiparous and uniparous goat were reverse transcribed, respectively, into cDNA using the AceQ qPCR SYBR Green Master Mix lncRNA RT-PCR Kit (Vazyme Biotech Co., Ltd.) for qRT-PCR verification.

Basic data statistics
Before further analysis, raw data filtering was needed to decrease data noise. "Raw reads" were defined as reads containing the adapter sequence, a high content of unknown bases (reads with more than 10% unknown bases), and low-quality reads (>50% low-quality bases in a read, with a low-quality base defined as a base whose sequencing quality was no more than 10). Then clean reads were mapped to the rRNA reference sequence using SOAP aligner/ SOAP2 short-read alignment software [18] to remove the remaining rRNA reads, and leave the reads that were used for transcriptome assembly and quantification. The statistical analysis of the alignment results is presented for each sample.

Assembly and novel lncRNA prediction
Reads were mapped to the genome and assembled using Cufflinks [19]. Faux-reads were then generated from reference transcripts to capture features in the reference sequence that, due to low coverage, might be missing from the sequencing data. We used a reference annotationbased assembly method to reconstruct the transcripts, while background noise was reduced by using the FPKM (fragments per kilobase of exon per million fragments mapped) value and coverage threshold. These reads were then merged with the sequenced reads for assembly. The set of transcripts generated in the last step was then compared with the reference transcripts to remove the transcripts that were approximately equivalent to the whole or a portion of a reference transcript. Through comparison with the reference, we were able to detect the novel transcripts and calculate the coding potential of these transcripts to identify the novel lncRNA [20]. After the assembly, we compared the assembled transcripts to the reference annotation by Cuffcompare, and categorized the transcripts to 12 classes (=, c, j, e, i, o, p, r, u, x, s and.). Among them, only five candidate categories of transcripts will be extracted, including 'i', 'j', 'o', 'u', and 'x', which may contain novel transcripts. The "u" category meant unknown intergenic transcript. The "x" category meant lncRNAs with exonic overlap with known transcripts, but on the opposite strand. The "i" category contained lncRNAs with transcripts falling entirely within a reference intron. The "j" category meant potentially novel isoforms (fragments) with at least one splice junction shared with the reference transcripts.

Classification and annotation of lncRNAs
The lncRNAs were also classified by their position compared with the reference gene, with different strategies used for function prediction. For antisense lncRNA, we performed free energy calculation to detect possible hybridization sites for lncRNA and mRNA, which might be RNA-RNA interactions. To further reveal potential antisense lncRNA-mRNA interactions, we searched for all antisense lncRNA-mRNA complementary base pair duplexes using RNAplex [21], a tool specially created to rapidly search for short interactions between two long RNAs. If such interactions are found upstream or downstream of a gene, we will search for lncRNAs located in unknown regions. Because lncRNA can be processed to yield small RNAs, small RNA precursors can be predicted. Recent genome-wide studies suggested that the function of a significant fraction of long non-coding transcripts may be to serve as precursors for micro-RNAs (miRNAs) that are usually generated via the sequential cleavage of long transcripts by the enzymes Drosha and Dicer [22,23]. Thus, we aligned lncRNAs to miRBase [24] to detect potential pre-miRNAs by selecting those with a hit coverage greater than 90%.
The SVM-based software miRPara [25] was also used to predict probable miRNAs. In addition, Rfam divides non-coding RNAs into families based on evolution from a common ancestor. To better annotate novel lncRNAs from an evolutionary standpoint, we classified all predicted novel lncRNAs into different non-coding RNA families using INFERNAL, which categorizes non-coding RNAs and their conserved primary sequence and RNA secondary structure through the use of multiple sequence alignments, consensus secondary structure annotation, and covariance models [26]. These results will classify lncRNA families by their RNA structure and sequence similarities, helping us to reveal the potential functions of the lncRNAs.

Screening of differentially expressed lncRNAs and real-time PCR validation
Cuffdiff was used to calculate expression levels of lncRNAs in more than one condition and test them for significant differences on the basis of FPKM values. To screen differentially expressed genes from the Cuffdiff results, we followed the priority rule: STATUS = OK and P 0.05. The differentially expressed lncRNAs were screened using NOIseq [27] following the priority rule: filtering condition, log2 (ratio) ! 1; probability ! 0.8.
Real-time PCR (qPCR) was performed to validate the Solexa sequencing data. Eight differentially expressed lncRNAs between the two libraries were randomly selected for qPCR verification. The primers were designed based on the lncRNA sequences (S1 Table). After incubation at 37˚C for 1 h and deactivation at 95˚C for 5 min, we obtained the template for qPCR. The reaction solution of qPCR contained 2.0 μl cDNA, 10 μl AceQ qPCR SYBR Green Master Mix and ROX Reference Dye 1, 0.4 μl of each primer, and 6.8 μl ddH 2 O. GAPDH was used as the internal reference gene control. qPCR was performed using standard protocols on the StepOnePlus™ Real-Time PCR System (Thermo Fisher Scientific). The threshold cycle (CT) was obtained from each reaction, and the relative expression level of each lncRNA was evaluated using the 2 -ΔΔCt method. Three biological replicates for each selected lncRNA were used. The expression levels of lncRNA in the ovaries of Mul and Uni were determined individually and the t-test was used to examine the significant of expression differences between the two samples using SAS v8.0 software.

Target genes prediction of differentially expressed lncRNAs and function analyses
To explore the functions of lncRNAs, we predicted the target genes of the differentially expressed lncRNAs through Cis role. Cis role is lncRNA action on neighboring target genes. The basic principles of Cis role of target gene is prediction the lncRNA function and the protein coding gene which near the lncRNA's coordinates. We searched coding genes 100 kb upstream and downstream of lncRNA and then analyzed their function by functional enrichment analysis.
Gene Ontology (GO) and KEGG pathway analyses can help us to understand the biological functions of genes [28]. GO enrichment analysis of lncRNA target genes were implemented by the GOseq R package (http://www.bioconductor.org/packages/release/bioc/html/goseq.html), in which gene length bias was corrected. GO terms with corrected P-value less than 0.05 were considered significantly enriched by differential expressed genes. KEGG is database resource for understand high-level functions and utilities of the biological system, such as the cell, the organism and the ecosystem, from molecular-level information. We used KOBAS software (http://kobas.cbi.pku.edu.cn/download.php) to test the statistical enrichment of lncRNA target genes in KEGG pathways. Pathway with corrected FDR-value less than 0.05 were considered significantly enriched by differential expressed genes.

Basic data statistics
In the base composition of the reads, on the X-axis, 1-90 bp represents "read 1" and 91-180 bp represents "read 2". The A curve should overlap with the T curve, while the G curve should overlap with the C curve. In this case, the figure showed a balanced composition (S1 Fig). By comparing with the rRNA reference sequence through SOAP aligner/ SOAP2 short-read alignment software, 52,178,296 and 55,077,126 clean reads in the Uni and Mul libraries, respectively, were mapped to the genome. Only 0.02% and 0.03% of the total sequence reads were mapped to single-end in the two libraries, respectively. Additionally, only 0.01% of the total sequence reads were mapped to paired-end mapping reads in the two libraries (S2 Table). Otherwise, about 80% of the total mapped reads were mapped to the reference genome, indicating that the quality of the two libraries is good (S3 Table).

Assembly and novel lncRNA prediction
The clean reads were mapped to the genome using TopHat (Fig 1). There were 663,188 genes were identified in the gene sequencing coverage statistics of the Uni. Among them, transcript coverage on the genome of 63,140 genes was more than50%, which accounted for 9.24% of the identified genes. Meanwhile, 675,422 genes were identified in the Mul library, and 383,043 genes (accounting for 56.71%), whose transcript coverage on the genome was more than 50%. In total, 511,721 lncRNAs and 13,875 mRNAs were identified in the two libraries. Among them, 271,584 lncRNA and 13,224 mRNA were co-expressed in the two libraries, and 240,137 lncRNAs and 651 mRNAs were respectively expressed in the certain library (Table 1).
Cuffcompare categories were used to distinguish transcripts according to 12 class codes (S4 Table). Only five candidate categories of transcripts will be extracted, including "i", "j", "o", "u" and "x" which may contain novel transcripts. The "u" category contained 525,553 transcripts, which meant unknown intergenic transcript lncRNAs (lincRNAs). The "x" category contained 4,053 transcripts, which meant lncRNAs with exonic overlap with known transcripts, but on the opposite strand. The "i" category contained lncRNAs with transcripts falling entirely within a reference intron. The "j" category meant potentially novel isoforms (fragments) with at least one splice junction shared with the reference transcripts, which contained 38,842 transcripts, and 102 transcripts in "o" category. Due to the complexity of the repeat region, the probability of transcript misassembly is also increased. Thus, we filtered transcripts classified with "r", which could possibly be artificial fragments. As determined by their genomic locations with respect to nearby genes, the 568,550 lncRNAs included members of all four classes of lncRNAs, with most in the u class, which overlap unknown intergenic transcript lncRNAs. The results of the lncRNA triage were used for further analysis.

Classification and annotation of lncRNA
By examining the up/downstream lncRNAs of a gene, we found 27,228 lncRNAs had up/downstream gene. Of these, 13,069 lncRNAs were located in the antisense strand and 14,159 were located in the sense strand. Therefore, 49% of lncRNAs were downstream of the bracketing genes, and the others were upstream. These lncRNAs could overlap with cisregulatory elements probably involved in transcriptional regulation. TCONS_00076241, TCONS_00087806, and TCONS_00076240 were upstream of XM_005686239.1 (BUB1) in the oocyte meiosis pathway. TCONS_00248477 was downstream of XM_005676195.1 (MAP3K19) in the GnRH signaling pathway. TCONS_00136407 and TCONS_00146968 were downstream of XM_005689083.1 (MOS) in the progesterone-mediated oocyte maturation pathway. TCONS_00136407 and TCONS_00146968 were downstream of XM_005689083.1 (MOS) in the oocyte meiosis pathway. MOS, BUB1, and MAP3K19 may regulate oocyte meiosis ( Table 2).
To better annotate the novel lncRNA from an evolutionary standpoint, we classified all of the predicted novel lncRNAs into different non-coding RNA families using INFERNAL. Fifty families were obtained via the lncRNA family prediction. The C/G content was between 25% and 70%. Xist_exon1, TUG1_3, TUG1_4, SOX2OT_exon3, and SMAD5 may regulate reproduction and embryonic development [29][30][31] (Table 3).

Differentially expressed lncRNAs and real-time PCR Validation
Compared to the Uni library, 26,187 up-regulated and 157,567 down-regulated lncRNAs were identified in the Mul library (Fig 3). A total of 183,754 lncRNAs were significantly differently expressed between the Mul and Uni with |log2Ratio|2 ! 1 and p-value 0.05. Among them, 455 lncRNAs, which included 188 known genes and 267 predicted transcripts, were coexpressed in the two samples. In the Mul, the highest expressed lncRNA from the co-expressed lnRNAs was XLOC-013414, which FPKM value was 5052.08, whereas, its FPKM in the Uni was 572.03. In the Uni, the highest expressed lncRNA from the co-expressed lnRNAs was 102188530, which FPKM value was 741.67, whereas, its FPKM in the Mul was 109.21. In addition, 157,523 lncRNAs were uniquely expressed in the Uni, and 25,776 uniquely lncRNAs were expressed in the Mul. The expression of uniquely expressed lncRNA was relatively low. Such as the expression of XLOC-505751, which was highest uniquely expressed in the Uni, was 60.29. And the expression of XLOC-263636, which was highest uniquely expressed in the Mul, was 70.38 (S7 Table).
To ensure the accuracy and reliability of the RNA-seq results, eight differentially expressed lncRNAs were randomly selected for qPCR to confirm the expression changes. The results showed good consistency between the RNA-seq results and the qPCR results (Fig 4).

Function analyses of target genes of differentially expressed lncRNAs
The target genes of the differentially expressed lncRNAs were predicted through Cis role. After searching the coding genes 100 kb upstream and downstream of lncRNA, we analyzed the function of the target genes through GO and KEGG enrichment analysis. In total, there Differentially expressed long non-coding RNAs between multiparous and uniparous goat  Start, start site of the alignment; End, end site of the alignment; E-value, expect value for the alignment; Score, score of the alignment given by Infernal; GC %, GC content of the conserved region; Rank, Rank of this lncRNA in the corresponding family.
https://doi.org/10.1371/journal.pone.0183163.t003  Differentially expressed long non-coding RNAs between multiparous and uniparous goat were 32,627 genes were predicted as the target genes of differentially lncRNAs between the Uni and Mul (S8 and S9 Tables). Based on the gene background, there were 439, 432, and 407 target genes clustered into 141, 198, and 868 GO terms in component, function and process ontology, respectively. The GO terms "Cell part", "Binding" and "Cellular process" were the biggest terms, which included the most target genes, in the above ontologies, respectively. However, compared with the reference gene background, no GO term was significantly enriched (p-value 0.05). The "Phosphatidylinositol binding" and "Cyclic-nucleotide phosphodiesterase activity" in function ontology, and "Cholesterol metabolic process", "Cell aging", "Aging", "Steroid metabolic process" in process ontology were the only GO terms with corrected p-value < 1 (S10 Table).

Discussion
In the study, we used RNA-seq to study the lncRNAs expressed in the ovaries of multiparous (Mul) and uniparous (Uni) goat to explore the function of lncRNAs in ovulation and lambing. In total, 107,255,422 clean reads were obtained from the two libraries, and 52,178,296 and 55,077,126 clean reads in the Uni and Mul libraries, respectively, were mapped to the genome. However, a small part of the lncRNA sequence still could not be aligned to the target genome and located to a specific chromosome. Because the sequence was broken up into smaller pieces called segments, with the software inferring that the read spans a splice junction and estimating the location of the splice sites of the junction. On one hand, we believed that this shortcoming may be because the reference goat genome is still not perfect, with more pieces to be spliced. On the other hand, there may be post-transcriptional RNA editing, resulting in inconsistencies between the lncRNA sequence and the corresponding DNA sequence. Some lncRNAs highly expressed in the ovary may lack a classical poly (A) tail [35]. In the renal medulla, 28% of lncRNAs were found to lack a classical poly (A) tail [12]. Strand-specific sequences rather than 3'-directed cDNA sequences were used in this study. Strand-specific sequences help to determine the polarity of transcripts and improve the accuracy of the quantification of antisense transcripts. When lncRNA and mRNA were isolated from total RNA, the ribo-depleted method was used instead of the poly (A) enrichment method. Large-scale bioinformatic studies suggest that a significant fraction (>24%) of long non-coding transcripts present in cells may lack a classical poly (A) tail.
In this study, we searched for lncRNAs located in unknown regions if they could be found upstream or downstream of a gene to predict the function of the lncRNA. TCONS_00076241, TCONS_00087806, and TCONS_00076240 were upstream of XM_005686239.1 (BUB1). TCONS_00248477 was downstream of XM_005676195.1 (MAP3K19). TCONS_00136407 and TCONS_00146968 were downstream of XM_005689083.1 (MOS). Studies have shown that XM_005689083.1 and XM_005686239.1 play a role in meiosis I in the oocyte [36]. XM_005676195.1 may regulate GnRH in the GnRH signaling pathway. Bta-mir-365-1, cfa-mir-493 [37], and bta-mir-493 [38] are precursors of TCONS_00284825 and oar-mir-493 [39] is a precursor of TCONS_00320849 in pre-miRNA prediction. TCONS_00320849 has a mature miRNA sequence, miR-365a, with the same mature sequence in the miRBase database. miR-365 can regulate proliferation by mediating the expression of KRAS [40]. miR-365 regulates cell cycle progression and apoptosis in various types of tumors [41]. Thus TCONS_00136407, TCONS_00146968, and TCONS_00320849 may also play roles in meiosis in the oocyte.
To further detect the lncRNAs, which may participate in the ovulation and lambing of goat, we analyzed the differentially expressed lncRNAs between the Uni and Mul. A total of 183,754 lncRNAs were significantly differently expressed between the two libraries. Among them, 455 lncRNAs, which included 188 known genes and 267 predicted transcripts, were co-expressed between the two samples. In addition, 157,523 lncRNAs were uniquely expressed in the Uni, and 25,776 lncRNAs were uniquely expressed in the Mul. Accordingly, we analyzed the lncRNAs which have high expression levels and were significantly difference expressions between the Uni and Mul libearies. XLOC_013414 was differentially expressed lncRNA with the highest expression level in the Uni library, whereas 102188530 (COL1A1) had the highest expression in the Mul library. The function of XLOC_013414 is unknown but COL1A1 polymorphisms may individually play minor roles in osteoporosis and fracture [42]. Meanwhile, XLOC_500835, XLOC_075014, XLOC_523784 were significantly differentially expressed novel lncRNA between the Uni and the Mul. And their target genes were KIAA0368, CNGA3 and MID2, respectively. KIAA0368 as selectively up-regulated by contact allergene, may be potential markers of skin irritation and allergy [43]. Achromatopsia can be caused by mutations in the CNGA3 gene [44]. MID2 leads to misregulation of microtubule organization and the downstream disease pathology associated with X-linked intellectual disabilities [45]. In consideration of the specific and higher expression levels of XLOC_013414 and COL1A1, some further research of their roles in the ovary is warranted.
According to the KEGG pathway analysis, a total of 875 target genes, which showed 22,885 background genes, of the differentially expressed lncRNAs between the Uni and Mul were annotated to 236 pathways. Among them, 12 pathways were significantly enriched (p-value 0.05). In addition, five pathways, "Progesterone-mediated oocyte maturation", "Steroid biosynthesis", "Oocyte meiosis", "GnRH signaling pathway", and "Steroid hormone biosynthesis" were related to the animal reproduction. An increased FSH secretion level will lead to multiple follicular development and ovulation in the follicular phase [46][47][48]. In addition, Neat1 knockout mice fail to become pregnant despite normal ovulation [49]. Ovarian hormones can regulate the expression of Foxl2 to expand the number of genes controlled by the hypothalamicpituitary-gonadal axis, ultimately dictating reproductive fitness [16].
In summary, the differentially expressed profile of lncRNAs was successfully established in the Uni and Mul phase libraries. Through Cis role analysis, 24 lncRNAs were predicted to overlap with cis-regulatory elements, which involved in Progesterone-mediated oocyte maturation, Steroid biosynthesis, Oocyte meiosis, and gonadotropin-releasing hormone (GnRH) signaling pathway. Meanwhile, the KEGG pathway analysis on target genes of the differentially expressed lncRNAs confirmed this results. In addition, 10 lncRNAs harbored precursors of 40 miRNAs, which was reported to be related to proliferation, were annotated in the precursor analysis of miRNAs. The present expand the understanding of lncRNA biology and will provide a resource for lncRNA studies of ovulation and lambing.   Table. The protein-coding genes within 100 k upstream of an lncRNA. (XLSX) S10 Table. 10 GO analysis for the target genes of the differentially expressed lncRNAs. (XLSX) S11 Table. KEGG pathways for the target genes of the differentially expressed lncRNAs. (XLS)