Transcriptome Analysis Reveals the Potential Role of Long Non-coding RNAs in Mammary Gland of Yak During Lactation and Dry Period

The mammary gland is a remarkably dynamic organ of milk synthesis and secretion, and it experiences drastic structural and metabolic changes during the transition from dry periods to lactation, which involves the expression and regulation of numerous genes and regulatory factors. Long non-coding RNA (lncRNA) has considered as a novel type of regulatory factors involved in a variety of biological processes. However, their role in the lactation cycle of yak is still poorly understood. To reveal the involved mechanism, Ribo-zero RNA sequencing was employed to profile the lncRNA transcriptome in mammary tissue samples from yak at two physiological stages, namely lactation (LP) and dry period (DP). Notably, 1,599 lncRNA transcripts were identified through four rigorous steps and filtered through protein-coding ability. A total of 59 lncRNAs showed significantly different expression between two stages. Accordingly, the results of qRT-PCR were consistent with that of the transcriptome data. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses indicated that target genes of differentially expressed lncRNAs (DELs) were involved in pathways related to lactation, such as ECM-receptor interaction, PI3K-Akt signaling pathway, biosynthesis of amino acids and focal adhesion etc. Finally, we constructed a lncRNA-gene regulatory network containing some well known candidate genes for milk yield and quality traits. This is the first study to demonstrate a global profile of lncRNA expression in the mammary gland of yak. These results contribute to a valuable resource for future genetic and molecular studies on improving milk yield and quality, and help us to gain a better understanding of the molecular mechanisms underlying lactogenesis and mammary gland development of yak.

Introduction

In the past decade, human whole-genome transcriptional study discovered that approximately two-thirds of genomic DNA is pervasively transcribed (Djebali et al., 2012). Still less than 2% of the mammalian genome translates into the proteins, indicating that a considerable portion of the mammalian genome is transcribed into non-coding RNAs (ncRNAs) and do not function in protein-coding (Kopp and Mendell, 2018). The function-related non-protein-coding genome is particularly evident for a kind of small non-coding RNAs (ncRNAs) such as microRNAs (miRNAs), PIWI-interacting RNAs (piRNAs) and small interfering RNAs (siRNAs) (Pauli et al., 2011). Aside from small RNA, long non-coding RNAs (lncRNAs) have been widely concerned in recent years as potentially novel and pivotal regulators of biological process. LncRNAs are a heterogeneous group of non-protein-coding transcripts of greater than 200 nucleotides, which expression levels are usually lower than that of protein-coding genes (Hezroni et al., 2015; Ulitsky, 2016). Based on their position in the genome, lncRNAs can be categorized as long intergenic ncRNA (lincRNAs), enhancer RNAs (eRNAs), intronic lncRNAs, anti-sense lncRNAs and promoter-associated short RNAs (pasRNA) (Weikard et al., 2017). The regulatory roles of lncRNAs are gradually being revealed in a diverse spectrum of biological processes such as genomic imprinting (Sleutels et al., 2002), RNA splicing (Gonzalez et al., 2015), chromosome conformation (Rinn and Chang, 2012), epigenetic regulation (Vizoso and Esteller, 2012), transcriptional control (Salmena et al., 2011), and allosterically regulating enzymatic activity (Quinn and Chang, 2016).

With the increasing interest in lncRNA, a growing number of studies have discovered numerous lncRNAs in mammals, including Homo sapiens (Khalil et al., 2009), Mus musculus (Guttman et al., 2010), Bos taurus (Kern et al., 2018), Bubalus bubalis (Huang et al., 2019), and Sus scrofa (Gao et al., 2017). Recent studies have showed that several lncRNAs serve a role in the developing mammary gland and lactation. Pregnancy-induced non-coding RNA (PINC), a mammalian specific and alternatively spliced lncRNA earliest found in rat mammary gland, which has a function in regulating survival and cell cycle progression of mammary epithelial cells (Ginger et al., 2006). Increased expression of PINC in HC11 cells inhibits lactogenic differentiation, while suppressed expression of PINC could promote differentiation (Shore et al., 2012). LncRNA Zfas1 is localized in the epithelial cells of the mammary gland duct and alveoli. Loss of Zfas1 could promote the cell proliferation and differentiation in mammary epithelial cell line (Askarian-Amiri et al., 2011). Nuclear-enriched abundant transcript 1 (Neat1) has been reported to play a critical role in corpus luteum (Nakagawa et al., 2014) and placenta (Gremlich et al., 2014) development process. It also involves in mouse mammary gland morphogenesis and lactation. Knockout Neat1 in mouse alveolar cells could decrease the rates of cell proliferation (Standaert et al., 2014). LncRNA X inactivate−specific transcript (XIST) plays a vital role in the X-inactivation process in female mammals (Chureau et al., 2011). Previous observation demonstrated that XIST inhibits cell viability and promotes apoptosis in bovine mammary alveolar cell−T under inflammatory conditions (Ma et al., 2019).

Yak (Bos grunniens) is a multipurpose livestock, providing milk, meat, wool and transportation for indigenous people on the Qinghai-Tibetan Plateau. Yak milk and milk products are the main part of the daily diets of Tibetan nomads and play a key role in the health maintenance of Tibetans in the hypoxic environment. Compared with milk of dairy cattle, yak milk contains a higher content of protein, antioxidant vitamins, specific enzymes, and biologically active fatty acids (Ding et al., 2013; Guo et al., 2014). However, industrialized production of yak milk production is limited by the low milk yield of yak (150–500 kg of fresh milk per lactation) (Li et al., 2011). Therefore, improving milk yield is one of the significant breeding objectives in the yak industry. The mammary gland is a crucial secretory gland for milk synthesis and secretion, which provides necessary nutrients for human and neonatal offspring (McManaman and Neville, 2003; Tong et al., 2017). Post-pubertal mammary gland experiences a cycle of cell proliferation, differentiation, de-differentiation and apoptosis during its lactation process that is under the regulation of distinct hormones and regulating genes (Hennighausen and Robinson, 2005; Watson and Khaled, 2008). Detection of expression and function of coding genes and ncRNAs at different lactation stages of post-pubertal mammary gland development is not only essential for investigating the molecular mechanism of lactation, but also for improving milk yield and quality. A previous transcriptome study has identified numerous genes related to lactogenesis, milk secretion, and mammary gland development of yak during lactation and dry period (Fan et al., 2018). However, biological function of lncRNA have not been systematically identified in mammary glands of yak during the different lactation stages. Here, comparative analysis of lncRNAs in yak mammary gland was performed at two physiology stages (lactation period and dry period) by constructing five RNA sequencing (RNA-seq) libraries and sequencing. Our study will provide precious resources for yak lncRNA studies and contribute to a better understanding of the regulatory mechanisms of lactation.

