Transcriptome analysis in the shell gland of laying hens affecting eggshell qualities

Background Eggshell plays an important role in protecting against physical damage and microorganic invasion. It is subject to quality loss with increasing hen age, and fragile eggshells result in huge economic losses to the poultry industry. Therefore, improving eggshell quality is particularly important. However, little is known about the potential molecular mechanisms regulating eggshell quality in chickens. In this study, we aimed to compare differential expression of long non-coding RNAs (lncRNAs) and mRNAs between old and young laying hens to identify related candidate genes for chicken shell gland development by the method of high-throughput RNA sequencing (RNA-seq). In total, we detected 176 and 383 differentially expressed (DE) lncRNAs and mRNAs, respectively. Moreover, functional annotation analysis based on the Gene Ontology (GO) and Kyoto encyclopedia of genes and genomes (KEGG) databases revealed that DE-lncRNAs and DE-mRNAs were signicantly enriched in “phosphate-containing compound metabolic process”, “mitochondrial proton-transporting ATP synthase complex”, “inorganic anion transport”, and other terms related to eggshell calcication and cuticularization. Through integrated analysis, we found that some important genes such as FGF14, COL25A1, GPX8, and GRXCR1 and their corresponding lncRNAs were expressed differentially between two groups, and the results of quantitative real-time polymerase chain reaction (qPCR) among these genes were also in excellent agreement with the sequencing data. In addition, our research indicates that FGF14, COL25A1, GPX8, and the members of the SLC family may be key genes that affect eggshell quality in hens. This study provides a catalog of lncRNAs and mRNAs of the laying hen eggshell gland and will contribute to a fuller understanding of the molecular mechanisms of the function of the shell gland in poultry. Our ndings will provide a valuable reference for the development of breeding programs aimed at breeding excellent poultry with high eggshell quality or regulating dietary nutrient levels to improve eggshell quality. holoenzyme bridging (GO:0030674)”. Several genes were in most notably FGF14 and Another important group of DEGs were involved in membrane ber formation and/or encoded extracellular matrix proteins; “extracellular space (GO:0005615)”, “membrane (GO:0016020)”, “brinogen complex (GO:0005577)” were implicated in this.

in laying hens. Therefore, in this study, we performed a transcriptomic analysis of the shell gland among old and young laying hens to identify related mRNAs, lncRNAs and pathways. We captured both lncRNAs and mRNAs from fragmented or intact RNA samples to compare whole transcriptomes of old and young chicken shell glands at unprecedented depth. Then, the differentially expressed (DE)-lncRNAs were used in bioinformatics analyses to predict cisand transtarget genes and to construct lncRNA-mRNA co-expression interaction networks. Next, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional analyses were performed to investigate the related roles of differentially expressed genes (DEGs). We hypothesized that this approach would lead to the identi cation of only critical genes and pathways relevant to the two age groups examined for shell gland development. In this way, the present study provides predictions about the interactions of lncRNAs and mRNAs, and the information generated from these predictions can be utilized in future studies of lncRNA function during chicken shell gland development at different ages.

Reads Mapping
In total, we obtained 82,871,160 − 86,696,604 and 86,622,130 − 86,688,990 raw reads from the libraries from shell gland tissues of old chickens (n = 4) and young chickens (n = 4), respectively. Correspondingly, we ultimately obtained 80,510,552 − 138,847,948 and 84,918,798 − 85,469,778 clean reads by ltering and removing sequence reads with adapters and low quality, respectively. In addition, the Q30 of each sample was not < 90.85%. Subsequently, we found that > 78.77% of the clean reads completely mapped to the chicken reference genome. The unique mapped reads ranged from 66. .46% of the total mapped reads ( Table 2).