Materials and Methods

Animal and Sample Collection

Five healthy and mastitis-free female Ashidan yaks in their second parity were selected for this study. All individuals were freely grazed on natural pasture with free drinking at Datong yak farm of Qinghai Province. Due to the limitation of sampling condition, we only selected two yaks at peak lactation (LP; 120 days postpartum; non-pregnant period), and three yaks at the dry period (DP; calved in the previous year, non-pregnant period). After slaughter, the rear mammary gland of each yak was cut immediately, and alveolar tissue from the middle of the upper one-third of the gland was harvested and stored in liquid nitrogen. All of the experimental protocols and procedures were approved by Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences of CAAS (Permit No. SYXK-2014-0002).

Total RNA Isolation

Total RNA was isolated from mammary gland tissues with the TRIzol reagent following the manufacturer’s protocol (Invitrogen, CA, United States). The quality and purity of total RNA were evaluated by a NanoPhotometer^® spectrophotometer (IMPLEN, CA, United States), and integrity was assessed by the RNA Nano 6000 Assay Kit of the Bioanalyzer 2100 system (Agilent Technologies, CA, United States). The 260/280 ratio for each sample was approximately 2.0, and the RNA integrity number (RIN) was greater than 8.0.

Library Preparation and Sequencing

The same amount of total RNA (3 μg) from each sample was used to construct RNA-Sequencing library, and ribosomal RNA was depleted from total RNA by Epicentre Ribo-zero^TM rRNA Removal Kit (Epicentre, United States). Subsequently, sequencing libraries were separately generated using the rRNA-depleted RNA by NEBNext^® Ultra^TM Directional RNA Library Prep Kit for Illumina^® (NEB, United States) following the manufacturer’s instructions. The quality of all libraries was assessed on the Agilent Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, CA, United States). Sequencing libraries were then sequenced on an Illumina Hiseq 2500 platform (Novogene Bioinformatics Institute, Beijing, China) to generate 125 bp paired-end reads. Quality control of RNA-seq reads was performed using FastQC (Andrews, 2012).

Data Analysis

Raw data (raw reads) of fastq format were first filtered through in-house Perl scripts. In this step, clean reads were obtained by removal of the reads containing adapter molecules, reads containing unknown bases (>10%), and low quality reads (more than 50% of the bases had quality scores ≤10) from raw reads. Meanwhile, the Q30 and GC contents of the clean data were calculated. Genome sequence of yak (BosGru_v2.0) and annotation file were downloaded from NCBI¹. Paired-end clean reads were mapped to the yak genome using Tophat2 (v2.0.9) (Kim et al., 2013). Scripture (beta2 version) (Guttman et al., 2010) and Cufflinks (v2.1.1) software (Trapnell et al., 2010) were used to assemble the transcripts. Cuffdiff (version 2.1.1) (Trapnell et al., 2010) was used to calculate thefragments per kilobase of exon per million fragments mapped (FPKM) value of transcripts with default parameters, including lncRNAs and mRNAs in each sample.

LncRNAs Identification

To identify the putative lncRNA, we filtered the assembled transcripts through the following exclusion criterion (Figure 1A). (1) transcripts with a single exon were removed. (2) transcripts with a length of more than 200 nt were kept. (3) transcripts with FPKM ≥ 0.5 were kept. (4) transcripts that overlap with the exon region of the database annotation were removed. (5) CNCI (Coding-Non-Coding-Index, v2) (Sun et al., 2013), CPC (encoding potential calculator, 0.9-r2) (Kong et al., 2007) and Pfam Scan (v1.3) (Punta et al., 2012) were used to predict the protein-coding ability of the transcripts. Transcripts without protein-coding ability were selected as a candidate set of lncRNAs. The expression levels of mRNA and lncRNA were calculated based on assembled transcript files. The edgeR (exact test for negative binomial distribution) Bioconductor package (Robinson et al., 2010) in R software was used to identify differentially expressed lncRNAs (DELs) and differentially expressed genes (DEGs) between LP and DP groups. For biological replicates, genes or lncRNAs with | log2 (fold change)| ≥ 1 and P-value < 0.05 were considered as differentially expressed according to the previous study (Zheng et al., 2018).

FIGURE 1

Figure 1. Generation and characterization of lncRNAs in mammary gland of yak. (A) Pipeline for identification of lncRNAs. (B) Classification of lncRNAs. (C) The length distribution of lncRNAs. (D) ORF length distribution of lncRNAs. (E) Exon number distribution of lncRNAs. (F) Expression level of lncRNAs and mRNAs. Red and blue represent lncRNAs and mRNAs, respectively.

LncRNA Target Gene Prediction and Functional Enrichment Analysis

The cis- and trans-acting of lncRNAs can modulate the interactions with target genes. The cis role of lncRNAs is to interact with neighboring target genes. In this study, the protein-coding genes located in 100 kb upstream and downstream of DELs were searched as described previously (Gao et al., 2019). The trans role refers to lncRNAs that regulate other genes on different positions of the genome. To identify the trans-target genes of the DELs, the expressed correlation of DELs and DEGs was estimated with Pearson’s correlation coefficient (r) using the cor.test function in R software. Those DEGs with a P value <0.05 and | r-value| > 0.95 were considered as trans-target genes. To infer the function of the lncRNA, we used the G:Profile web tool (Raudvere et al., 2019) to implement gene ontology (GO) analyses of target genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of lncRNA target genes was conducted using KOBAS 3.0 software (Xie et al., 2011) by using a hypergeometric test. Pathway with gene count >5 and P-value <0.05 were considered signifcantly enriched by target genes. Based on the enrichment analysis results, function and quantitative trait locus (QTL) information of the differentially expressed cis- and trans-target genes, DELs and target genes associated with milk traits and lactation were selected to construct the regulatory network. We used Cytoscape v.3.7.2 software (Cline et al., 2007) to visualize the regulatory networks.

Verification of Sequencing Data Using qRT-PCR

To verify the accuracy of sequencing results, quantitative real-time PCR (qRT-PCR) method was used to detect the relative expression of ten randomly selected DELs. One microgram of total RNA from each sample was transcribed into cDNA using Transcriptor First Strand cDNA Synthesis Kit (Roche Diagnostics, Mannheim, Germany) following the manufacturer’s protocols. The lncRNAs expressions were then measured by qRT-PCR using the LightCycler 96 Real-Time PCR system (Roche Diagnostics, Mannheim, Germany). The 20 μL PCR reactions included 1 μL of diluted cDNA, 1 μL of forward primer (10 μM), 1 μL of reverse primer (10 μM), 10 μL of SYBR TB Green mix (TaKaRa, Dalian, China), and 7 μL of ddH₂O. The PCR program for amplification was as follows: an initial denaturation step of 30 s at 95°C, and 45 cycles of 5 s at 95°C, and 30 s at 60°C. A final melting program ranging from 65°C to 97°C with an increments of 0.5°C and acquiring fluorescence after each step. The relative expression of DELs was analyzed using the 2^–ΔΔCt method (Silver et al., 2006) and normalized by hydroxymethylbilane synthase (HMBS) and tyrosine 3-monooxygenase/tryptophan 5-monooxygenase activation protein, zeta polypeptide (YWHAZ) (Wu et al., 2019).

Results

Overview of RNA-seq Data

In order to identify novel lncRNAs in mammary gland of yak, five cDNA libraries were sequenced using Illumina Hiseq 2500 platform, and 503.07 and 299.23 million raw paired-end reads were obtained in lactation and dry stages, respectively. After removing the adaptor sequence and low-quality raw reads, we obtained 784.05 million clean reads with an average of 94.35 million reads (ranging from 135.20 to 181.06) for each sample. The GC contents were 50.51–53.23% while the Q30 ranged between 96.35 and 96.53%. Of the clean reads, 83.79–86.48% could be perfectly mapped to the yak reference genome by the Tophat2 software (Kim et al., 2013) (Table 1). Pearson’s correlation analysis showed that the correlation between the groups was strong and providing reliable data for the next steps (Supplementary Figure 1).

TABLE 1

Table 1. Summary of RNA-Seq data and mapping.

Characterization and Expression Profiling of LncRNAs in the Yak Mammary Gland

We identified 1,599 novel lncRNAs from the yak mammary gland tissues by using CNCI, CPC and Pfam software (Supplementary Table 1 and Figure 2). According to the genome alignment, these lncRNA transcripts were classified into 1,374 (85.9%) lincRNA and 225 (14.1%) anti-sense lncRNA (Figure 1B). The size of lncRNAs length in the mammary gland ranged from 219 to 26,695 nt, with 63.7% of the total number of lncRNAs were approximately 200 to 2,000 nt in length (Figure 1C). The lncRNAs whose ORF length ranged from 101 to 150 bp occupied the largest proportion of the total lncRNAs (Figure 1D). Most of the lncRNAs contained three or fewer exons (84.6%) and the average exon number of the lncRNA was 2.65 (Figure 1E). Expression levels of the lncRNAs were also lower than those of the genes (Figure 1F).

FIGURE 2

Figure 2. Number of upregulated and downregulated DEGs and DELs of LP vs DP.

Expression levels of the lncRNA transcripts were determined through FPKM method by using Cufflinks (v2.1.1) software. A total of 59 DELs were identified transcripts between LP and DP groups (| log2 (fold change)| ≥ 1 and P-value < 0.05). Of these DELs, 25 lncRNAs were significantly upregulated and 34 lncRNAs were downregulated in the mammary gland at lactation stage (Figure 2 and Supplementary Table 2). We also found 833 DEG between LP and DP groups | log2 (fold change)| ≥ 1 and P-value < 0.05) (Figure 2 and Supplementary Table 3).

LncRNAs Target Genes of Cis and Trans Regulated

To further explore the function of lncRNAs in yak mammary gland, we identified potential cis- and trans-target genes of DELs, and predicted the function of target genes through GO and KEGG pathways enrichment analysis. The cis-targets of lncRNAs were predicted by screening genes around 100 kb upstream and downstream of lncRNAs. The results showed that there were 42 lncRNAs correspond to 130 cis-target genes (Supplementary Table 4). GO enrichment analysis revealed cis-target genes of DELs were significantly enriched in 182 terms (131 under biological processes, 15 under cellular compartments, and 36 under molecular functions) (adjusted P value <0.05). Detailed information is listed in Supplementary Table 5. The most significantly enriched GO terms were cellular anatomical entity, cellular anatomical entity, binding, biological regulation and intracellular. However, there was no significant enrichment of pathway for cis-target genes of DELs.

We also predicted potential target genes of lncRNAs in a trans manner. A total of 3,228 interactions (including 2,549 positive and 679 negative correlations) were identified between 59 DELs and 758 DEGs (Supplementary Table 6). In our results, 23 cis-genes overlapped with trans-target genes. GO analysis showed that the trans-target genes were significantly enriched in 836 terms (654 biological processes, 88 cellular compartments, and 94 molecular functions), which involves a variety of biological processes (Supplementary Table 7). The top enriched terms were cellular process, cellular anatomical entity and bingding (Figure 3A). It is important that some of the terms were related to mammary gland development, such as cell surface receptor signaling pathway, cell differentiation, cell population proliferation, epithelial cell proliferation, and mammary gland development. KEGG analysis revealed the trans-target genes were enriched in 35 pathways (Supplementary Table 8). The top 30 KEGG pathways are showed in Figure 3B. Importantly, a few pathways were associated with lactation, such as ECM-receptor interaction, PI3K-Akt signaling pathway, PPAR signaling pathway, biosynthesis of amino acids, suggesting that lncRNAs act in the trans-regulation of genes related to lactation. Moreover, multiple DELs had a common target gene.

FIGURE 3

Figure 3. GO and KEGG enrichment of trans-target genes of DELs. (A) Top30 GO enrichment terms. (B) Top30 KEGG pathway enrichment terms. The horizontal axis represents a rich factor and the vertical axis represents the pathway. The size of the bubble indicates the number of target genes enriched in the pathway, and the color of the bubble represents a different P-value range. Rich Factor is the ratio of differentially expressed gene numbers annotated in this pathway terms to all gene numbers annotated in this pathway term.

Integrated Analysis

To better understand the relationship between DELs and lactation process, we chose 42 target genes related to lactation and mammary gland development to construct lncRNA-gene regulatory network. The regulatory network showed that 54 DELs interacted with at one target genes associated with lactation and mammary gland development (Figure 4 and Supplementary Table 9). Notably, multiple lncRNAs (LNC_001595, LNC_000894 and LNC_001600 and others) appear to regulate insulin like growth factor 1 (IGF1), indicating that they might regulate lactation process though the PI3K-Akt pathway.