Identi cation and Characterization of LncRNAs
We performed a comparative analysis of the structure of lncRNAs and mRNAs to study the basic features of lncRNAs in the chicken shell gland. This was not just to determine the difference between lncRNAs and mRNAs but also to verify if the predicted lncRNAs were consistent with general characteristics. In this study, the intersection of the Coding Potential Calculator (CPC), Coding-Non-Coding Index (CNCI), and Protein Families Database (PFAM) results yielded 5,334 lncRNA transcripts, including the identi ed conservative lncRNAs (Fig. 1A). Interestingly, previous reports indicate that protein coding transcripts are longer in length and more conserved than lncRNAs [13]. In agreement with this, we found that the predicated lncRNAs are shorter in length than protein coding transcripts (Fig. 1B) and tend to contain fewer exons (Fig. 1C). We also found that the average open reading frame (ORF) length of the predicted lncRNAs was 126 amino acids (aa), which was less than mRNA (687 aa, Fig. 1D).

Differential Expression of Predicted LncRNAs and mRNAs in the Eggshell Gland
The expression level of each lncRNA and mRNA was estimated by FPKM using Cuffdiff. To explore similarities and to compare the relationships between the different libraries, we measured the expression patterns of DE-lncRNAs and protein-coding genes by systematic cluster analysis (Fig. 2). As a result, we identi ed 176 lncRNA transcripts that were expressed differentially in the eggshell glands between the old group and young group (Supplementary Table S1), and the sequences could be found in the Supplementary le 1. Compared to the young group, 91 lncRNAs were up-regulated, and 85 lncRNAs were downregulated, in the old group. Among these, the 20 most signi cantly up-regulated or down-regulated lncRNAs are presented in Table 3 ( Fig. 2A, 2B and Table 3). Differential expression of mRNAs in shell gland tissues of the old group was also compared to that in the young group. A total of 383 mRNAs were found to be expressed differentially, with 204 up-regulated and 179 down-regulated (Fig. 2C, 2D, and Supplementary Table S2).

Construction of the LncRNA-mRNA Co-expression Network
To investigate the questions of whether the functions of DE-lncRNAs are in agreement with their target genes in regulating the chicken eggshell gland, and how do lncRNAs and their target genes interact (cis or trans), we constructed a co-expression network between DE-lncRNAs and their signi cantly correlated DE cisand trans-target genes using Cytoscape (Fig. 3). For the old chicken vs. young chicken comparison, the lncRNA-mRNA co-expression interaction network comprised 37 network nodes and 48 lncRNA-gene connections among 13 DE-lncRNAs and 24 DE-mRNAs. In addition, both TCONS_00181492 and TCONS_03123639 regulated their target genes FGF14 and GRXCR1 in cis, as shown in the co-expression network. As seen in Fig. 3, one mRNA may correlate with one to four lncRNAs, and one lncRNA may correlate with one to six mRNAs.