FIGURE 4

Figure 4. Network plot of candidate lncRNAs and mRNAs. The orange triangles, blue triangles, orange round and blue round represent upregulated lncRNAs, downregulated lncRNAs, upregulated mRNAs and downregulated mRNAs, respectively.

Verification of Gene Expression Profiles Using qRT-PCR

To validate the reproducibility of the DELs obtained from RNA-seq, ten DELs were randomly selected for confirmation by qRT-PCR. The primer details are presented in Supplementary Table 10. As showed in Figure 5, the expressions of LNC_001445, LNC_001581, LNC_000118 and LNC_000472 were upregulated, whereas LNC_001571, LNC_000408, LNC_001114, LNC_001248, LNC_001049 and LNC_001556 were downregulated in LP compared to DP. All ten DELs had similar expression patterns in comparison to the RNA-seq data, indicating the reliability of our RNA-seq data.

FIGURE 5

Figure 5. Validation and comparison of log2(fold change) in ten DELs between qRT-PCR and RNA-Seq.

Discussion

The mammary gland is a complex organ for synthesis and secretion of milk, which undergoes functional differentiation during the female reproductive cycle (Hurley, 1989). During pregnancy and lactation, epithelial cells of the mammary gland start to proliferate and differentiate into functional alveolar cells, which produce and secrete milk after parturition (Macias and Hinck, 2012). In the dry period, apoptosis of the majority of differentiated alveolar cells was induced, and gland restores to a pre-pregnancy state (Watson, 2006). Deciphering these physiological changes is essential to explore the complex molecular mechanisms underlying the transition from the dry period to lactation. One of the important aspects was the transcript originated from the non-coding regions of the genome. However, systematic prediction and regulatory function of lncRNAs in yak mammary gland have rarely been reported. In this study, we comprehensively investigated lncRNAs in the mammary gland of yak at lactation and dry period.

In the present study, more than 83.79% of the clean reads were mapped to the yak genome, which was similar to those of other RNA-seq studies (Zi et al., 2018). Similar to lncRNA features of previous studies, the average spanning exon number and expression levels of lncRNA was much less than that of the genes, indicating the reliability of lncRNA identification (Zheng et al., 2018; Li et al., 2020). Furthermore, we identified a total of 59 DELs from 5 samples, comprising 25 upregulated and 34 downregulated DELs in LP compared with the DP (| fold change| > 2, P-value <0.05). Remarkably, in our study design, only two biological replicates in LP were used due to the limited sample availability. Rapaport et al. (2013) demonstrated that over 90% of differently expressed genes at the top expression levels could be detected with using two replicates and 5% of the reads. However, higher biological replicate numbers are still recommended to improve the detection power.

Recent studies have demonstrated that lncRNAs control gene expression through both cis- and trans-acting mechanism, which play an important role in a wide range of biological processes (Bhat et al., 2016). From this view, the function of lncRNA can be inferred on the basis of knowledge of their target genes. Cis-regulatory lncRNAs have enhancer-like activity and promote expression of neighboring genes. Accordingly, the cis prediction analysis showed that there were 42 DELs near to 130 genes with less than 100 kb distance. It is worth noting that some of the cis-target genes were involved in lactation and mammary gland development. For example, AKT Serine/Threonine Kinase 1 (AKT1) was predicted to be a target gene of LNC_000763, which regulates of milk synthesis and metabolism in the lactating mammary gland (Boxer et al., 2006). Constitutively active Akt1 in mice delayed mammary gland involution and epithelial cell apoptosis (Ackler et al., 2002). In contrast, depletion of Akt1 in mice resulted in a few mammary gland developmental defects, such as ductal outgrowth and deficiency in terminal end bud growth and alveolar bud development (LaRocca et al., 2011). LIM homeobox 3(LHX3) is a predicted cis-target of LNC_000108, which regulates the release of hormones related to lactation and metabolism (Pfäffle and Klammt, 2011). Liu et al. (2011) found that variations of LHX3 gene were significantly associated with milk performance in goat. Secreted frizzled-related protein-2 (SFRP2) was predicted to be a target gene of LNC_001049. A previous study indicated that the interaction of SFRP2 with extracellular matrix (ECM) could prevent apoptosis of mammary epithelial cells (Lee et al., 2004). LNC_000118 was predicted to act on the target gene activating transcription factor 3 (ATF3). It is reported that the expression of ATF3 was significantly different in the mammary glands of Holstein cows with extremely high and low content of milk protein and fat (Cui et al., 2014). During the lactation cycle, the expression of ATF3 increased to a peak at day 15, followed by a decrease through day 240 relative to parturition in bovine mammary tissue (Invernizzi et al., 2012). These findings suggest that lncRNA might regulate lactation, milk synthesis and mammary gland development through their actions on neighboring genes.

Still, many lncRNAs act in trans mode to regulate the target genes, which are distant from the transcription sites of the lncRNAs. Co-expression analysis identified 59 DELs were interacted with 758 genes based on the expression correlation coefficient (r-value > 0.95 or < −0.95 and P-value <0.05). The results of KEGG analysis of the trans-target genes were involved in several keys signaling pathways, such as ECM-receptor interaction, PI3K-Akt signaling pathway, biosynthesis of amino acids and focal adhesion. These above KEGG pathways have been suggested to participate in the lactation process and mammary gland development (Watson, 2006; Nagy et al., 2007; Kim and Wu, 2009; Macias and Hinck, 2012). PPARG Coactivator 1 Alpha (PPARGC1A) is predicted to be a target gene of LNC_000883, LNC_001601 and LNC_000771. PPARGC1A, also known as PGC-1a, is a transcriptional coactivator that regulates genes involved in eglucose and lipid transportation and oxidation, which plays a crucial role in mammary gland development and lactation (Qiu et al., 2020). The polymorphism of the PPARGC1A gene has been proved to be significantly associated with milk fat yield, milk fat composition, milk protein percentage (Weikard et al., 2005; Khatib et al., 2007; Schennink et al., 2009). Ten DELs (such as LNC_001595, LNC_000894 and LNC_001600) were co-expressed with insulin like growth factor 1 (IGF1), a master regulator of mammary gland development, which considered as a candidate gene for milk protein and fat yields in Polish Holstein Friesian cows (Szewczuk et al., 2011). It is also reported that overexpression of IGF1 gene in mice during pregnancy and lactation caused a delay of the apoptotic loss of mammary epithelial cells during the declining phase of lactation (Hadsell et al., 2002). In addition, 6 DELs were predicted to act on the target gene beta1 integrin (ITGB1), which plays a crucial role in the maintenance of functional mammary stem cells and mammary morphogenesis (Taddei et al., 2008). In the mammary luminal cell population, conditional deletion of ITGB1 could inhibit the alveologenesis and lactation (Naylor et al., 2005). Furthermore, our results showed that a cluster of DELs regulates well-known genes affecting milk traits; e.g., lactotransferrin (LTF) (El-Domany et al., 2019) and lactalbumin alpha (LALBA) (García-Gámez et al., 2012). These observations suggested that these lncRNAs could be closely related to milk secretion due to their trans-targeted genes affecting lactation and mammary gland development. Although these lncRNAs need further experimental investigation, this information may help us to explore the potential regulatory mechanisms of lncRNA involved in lactation and mammary gland development in the yak.

Conclusion

In the current study, we analyzed lncRNA expression profiles in the mammary gland of yak at two different physiological stages (LP and DP) for the first time. A total of 59 DELs were identified between two stages. Furthermore, several DELs related to lactation and mammary gland were annotated based on GO and KEGG analysis of cis- and trans-target genes. The DELs in the regulatory network may be considered as promising targets for the exploration of functional lncRNAs in yak and thus could be used to improve the milk traits. Considering the fact that the lncRNAs annotated on yaks are not complete, this work provides a valuable resource for improving the lncRNA database of yak and ideal candidates for future studies to illustrate the molecular mechanisms linking lncRNA to the regulation of lactation and mammary gland development.

Data Availability Statement

The datasets generated for this study can be found in the Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) at NCBI, with the BioProject ID: PRJNA626061.

Ethics Statement

The animal study was reviewed and approved by Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences of CAAS (Permit No. SYXK-2014-0002).

Author Contributions

XW, PY, and XG conceived the research and drafted the manuscript with comments from all authors. XW, XZ, JP, LX, and XY conducted the experiments and analyses. CL, PB, and MC participated in experiments and revised the manuscript. All authors approved the final version.

Funding

This study was supported by the Agricultural Science and Technology Innovation Program (CAAS-ASTIP-2014-LIHPS-01), National Beef Cattle Industry Technology & System (CARS-37), and the Special Fund of the Chinese Central Government for Basic Scientific Research Operations in Commonweal Research Institutes (No. 1610322020001).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2020.579708/full#supplementary-material

Supplementary Figure 1 | Pearson’s correlation analysis between samples.

Supplementary Figure 2 | Candidate lncRNA identification with coding potential assessment of CNCI, CPC2, and PfamScan.

Footnotes

1. ^ http://www.ncbi.nlm.nih.gov/

References

Ackler, S., Ahmad, S., Tobias, C., Johnson, M. D., and Glazer, R. I. (2002). Delayed mammary gland involution in MMTV-AKT1 transgenic mice. Oncogene 21, 198–206. doi: 10.1038/sj.onc.1205052

PubMed Abstract | CrossRef Full Text | Google Scholar

Andrews, S. (2012). FastQC: A Quality Control Application for High Throughput Sequence Data. Available online at: https://www.bioinformatics.babraham.ac.uk/projects/fastqc/ (accessed October 07, 2019).

Google Scholar

Askarian-Amiri, M. E., Crawford, J., French, J. D., Smart, C. E., Smith, M. A., Clark, M. B., et al. (2011). SNORD-host RNA Zfas1 is a regulator of mammary development and a potential marker for breast cancer. RNA 17, 878–891. doi: 10.1261/rna.2528811

PubMed Abstract | CrossRef Full Text | Google Scholar

Bhat, S. A., Ahmad, S. M., Mumtaz, P. T., Malik, A. A., Dar, M. A., Urwat, U., et al. (2016). Long non-coding RNAs: mechanism of action and functional utility. Noncod. RNA Res. 1, 43–50. doi: 10.1016/j.ncrna.2016.11.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Boxer, R. B., Stairs, D. B., Dugan, K. D., Notarfrancesco, K. L., Portocarrero, C. P., Keister, B. A., et al. (2006). Isoform-specific requirement for Akt1 in the developmental regulation of cellular metabolism during lactation. Cell Metab. 4, 475–490. doi: 10.1016/j.cmet.2006.10.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Chureau, C., Chantalat, S., Romito, A., Galvani, A., Duret, L., Avner, P., et al. (2011). Ftx is a non-coding RNA which affects Xist expression and chromatin structure within the X-inactivation center region. Hum. Mol. Genet. 20, 705–718. doi: 10.1093/hmg/ddq516

PubMed Abstract | CrossRef Full Text | Google Scholar

Cline, M. S., Smoot, M., Cerami, E., Kuchinsky, A., Landys, N., Workman, C., et al. (2007). Integration of biological networks and gene expression data using Cytoscape. Nat. Protoc. 2, 2366–2382. doi: 10.1038/nprot.2007.324

PubMed Abstract | CrossRef Full Text | Google Scholar

Cui, X., Hou, Y., Yang, S., Xie, Y., Zhang, S., Zhang, Y., et al. (2014). Transcriptional profiling of mammary gland in Holstein cows with extremely different milk protein and fat percentage using RNA sequencing. BMC Genomics 15:226. doi: 10.1186/1471-2164-15-226

PubMed Abstract | CrossRef Full Text | Google Scholar

Ding, L., Wang, Y., Kreuzer, M., Guo, X., Mi, J., Gou, Y., et al. (2013). Seasonal variations in the fatty acid profile of milk from yaks grazing on the Qinghai-Tibetan plateau. J. Dairy Res. 80, 410–417. doi: 10.1017/S0022029913000496

PubMed Abstract | CrossRef Full Text | Google Scholar