Enrichment Analysis of the Nearest Neighbor Genes of the LncRNAs
To investigate the functions of the lncRNAs, we predicted their potential cis targets. We searched for protein-coding genes 10 kb and 100 kb upstream and downstream of all of the identi ed lncRNAs. We found 176 lncRNAs that were transcribed close to (< 10 kb) 206 neighboring protein-coding genes, and 176 lncRNAs that were transcribed close to (< 100 kb) 154 neighboring protein-coding genes (Supplementary Tables S3 and S4). To explore the functions between lncRNAs and their cis-regulated target genes, we performed GO analysis. We found 90 GO terms (< 10 kb) that were signi cantly enriched (p < 0.05) ( Supplementary Table S5), and most of these terms were associated with biological processes and molecular functions ( Supplementary Fig. S1). In addition, we found 140 GO terms (< 100 kb) that were signi cantly enriched (p < 0.05) (Supplementary Table S6), and most of these terms were associated with biological processes, molecular functions, and cellular components ( Supplementary Fig. S2). For example, the main enriched terms included "protein phosphorylation (GO:0006468)", "phosphate-containing compound metabolic process (GO:0006796)", "phosphorus metabolic process (GO:0006793)", "protein modi cation process (GO:0006464)", "ATP binding (GO:0005524)", "ATP-dependent helicase activity (GO:0008026)", and "mitochondrial protontransporting ATP synthase complex (GO:0005753)" (Tables 4 and 5). Most of them were closely related to the formation of eggshells, which suggests that one of the principal roles of lncRNAs may be to regulate the synthesis and metabolism of organics and minerals. Pathway analysis indicated that co-location genes were signi cantly enriched in four (< 10 kb) and six (< 100 kb) KEGG pathways (p < 0.05), respectively (Tables 6 and 7). These data suggest that the function of the shell gland may be regulated by the action of lncRNAs on these neighboring protein-coding genes.  Note: "Input number" represent DE-lncRNAs corresponds to the gene number associated with the PATHWAY. This notation is also applicable to Table 7, 9 and11. We also predicted the potential targets of lncRNAs in trans-regulatory relationships using co-expression analysis. The COR method was used to analyze the correlation between the lncRNAs and mRNAs in samples, and the main functions of the lncRNAs were predicted using mRNA, with a correlation absolute value > 0.95. We found 176 lncRNAs that were transcribed close to 791 protein-coding genes (Supplementary Table S7). Functional analysis indicated that the co-expressed genes were signi cantly enriched in 174 GO terms (95 under biological process, 43 under cellular component, and 36 under molecular function) that encompassed a variety of biological processes (p < 0.05) (Supplementary Table S8 and Fig. S3). Importantly, some of the terms were related to organic metabolism and genetic development, including "cellular protein metabolic process (GO:0044267)", "macromolecule biosynthetic process (GO:0009059)", "protein metabolic process (GO:0019538)", "Ras GTPase binding (GO:0017016)", and "GTPase binding (GO:0051020)" ( Table 8). Most of them were associated with organic synthesis and metabolism. The co-expressed genes were signi cantly enriched in nine KEGG pathways (p < 0.05) ( Table 9), where the pathways "Salmonella infection" and "AGE-RAGE signaling pathway in diabetic complications" affected the function of the shell gland of aging laying hens.
As the disease resistance of aging hens is weakened, the metabolism and synthesis ability of the body are reduced, which leads to a decline in eggshell quality.   Table S9), and most of these terms were associated with biological process, molecular function, and cellular component (Table 10 and Supplementary Fig. S4). The majority of DEGs were categorized as ion transport in the eggshell gland during formation of the eggshell. The GO terms included "inorganic anion transport (GO:0015698)", "inorganic anion transmembrane transporter activity (GO:0015103)", "anion transmembrane transporter activity (GO:0008509)", "electron carrier activity (GO:0009055)", and "calcium ion binding (GO:0005509)". Thirty-one genes were categorized under these terms, including Glutaredoxin cysteine-rich 1 (GRXCR1) and the members (SLC1A3, SLC6A4, SLC20A1, SLC22A13, SLC26A3, SLC30A8, SLC39A2, SLC43A3, and SLC45A2) of the sodium-dependent phosphate transporter family. GO term analysis also revealed some DEGs with possible roles in protein translation and binding. The terms included "protein polymerization (GO:0051258)", "regulation of G-protein coupled receptor protein signaling pathway (GO:0008277)", "transcription factor complex (GO:0005667)", "DNA-directed RNA polymerase II, holoenzyme (GO:0016591)", and "protein binding, bridging (GO:0030674)". Several genes were enriched in these terms, most notably FGF14 and COL5A2. Another important group of DEGs were involved in membrane ber formation and/or encoded extracellular matrix proteins; "extracellular space (GO:0005615)", "membrane (GO:0016020)", and " brinogen complex (GO:0005577)" were implicated in this. In addition, we also found 10 KEGG pathways that were signi cantly enriched (p < 0.05) ( Table 11), several of which were related to the function of the shell gland, including "Glycine, serine, and threonine metabolism", "ABC transporters", and "Toll-like receptor signaling pathway". The D-3-phosphoglycerate dehydrogenase (PHGDH) gene was signi cantly enriched in the serine metabolism pathway, and the osteopontin (SPP1) gene is a matrix protein that was signi cantly enriched in the "Toll-like receptor signaling pathway".