Djebali, S., Davis, C. A., Merkel, A., Dobin, A., Lassmann, T., Mortazavi, A., et al. (2012). Landscape of transcription in human cells. Nature 489, 101–108. doi: 10.1038/nature11233

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Domany, W. B., Radwan, H. A., Ateya, A. I., Ramadan, H. H., Marghani, B. H., and Nasr, S. M. (2019). Genetic Polymorphisms in LTF/EcoRI and TLR4/AluI loci as candidates for milk and reproductive performance assessment in Holstein cattle. Reprod. Domest. Anim. 54, 678–686. doi: 10.1111/rda.13408

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, J. F., Luo, Y. Z., Yu, S. J., Cui, Y., Xu, G. Q., Wang, L. B., et al. (2018). Transcriptional profiling of two different physiological states of the yak mammary gland using RNA sequencing. PLoS One 13:e201628. doi: 10.1371/journal.pone.0201628

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, P. F., Guo, X. H., Du, M., Cao, G. Q., Yang, Q. C., Pu, Z. D., et al. (2017). LncRNA profiling of skeletal muscles in large white pigs and mashen pigs during development. J. Anim. Sci. 95, 4239–4250. doi: 10.2527/jas2016.1297

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Li, S., Lai, Z., Zhou, Z., Wu, F., Huang, Y., et al. (2019). Analysis of long non-coding RNA and mRNA expression profiling in immature and mature bovine (Bos taurus) testes. Front. Genet. 10:646. doi: 10.3389/fgene.2019.00646

PubMed Abstract | CrossRef Full Text | Google Scholar

García-Gámez, E., Gutiérrez-Gil, B., Sahana, G., Sánchez, J. P., Bayón, Y., and Arranz, J. J. (2012). GWA analysis for milk production traits in dairy sheep and genetic support for a QTN influencing milk protein percentage in the LALBA gene. PLoS One 7:e47782. doi: 10.1371/journal.pone.0047782

PubMed Abstract | CrossRef Full Text | Google Scholar

Ginger, M. R., Shore, A. N., Contreras, A., Rijnkels, M., Miller, J., Gonzalez-Rimbau, M. F., et al. (2006). A noncoding RNA is a potential marker of cell fate during mammary gland development. Proc. Natl. Acad. Sci. U.S.A. 103, 5781–5786. doi: 10.1073/pnas.0600745103

PubMed Abstract | CrossRef Full Text | Google Scholar

Gonzalez, I., Munita, R., Agirre, E., Dittmer, T. A., Gysling, K., Misteli, T., et al. (2015). A lncRNA regulates alternative splicing via establishment of a splicing-specific chromatin signature. Nat. Struct. Mol. Biol. 22, 370–376. doi: 10.1038/nsmb.3005

PubMed Abstract | CrossRef Full Text | Google Scholar

Gremlich, S., Damnon, F., Reymondin, D., Braissant, O., Schittny, J. C., Baud, D., et al. (2014). The long non-coding RNA NEAT1 is increased in IUGR placentas, leading to potential new hypotheses of IUGR origin/development. Placenta 35, 44–49. doi: 10.1016/j.placenta.2013.11.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Guo, X., Long, R., Kreuzer, M., Ding, L., Shang, Z., Zhang, Y., et al. (2014). Importance of functional ingredients in yak milk-derived food on health of Tibetan nomads living under high-altitude stress: a review. Crit. Rev. Food Sci. 54, 292–302. doi: 10.1080/10408398.2011.584134

PubMed Abstract | CrossRef Full Text | Google Scholar

Guttman, M., Garber, M., Levin, J. Z., Donaghey, J., Robinson, J., Adiconis, X., et al. (2010). Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs. Nat. Biotechnol. 28, 503–510. doi: 10.1038/nbt.1633

PubMed Abstract | CrossRef Full Text | Google Scholar

Hadsell, D. L., Bonnette, S. G., and Lee, A. V. (2002). Genetic manipulation of the IGF-I axis to regulate mammary gland development and function. J. Dairy Sci. 85, 365–377. doi: 10.3168/jds.S0022-0302(02)74083-6

CrossRef Full Text | Google Scholar

Hennighausen, L., and Robinson, G. W. (2005). Information networks in the mammary gland. Nat. Rev. Mol. Cell Biol. 6, 715–725. doi: 10.1038/nrm1714

PubMed Abstract | CrossRef Full Text | Google Scholar

Hezroni, H., Koppstein, D., Schwartz, M. G., Avrutin, A., Bartel, D. P., and Ulitsky, I. (2015). Principles of long noncoding RNA evolution derived from direct comparison of transcriptomes in 17 species. Cell Rep. 11, 1110–1122. doi: 10.1016/j.celrep.2015.04.023

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J., Zheng, Q., Wang, S., Wei, X., Li, F., and Ma, Y. (2019). High-throughput RNA sequencing reveals NDUFC2-AS lncRNA promotes adipogenic differentiation in Chinese buffalo (Bubalus bubalis L). Genes 10:689. doi: 10.3390/genes10090689

PubMed Abstract | CrossRef Full Text | Google Scholar

Hurley, W. L. (1989). Mammary gland function during involution. J. Dairy Sci. 72, 1637–1646. doi: 10.3168/jds.S0022-0302(89)79276-6

CrossRef Full Text | Google Scholar

Invernizzi, G., Naeem, A., and Loor, J. J. (2012). Endoplasmic reticulum stress gene network expression in bovine mammary tissue during the lactation cycle. J. Dairy Sci. 95, 2562–2566. doi: 10.3168/jds.2011-4806

PubMed Abstract | CrossRef Full Text | Google Scholar

Kern, C., Wang, Y., Chitwood, J., Korf, I., Delany, M., Cheng, H., et al. (2018). Genome-wide identification of tissue-specific long non-coding RNA in three farm animal species. BMC Genomics 19:684. doi: 10.1186/s12864-018-5037-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Khalil, A. M., Guttman, M., Huarte, M., Garber, M., Raj, A., Rivea, M. D., et al. (2009). Many human large intergenic noncoding RNAs associate with chromatin-modifying complexes and affect gene expression. Proc. Natl. Acad. Sci. U.S.A. 106, 11667–11672. doi: 10.1073/pnas.0904715106

PubMed Abstract | CrossRef Full Text | Google Scholar

Khatib, H., Zaitoun, I., Wiebelhaus-Finger, J., Chang, Y. M., and Rosa, G. J. (2007). The association of bovine PPARGC1A and OPN genes with milk composition in two independent Holstein cattle populations. J. Dairy Sci. 90, 2966–2970. doi: 10.3168/jds.2006-812

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14:R36. doi: 10.1186/gb-2013-14-4-r36

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, S. W., and Wu, G. (2009). Regulatory role for amino acids in mammary gland growth and milk synthesis. Amino Acids 37, 89–95. doi: 10.1007/s00726-008-0151-155

CrossRef Full Text | Google Scholar

Kong, L., Zhang, Y., Ye, Z. Q., Liu, X. Q., Zhao, S. Q., Wei, L., et al. (2007). CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 35, W345–W349. doi: 10.1093/nar/gkm391

PubMed Abstract | CrossRef Full Text | Google Scholar

Kopp, F., and Mendell, J. T. (2018). Functional classification and experimental dissection of long noncoding RNAs. Cell 172, 393–407. doi: 10.1016/j.cell.2018.01.011

PubMed Abstract | CrossRef Full Text | Google Scholar

LaRocca, J., Pietruska, J., and Hixon, M. (2011). Akt1 is essential for postnatal mammary gland development, function, and the expression of Btn1a1. PLoS One 6:e24432. doi: 10.1371/journal.pone.0024432

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, J. L., Lin, C. T., Chueh, L. L., and Chang, C. J. (2004). Autocrine/paracrine secreted frizzled-related protein 2 induces cellular resistance to apoptosis: a possible mechanism of mammary tumorigenesis. J. Biol. Chem. 279, 14602–14609. doi: 10.1074/jbc.M309008200

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Ma, Y., Li, Q., Wang, J., Cheng, J., Xue, J., et al. (2011). The chemical composition and nitrogen distribution of Chinese yak (Maiwa) milk. Int. J. Mol. Sci. 12, 4885–4895. doi: 10.3390/ijms12084885

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Li, C., Wureli, H., Ni, W., Zhang, M., Li, H., et al. (2020). Screening and evaluating of long non-coding RNAs in prenatal and postnatal pituitary gland of sheep. Genomics 112, 934–942. doi: 10.1016/j.ygeno.2019.06.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J. B., Lan, X. Y., Xu, Y., Li, Z. J., Lei, C. Z., and Chen, H. (2011). Combined effects of three novel SNPs within goat LHX3 gene on milk performance. Genes Genom. 33, 549–556. doi: 10.1007/s13258-011-0026-0

CrossRef Full Text | Google Scholar

Ma, M., Pei, Y., Wang, X., Feng, J., Zhang, Y., and Gao, M. Q. (2019). LncRNA XIST mediates bovine mammary epithelial cell inflammatory response via NF-kappaB/NLRP3 inflammasome pathway. Cell Prolif. 52:e12525. doi: 10.1111/cpr.12525

PubMed Abstract | CrossRef Full Text | Google Scholar

Macias, H., and Hinck, L. (2012). Mammary gland development. Wiley Interdiscip. Rev. Dev. Biol. 1, 533–557. doi: 10.1002/wdev.35

PubMed Abstract | CrossRef Full Text | Google Scholar

McManaman, J. L., and Neville, M. C. (2003). Mammary physiology and milk secretion. Adv. Drug Deliver. Rev. 55, 629–641. doi: 10.1016/S0169-409X(03)00033-4

CrossRef Full Text | Google Scholar

Nagy, T., Wei, H., Shen, T. L., Peng, X., Liang, C. C., Gan, B., et al. (2007). Mammary epithelial-specific deletion of the focal adhesion kinase gene leads to severe lobulo-alveolar hypoplasia and secretory immaturity of the murine mammary gland. J. Biol. Chem. 282, 31766–31776. doi: 10.1074/jbc.M705403200

PubMed Abstract | CrossRef Full Text | Google Scholar

Nakagawa, S., Shimada, M., Yanaka, K., Mito, M., Arai, T., Takahashi, E., et al. (2014). The lncRNA Neat1 is required for corpus luteum formation and the establishment of pregnancy in a subpopulation of mice. Development 141, 4618–4627. doi: 10.1242/dev.110544

PubMed Abstract | CrossRef Full Text | Google Scholar

Naylor, M. J., Li, N., Cheung, J., Lowe, E. T., Lambert, E., Marlow, R., et al. (2005). Ablation of beta1 integrin in mammary epithelium reveals a key role for integrin in glandular morphogenesis and differentiation. J. Cell Biol. 171, 717–728. doi: 10.1083/jcb.200503144

PubMed Abstract | CrossRef Full Text | Google Scholar

Pauli, A., Rinn, J. L., and Schier, A. F. (2011). Non-coding RNAs as regulators of embryogenesis. Nat. Rev. Genet. 12, 136–149. doi: 10.1038/nrg2904

PubMed Abstract | CrossRef Full Text | Google Scholar

Pfäffle, R., and Klammt, J. (2011). Pituitary transcription factors in the aetiology of combined pituitary hormone deficiency. Best Pract. Res. Clin. Endocrinol. Metab. 25, 43–60. doi: 10.1016/j.beem.2010.10.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Punta, M., Coggill, P. C., Eberhardt, R. Y., Mistry, J., Tate, J., Boursnell, C., et al. (2012). The Pfam protein families database. Nucleic Acids Res. 40, D290–D301. doi: 10.1093/nar/gkr1065

PubMed Abstract | CrossRef Full Text | Google Scholar

Qiu, L., Fan, X., Zhang, Y., Teng, X., and Miao, Y. (2020). Molecular characterization, tissue expression and polymorphisms of buffalo PPARGC1A gene. Arch. Tierzucht. 63, 249–259. doi: 10.5194/aab-63-249-2020

PubMed Abstract | CrossRef Full Text | Google Scholar

Quinn, J. J., and Chang, H. Y. (2016). Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet. 17, 47–62. doi: 10.1038/nrg.2015.10

PubMed Abstract | CrossRef Full Text | Google Scholar

Rapaport, F., Khanin, R., Liang, Y., Pirun, M., Krek, A., Zumbo, P., et al. (2013). Comprehensive evaluation of differential gene expression analysis methods for RNA-seq data. Genome Biol. 14:R95. doi: 10.1186/gb-2013-14-9-r95

PubMed Abstract | CrossRef Full Text | Google Scholar

Raudvere, U., Kolberg, L., Kuzmin, I., Arak, T., Adler, P., Peterson, H., et al. (2019). g:Profiler: a web server for functional enrichment analysis and conversions of gene lists (2019 update). Nucleic Acids Res. 47, W191–W198. doi: 10.1093/nar/gkz369

PubMed Abstract | CrossRef Full Text | Google Scholar

Rinn, J. L., and Chang, H. Y. (2012). Genome regulation by long noncoding RNAs. Annu. Rev. Biochem. 81, 145–166. doi: 10.1146/annurev-biochem-051410-092902

PubMed Abstract | CrossRef Full Text | Google Scholar