Validation of DE-lncRNAs and -mRNAs
To further validate the reliability and reproducibility of our RNA-seq data, four DE-lncRNAs (TCONS_00181492, TCONS_03234147, TCONS_03123639, and TCONS_01464392) and their corresponding target genes (FGF14, COL25A1, GRXCR1, and GPX8) related to eggshell quality were randomly selected for qPCR validation. The analysis showed that the expression tendencies of all four lncRNAs and their target genes were extremely concordant with the RNA-seq data, though the absolute fold changes differed between qPCR and RNA-seq ( Fig. 4 and Supplementary Table S10). Appreciably, TCONS_00181492, TCONS_03234147, and TCONS_03123639 were up-regulated their corresponding target genes, but TCONS_01464392 down-regulated GPX8. These results are consistent with that of the co-expression interaction network, especially for TCONS_00181492 and TCONS_03123639 regulating FGF14 and GRXCR1, respectively.

Disussion
Comparative transcriptome analyses of organ or tissues at different developmental stages can provide valuable insights into the question of how regulatory gene interaction networks control speci c biological processes and how diseases can arise [14]. Recently, increasing evidence has con rmed that lncRNAs are important regulatory factors of gene expression, regulating target genes by cis-acting (neighboring genes) or trans-acting (distant genes) mechanisms [15]. Furthermore, RNA-seq has been performed to provide an extensive lncRNA and gene expression pro le in different tissues of livestock and poultry (e.g., cell differentiation and development [16], cancer [17], and skeletal muscle development [18]). Previous studies of the hen uterus transcriptome and gene expression pro ling during formation of the eggshell demonstrate a large number of DEGs that participate in ion transport for eggshell mineralization and the secretion of matrix proteins [19][20][21][22][23]. Most of the previous studies report the roles of mRNAs in the avian eggshell gland, but systematic identi cation of the functions of lncRNAs remained unclear in the development of the chicken shell gland. Therefore, in this study, we performed transcriptome sequencing of the shell gland of laying hens in the peak and late laying periods and analyzed the DE-lncRNAs and DEGs to reveal their roles in eggshell quality. To the best of our knowledge, this study represents the rst systematic genome-wide analysis of lncRNAs and mRNAs in the chicken shell gland, providing a valuable catalog of functional lncRNAs and mRNAs associated with eggshell quality.
In the present study, we developed a highly stringent ltering pipeline to minimize the selection of false positive lncRNAs, which aimed to remove transcripts with evidence of protein-coding potential, and performed co-location mRNA prediction and co-expression mRNA prediction for the lncRNAs obtained from the chicken eggshell gland. Ultimately, we identi ed 176 DE-lncRNAs and 383 DE-mRNAs. To gain insight into how interactions between DE-lncRNAs and their corresponding target genes regulate shell gland development, we constructed co-expression interaction networks between DE-lncRNAs and their predicted cisand trans-target genes. Then, four DE-lncRNAs and their target genes related to eggshell quality were selected for qPCR validation, and the results were consistent with the RNA-seq data, which demonstrated that lncRNA TCONS_01464392 can target the GPX8 gene, and they are all down-regulated. LncRNAs TCONS_00181492, TCONS_03234147, and TCONS_03123639 target FGF14, COL25A1, and GRXCR1, respectively, and these six genes are up-regulated. Together, these results con rmed that the identi ed lncRNAs and mRNAs were of high quality.
The oviduct of hens is composed of the infundibulum, magnum, isthmus, shell gland, and vagina. Especially, the shell gland is the place where the eggshell is deposited [24]. The formation of the eggshell is a complex process involving the precipitation of calcium carbonate [25]. Mature follicles reach the shell gland and calcify layer by layer. After the mature follicles reach the shell gland, they need to go through the calci cation process, eventually form the eggshell, and the whole process takes about 15-16 h. Approximately 94% of minerals in the eggshell are calcium carbonate, with other inorganic minerals being calcium phosphate, magnesium phosphate, and magnesium carbonate [25]. Previous studies suggest that eggshell calci cation requires the interaction of numerous processes, including transcellular and paracellular transport of minerals and the secretion of different matrix proteins [26][27][28]. Particularly, ion transportation plays a crucial role in the process of eggshell formation. The ion channels contribute to the transportation of Ca 2+ from the plasma to the uterine lumen, which includes Na + , Ca 2+ , and K + channels [29]. Moreover, the characteristics of egg shell calci cation in poultry are that the body rapidly and massively transports Ca 2+ from blood to the lumen of the eggshell gland, and a calcium ATPase (calcium pump) is a key enzyme involved in Ca 2+ transport in the uterus during eggshell formation [30]. Apart from Ca 2+ , inorganic phosphate (Pi) is also essential in the formation of eggshells. Pi is involved in many biological processes, including nucleic acid synthesis, skeletal development, signaling cascades, and tooth mineralization [31][32][33]. More meaningfully, phosphorus participates in the transport mechanism of the calcium pump (calcium ATPase).
In the present study, we conducted GO and KEGG pathway enrichment analyses on DE-mRNAs and DE-lncRNAs and found that the most of identi ed DEGs were involved in eggshell calci cation and cuticularization pathways, such as "inorganic anion transport", "inorganic anion transmembrane transporter activity", "phosphate-containing compound metabolic process", "phosphorus metabolic process", "protein metabolic process", "mitochondrial protontransporting ATP synthase complex", "proton-transporting ATP synthase complex", and "calcium ion binding". Notably, SPP1 was signi cantly enriched in the "Toll-like receptor signaling pathway", and the authors of a previous study suggest that SPP1 is differentially expressed in the uterus between a low eggshell strength group and normal eggshell strength group during eggshell formation [23]. In addition, another study indicates that the PHGDH gene is highly overexpressed in the white isthmus during deposition of the eggshell membranes [19]. The PHGDH gene was also differentially expressed between two groups and enriched in the "Glycine, serine, and threonine metabolism pathway" in this study. Hence, all of these results indicate that the formulation of the eggshell is signi cantly affected by the shell gland of laying hens with different ages.
Based on the lncRNA-mRNA co-expression interaction networks, the predicted target gene of lncRNA TCONS_00181492 is FGF14. Prior to this analysis, little was known concerning the association between FGF14 and lncRNA. FGF14 is a well-known growth factor belonging to the FGF family. FGF family members possess broad mitogenic and cell survival activities and are involved in a variety of biological processes, including cell growth, embryonic development, tissue repair, morphogenesis, tumor growth, and invasion [34,35]. Previous work demonstrates that FGF14 is a functionally relevant component of the neuronal voltage-gated Na + (Nav) channel complex [36], and FGF14 can also regulate members of the presynaptic Cav2 Ca 2+ channel family [37]. Simultaneously, there is evidence that the transfer and concentration of Na + can directly affect the transportation of Ca 2+ and HCO 3 − in the chicken uterus [38].
In the present study, we found that the expression of FGF14 is up-regulated in the shell gland of chickens in the old group as compared to the young group.
The aforementioned studies indicate that the FGF14 gene plays an important role in chicken growth [39]. The predicted regulatory lncRNA, TCONS_00181492, was signi cantly more highly expressed in the shell gland in the old group than in the young group and controlled the expression of FGF14 via cis-acting mechanisms. Furthermore, TCONS_00181492 and FGF14 were positively correlated. Therefore, we had reason to speculate that TCONS_00181492 may regulate shell gland development in the chicken via the cis-acting target gene FGF14. Additionally, we found that the old hens had a higher incidence of disease than the young hens in the long-term cultivation of layers. Previous studies indicate that inherited mutations in FGF14 are linked to disease [40][41][42], and a study hints that the pathogenic effects of mutant FGF14 are likely mediated by dysregulation of both Ca 2+ and Na + channels [37]. All of these results indicate the possible role of FGF14 in aging laying hens with deteriorated eggshell quality.
COL25A1 was a predicted cis-target of TCONS_03234147 that is related to the focal adhesion pathway. Collagen XXV alpha 1 (COL25A1), the extracellular matrix gene, is a collagenous type II transmembrane protein, which was rst puri ed from senile plaques of Alzheimer's disease (AD) brains [43]. In recent years, work on collagen genes has attracted the attention of many researchers. Previous studies of the hen oviduct transcriptome during eggshell membrane formation identify a large number of differentially expressed collagen genes, such as collagen X (COL10A1), collagen I (COL1A1), collagen II (COL2A1), and collagen III (COL3A1) [19]. Moreover, COL11A1 was also differentially expressed between the normal eggshell strength group and low eggshell strength group in the study integrating transcriptome and genome re-sequencing in the chicken uterus [23]. TCONS_03234147 and its target gene COL25A1 were differentially expressed between the two groups in the present study, and their expression was higher in aging hens compared to young hens.
The GRXCR1 gene is the putative cis-target of TCONS_03123639 in the lncRNAs-genes network. The GRXCR1 gene encodes an evolutionarily conserved cysteine-rich protein with sequence similarity to the glutaredoxin family of proteins [44]. Recently, research on the function of the GRXCR1 gene has mostly been focused on diseases [45,46]. However, the biological function of the GRXCR1 gene is still rarely reported in livestock and poultry research. Herein, we found that GRXCR1 was enriched in the ion transport pathway, implying that GRXCR1 may play an important role in the formation of eggshells. Remarkably, the members (SLC1A3, SLC6A4, SLC20A1, SLC22A13, SLC26A3, SLC30A8, SLC39A2, SLC43A3, and SLC45A2) of the sodium-dependent phosphate transporter (SLC) family were also enriched in ion transport pathways (Table 10). Previous studies show that zinc ion transporters include two major families, SLC30 (Solute-Linked Carrier30, also named ZnT) and SLC39 (Solute-Linked Carrier 39, also named ZIP). ZnT contains 10 transporters of SLC30A1-SLC30A10, and ZIP contains 14 transporters of SLC39A1-SLC39A14. In our study, the differentially expressed SLC30A8 and SLC39A2 genes belong to ZnT family and ZIP family, respectively. Carbonic anhydrase located the eggshell gland epithelial cells is an important enzyme in the process of eggshell formation, which can reversibly catalyze the hydrolysis of H 2 CO 3 , regulate the concentration of HCO − in the eggshell gland, and then affect the Ca 2+ transport process and the calcium deposition in the eggshell, changing the quality of the eggshell [47]. Zinc ions are necessary for the activity center of carbonic anhydrase, so zinc can affect the activity of carbonic anhydrase [47]. Moreover, zinc is also a component of alkaline phosphatase, which may regulate some phosphorylated proteins related to the mechanism of eggshell formation and affect the synthesis of calcium carbonate crystals [48]. This provides us a vision for adding appropriate zinc to the diet of aging laying hens, which may reduce the deterioration of eggshell quality.
Through integration analysis of bioinformatics, we found that the differentially expressed TCONS_01464392 could target the GPX8 gene, whose expression was extremely signi cant, and their expression levels were negatively correlated. Glutathione peroxidases (GPXs) are enzymes that are present in almost all organisms, with the primary function of limiting peroxide accumulation. In mammals, GPXs consist of eight isoforms, but only two members (GPX7 and GPX8) reside in the endoplasmic reticulum [49,50]. A previous study demonstrates that GPX8 is enriched in mitochondria-associated membranes and can regulate Ca 2+ storage and uxes [49]. This indicates that the decline in eggshell quality of aging laying hens may be closely related to down-regulated GPX8 expression levels.