Robinson, M. D., McCarthy, D. J., and Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics 26, 139–140. doi: 10.1093/bioinformatics/btp616

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmena, L., Poliseno, L., Tay, Y., Kats, L., and Pandolfi, P. P. (2011). A ceRNA hypothesis: the rosetta stone of a hidden RNA language? Cell 146, 353–358. doi: 10.1016/j.cell.2011.07.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Schennink, A., Bovenhuis, H., Léon-Kloosterziel, K. M., van Arendonk, J. A., and Visker, M. H. (2009). Effect of polymorphisms in the FASN, OLR1, PPARGC1A, PRL and STAT5A genes on bovine milk-fat composition. Anim. Genet. 40, 909–916. doi: 10.1111/j.1365-2052.2009.01940.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Shore, A. N., Kabotyanski, E. B., Roarty, K., Smith, M. A., Zhang, Y., Creighton, C. J., et al. (2012). Pregnancy-induced noncoding RNA (PINC) associates with polycomb repressive complex 2 and regulates mammary epithelial differentiation. PLoS Genet. 8:e1002840. doi: 10.1371/journal.pgen.1002840

PubMed Abstract | CrossRef Full Text | Google Scholar

Silver, N., Best, S., Jiang, J., and Thein, S. L. (2006). Selection of housekeeping genes for gene expression studies in human reticulocytes using real-time PCR. BMC Mol. Biol. 7:33. doi: 10.1186/1471-2199-7-33

PubMed Abstract | CrossRef Full Text | Google Scholar

Sleutels, F., Zwart, R., and Barlow, D. P. (2002). The non-coding Air RNA is required for silencing autosomal imprinted genes. Nature 415, 810–813. doi: 10.1038/415810a

PubMed Abstract | CrossRef Full Text | Google Scholar

Standaert, L., Adriaens, C., Radaelli, E., Van Keymeulen, A., Blanpain, C., Hirose, T., et al. (2014). The long noncoding RNA Neat1 is required for mammary gland development and lactation. RNA 20, 1844–1849. doi: 10.1261/rna.047332.114

PubMed Abstract | CrossRef Full Text | Google Scholar

Sun, L., Luo, H., Bu, D., Zhao, G., Yu, K., Zhang, C., et al. (2013). Utilizing sequence intrinsic composition to classify protein-coding and long non-coding transcripts. Nucleic Acids Res. 41, e166. doi: 10.1093/nar/gkt646

PubMed Abstract | CrossRef Full Text | Google Scholar

Szewczuk, M., Zych, S., and Czerniawska-Pia̧tkowska, E. (2011). Association between IGF1/TasI polymorphism and milk traits of Polish Holstein-Friesian cows. Arch. Anim. Breed. 54, 10–17. doi: 10.1556/ABiol.65.2014.1.3

PubMed Abstract | CrossRef Full Text | Google Scholar

Taddei, I., Deugnier, M. A., Faraldo, M. M., Petit, V., Bouvard, D., Medina, D., et al. (2008). Beta1 integrin deletion from the basal compartment of the mammary epithelium affects stem cells. Nat. Cell Biol. 10, 716–722. doi: 10.1038/ncb1734

PubMed Abstract | CrossRef Full Text | Google Scholar

Tong, C., Chen, Q., Zhao, L., Ma, J., Ibeagha-Awemu, E. M., and Zhao, X. (2017). Identification and characterization of long intergenic noncoding RNAs in bovine mammary glands. BMC Genomics 18:468. doi: 10.1186/s12864-017-3858-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Trapnell, C., Williams, B. A., Pertea, G., Mortazavi, A., Kwan, G., van Baren, M. J., et al. (2010). Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat. Biotechnol. 28, 511–515. doi: 10.1038/nbt.1621

PubMed Abstract | CrossRef Full Text | Google Scholar

Ulitsky, I. (2016). Evolution to the rescue: using comparative genomics to understand long non-coding RNAs. Nat. Rev. Genet. 17, 601–614. doi: 10.1038/nrg.2016.85

PubMed Abstract | CrossRef Full Text | Google Scholar

Vizoso, M., and Esteller, M. (2012). The activatory long non-coding RNA DBE-T reveals the epigenetic etiology of facioscapulohumeral muscular dystrophy. Cell Res. 22, 1413–1415. doi: 10.1038/cr.2012.93

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, C. J. (2006). Involution: apoptosis and tissue remodelling that convert the mammary gland from milk factory to a quiescent organ. Breast Cancer Res. 8:203. doi: 10.1186/bcr1401

PubMed Abstract | CrossRef Full Text | Google Scholar

Watson, C. J., and Khaled, W. T. (2008). Mammary development in the embryo and adult: a journey of morphogenesis and commitment. Development 135, 995–1003. doi: 10.1242/dev.005439

PubMed Abstract | CrossRef Full Text | Google Scholar

Weikard, R., Demasius, W., and Kuehn, C. (2017). Mining long noncoding RNA in livestock. Anim. Genet. 48, 3–18. doi: 10.1111/age.12493

PubMed Abstract | CrossRef Full Text | Google Scholar

Weikard, R., Kühn, C., Goldammer, T., Freyer, G., and Schwerin, M. (2005). The bovine PPARGC1A gene: molecular characterization and association of an SNP with variation of milk fat synthesis. Physiol. Genom. 21, 1–13. doi: 10.1152/physiolgenomics.00103.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, X., Zhou, X., Ding, X., Chu, M., Liang, C., Pei, J., et al. (2019). The selection of reference genes for quantitative real-time PCR in the Ashidan yak mammary gland during lactation and dry period. Animals 9:943. doi: 10.3390/ani9110943

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, C., Mao, X., Huang, J., Ding, Y., Wu, J., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. doi: 10.1093/nar/gkr483

PubMed Abstract | CrossRef Full Text | Google Scholar

Zheng, X., Ning, C., Zhao, P., Feng, W., Jin, Y., Zhou, L., et al. (2018). Integrated analysis of long noncoding RNA and mRNA expression profiles reveals the potential role of long noncoding RNA in different bovine lactation stages. J. Dairy Sci. 101, 11061–11073. doi: 10.3168/jds.2018-14900

PubMed Abstract | CrossRef Full Text | Google Scholar

Zi, X. D., Luo, B., Xia, W., Zheng, Y. C., Xiong, X. R., Li, J., et al. (2018). Characterization of transcriptional complexity during pre-implantation development of the yak (Bos grunniens) using RNA-Seq. Reprod. Domest. Anim. 53, 759–768. doi: 10.1111/rda.13167

PubMed Abstract | CrossRef Full Text | Google Scholar