Conclusions
In conclusion, the present study provides a systematic genome-wide analysis of lncRNAs and mRNAs in the chicken shell gland, and the data obtained represent a valuable resource for further investigation of the function of some of these lncRNAs and mRNAs associated with eggshell quality deterioration during the late period of laying eggs.

Methods
Animal and Sample Collection 8 Hy-Line Brown commercial laying hens used in this study were purchased from Zhuozhou Chicken Farm. These hens were randomly assigned to old (60week-old, n = 4) and young (31-week-old, n = 4) groups. All birds included in this study were raised on the same diet and managed conditions until slaughtered.
18 hours after laying egg, animals were euthanized by exsanguination of the carotidartery under CO2 inhalation (approximately 5 min in small container gassed with CO2 from a compressed gas cylinder). Then, we collected the eggshell glands of each hen from the same pre-determined site and immediately ash frozen in liquid nitrogen and then stored at -80 °C until RNA extraction.

Total RNA Extraction
Total RNA was extracted from shell gland tissues using TRIzol reagent according to the manufacturer's instructions (Invitrogen, USA). RNA integrity was ascertained by 1.5% agarose gel electrophoresis, and the purity and concentration of the RNA were measured by spectrophotometer (ALLSHENG, China).
cDNA Library Construction and RNA Sequencing

Sequence Analysis Transcriptome Assembly
Quality control of the RNA-seq reads was performed using FastQC (http://www. bioinformatics.babraham.ac.uk/projects/fastqc/). Clean reads were obtained by removing empty reads, adapter sequences, reads with > 10% N sequences, and low quality reads, in which the number of bases with a quality value Q ≤ 10 was > 50%. At the same time, the Q30, GC content, and sequence duplication level of the clean data were calculated. Reads that passed the quality control were then mapped to the Gallus gallus reference genome (Ensembl v72, Galgal v4.0). Based on this, the mapped reads of each sample were assembled with StringTie (v1.3.1) [51] using a reference-based approach.

Screening and Prediction of DEGs and DE-lncRNAs
Fragments per kilo base of exon per million fragments mapped (FPKM), which was raised by Florea et al. [52], means the expected number of fragments per kilo base of transcript sequence per million reads sequenced. It takes into account the effects of sequencing depth and gene length on the fragment count and is currently the most commonly used method for estimating gene expression leve l [53]. In this study, transcript abundance was identi ed by FPKM using Cuffdiff (http://cu inks.cbcb.umd.edu/manual.html#cuffdiff) [13]. Here, FPKM was used to calculate the fold change of DEGs between the two groups, and the FPKM of the protein-coding genes in each sample was computed by summing the FPKMs of the transcripts in each gene group. Moreover, we analyzed DEGs by using the edgeR package to calculate the p-value that was obtained by multiple hypothesis testing calibration [54,55]. LncRNAs or protein-coding genes with p < 0.05 and log 2 (fold change) > 1 were assigned as DEGs.

Construction of the LncRNA-gene Interaction Network
Previous studies con rm that lncRNAs can regulate gene expression through cis-acting and trans-acting mechanisms [56]. For each lncRNA locus, the 10 k/100 k upstream and downstream protein-coding genes (without overlap) were rst identi ed as cis-target genes. However, the genes that overlapped with the lncRNAs predicted by Lnctar (http://www.cuilab.cn/lnctar) were selected as trans-target genes. To further investigate the interactions between the DE-lncRNAs and their corresponding differentially expressed cisor trans-target genes, we constructed an interactive lncRNA-gene network based on their FPKM using Cytoscape software (http://www.cytoscape.org). Moreover, we calculated the Pearson correlation coe cient (COR) of each lncRNA and DEG expression value.
GO and Pathway Analysis GO enrichment analysis of DEGs or lncRNA target genes was implemented using the Molecule Annotation System (MAS) 3.0 (http://bioinfo.capitalbio.com/mas3), which is based on the KEGG database (Capital Bio, Beijing). GO terms with p < 0.05 were considered signi cantly enriched by DEGs.
KEGG is a database resource for understanding high-level functions and utilities of a biological system [57], such as the cell, the organism, and the ecosystem, from molecular-level information, especially large-scale molecular datasets generated by genome sequencing and other high-throughput experimental technologies (http://www.genome.jp/kegg/). We used KOBAS software to test the statistical enrichment of DEGs or lncRNA target genes in KEGG pathways [58]. control. The qPCR primer sequences are presented in Table 1.

Statistical Analysis
The results of quantitative expression are presented as the mean ± standard error (SEM), and the signi cance of the data was tested by two-tailed paired Student's t-test using SPSS version 20.0 (SPSS Inc., Chicago, IL, USA). The 2 −ΔΔCt method was used to analyze the results of qPCR as described [59], and βactin was used as an internal control to normalize all of the threshold cycle (Ct) values.