Transcriptome sequencing to detect the potential role of long non-coding RNAs in bovine mammary gland during the dry and lactation period
BMC Genomics volume 19, Article number: 605 (2018)
It is known that long non-coding RNAs (lncRNAs) play an important role in various biological processes, including cell proliferation, differentiation and apoptosis. However, their functions and profiles in lactation cycle of dairy cows are largely unknown. In this study, lncRNA-seq technique was employed to compare the expression profiles of lncRNAs and mRNAs from Chinese Holstein mammary gland in dry and lactation period.
Totally 3746 differentially expressed lncRNAs (DELs) and 2890 differentially expressed genes (DEGs) were identified from the dry and lactation mammary glands of Holstein cows. Functional enrichment analysis on target genes of lncRNAs indicated that these genes were involved in lactation-related signaling pathways, including cell cycle, JAK-STAT, cell adhesion, and PI3K-Akt signaling pathways. Additionally, the interaction between lncRNAs and their potential miRNAs was predicted and partly verified. The result indicated that the lactation-associated miR-221 might interact with lncRNAs TCONS_00040268, TCONS_00137654, TCONS_00071659 and TCONS_00000352, which revealed that these lncRNAs might be important regulators for lactation cycle.
This study provides a resource for lncRNA research on lactation cycle of bovine mammary gland. Besides, the interaction between lncRNAs and the specific miRNA is revealed. It expands our knowledge about lncRNA and miRNA biology as well as contributes to clarify the regulation of lactation cycle of bovine mammary gland.
Mammary gland is an important organ for the synthesis and secretion of milk, which provides essential nutrients for human and other animal offspring. The development and regression cycles of mammary gland include pregnancy, lactation, and involution, which is regulated by growth factors, hormones, and coding genes [1,2,3]. Janus kinases and signal transducers and activators of transcription ((JAK-STAT) have been shown to function downstream of several peptide hormones and cytokines that are required for postnatal development and secretory function of mammary gland . In mammary epithelial cells, the phosphorylated STAT5A and STAT5B form homodimers and heterodimers to modulate differentiation, survival and proliferation through alterations in cellular gene expression . The phosphatidylinositol 3-kinase-proteinkinase B/mammalian target of the rapamycin (PI3K-Akt/mTOR) signaling pathway regulates a wide range of cellular processes, such as cell proliferation, growth, survival and metastasis , and it is essential for mammary gland development . A conditional knockout of Akt1 averts the extended survival of mammary epithelial cells that express hyperactive STAT5, which indicates that the PI3K-Akt/mTOR pathway is a crucial downstream effector of JAK-STAT signaling .
In the past few years, many studies focused on the functions of protein-coding genes and microRNAs (miRNAs) for the development and regression cycles of mammary gland [9, 10]. ErbB3 played a crucial role in mammary epithelial survival and differentiation during pregnancy and lactation . Forty milk lipid synthesis- and secretion-associated genes, including FADS1, AGPAT6, GPAM, LPIN1, BTN1A1, LPL, CD36, FABP3, ACSL1, ACSS2, ACACA, FASN, SCD, XDH, BDH1, INSIG1, PPARG, and PPARGC1A were verified from dry period to the end of subsequent lactation period . During lactation, MAPK14, FRAP1, EIF4EBP2, GSK3A and TSC1 had an increased expression which revealed that these genes had important roles in milk protein synthesis and secretion . MiRNAs are a kind of non-coding RNA, which can silence or degrade gene expression by targeting the 3’UTR region of coding gene. An increasing number of studies had demonstrated that miRNAs were involved in lactation of mammary gland by regulating their target genes [14,15,16,17,18]. It had been found that miR-27a could regulate cellular triacylglycerol synthesis by targeting PPARG gene in bovine mammary epithelial cells (BMECs) . The overexpression of miR-206 changed the expression of Wnt, Tbx3 and Lef1 genes which were essential for mammary gland development, indicating that miR-206 might be a novel candidate for morphogenesis during the initiation of mammary gland formation . As a downstream regulator of PTEN, miR-486 expressed higher during bovine high quality lactation period and could regulate the secretion of β-casein, lactose and triglyceride in BMECs, which indicated that miR-486 was required for the development of cow mammary gland . In cattle, the expression of miR-221 was found to be higher in peak lactation than in early lactation, suggesting its role in the control of endothelial cell proliferation or angiogenesis . In mice, miR-221 regulated lipid metabolism in mammary epithelial cells and was expressed differentially at various stages during mammary gland development . MiR-212 and miR-132 were necessary for mice mammary gland development and growth by targeting MMP9 gene, especially for mammary epithelial ducts .
Long non-coding RNAs (lncRNAs) are transcripts longer than 200 nucleotides, which have uncovered new layers in the control of various biological processes, including cell proliferation, differentiation and apoptosis . H19 and SRA, two of the earliest identified regulatory lncRNAs, may play a role in the developing mammary gland . Initially, H19 was found to function to restrict embryonic growth, but later evidence showed that H19 had a role in long-term maintenance of adult hematopoietic stem cells . Long noncoding RNA mPINC (mouse pregnancy-induced non-coding RNA) and Zfas1 (Znfx1 antisense 1) had growth-suppressive roles in mammary epithelial cells [24, 25]. Previous observation demonstrated that lncRNA Neat1 could regulate mammary gland morphogenesis and lactation by investigating the proliferation of Neat1-mutant cells in mice . lncRNAs also could act as competing endogenous RNAs (ceRNAs) to control muscle differentiation and involve in goat lactation process [27, 28]. Pregnancy-induced noncoding RNA (PINC) could inhibit terminal differentiation of alveolar cells during pregnancy to prevent abundant milk production and secretion until parturition . The above-mentioned studies showed that lncRNAs had crucial roles in mammary gland cell proliferation and differentiation, which would be an important issue in lactation biology. Mammary gland development and regression was directly related with cow lactation. However, the functions and profiles of lncRNAs in bovine mammary gland during dry and lactation period were largely unknown. The objective of this study was to screen the lncRNAs associated with lactation by lncRNA-seq analysis in bovine mammary gland, which would provide a new perspective for lncRNAs in lactation biology and lay the foundation for their further function study in milk composition synthesis and secretion.
Animals and mammary gland tissue collection
Eight healthy Chinese Holstein dairy cows at the third or fourth parities used in this study were housed in free stall and had access to water and feed ad libitum at the Experimental Farm of Northwest A&F University (Yangling, Shaanxi, China). After intravenous injection of lidocaine hydrochloride, approximately 4 g of mammary gland tissues were harvested via repeated biopsies from four cows at the dry period and four cows at approximately 180 days during lactation period. The tissues were dissected, frozen in RNA later (TaKaRa, Dalian, China) and stored at − 80 °C for further analysis. All experimental and surgery procedures involved in this study were approved by the Experimental Animal Manage Committee of Northwest A&F University (2011–31101684).
Total RNA isolation and quality control for library construction
Total RNA of mammary gland tissues was isolated using Trizol reagent following the manufacturer’s instructions (Invitrogen, CA, USA). The integrity of RNA was detected using RNA Nano 6000 Assay Kit on the Bioanalyzer 2100 system (Agilent Technologies, Santa Clara, USA). RNA purity and concentration was measured using Nanodrop 2000 photometer spectrophotometer (Implen, Los Angeles, USA). The 260/280 ratio for all the samples from mammary gland tissues was about 2.0, and the RNA integrity number (RIN) was ≥8.0. After the determination of RNA purity and quality, 3.0 μg RNA per sample was used and ribosomal RNA was removed using Epicentre Ribo-zero™ rRNA Removal Kit (Epicentre, Madsion, WI, USA) for library construction. The RNA from three individuals in dry and three in lactation stage were pooled, respectively. Subsequently, the libraries were generated from the rRNA-depleted RNA pools using the NEBNext® Ultra™ Directional RNA Library Prep Kit for Illumina® (NEB, USA) following the manufacturer’s recommendations. In order to select cDNA fragments of preferentially 250~ 300 bp in length, the library fragments were purified with AMPure XP system (Beckman Coulter, Beverly, USA).
RNA sequencing, transcriptome assembly, and quantification of gene expression level
The coded samples were clustered using TruSeq PE Cluster reagent (Illumina, CA, USA) according to the manufacturer’s instructions on a cBot Cluster generation system. The libraries were sequenced on an Illumina Hiseq3000 platform and 100-bp paired-end reads were generated after cluster generation. The schematic of lncRNA-seq analysis was shown in Fig. 1. Raw data were first processed using in-house Perl scripts. In this step, clean data were obtained by trimming reads containing adapter, reads containing over 10% of ploy-N, and low-quality reads (> 20% of bases whose Phred scores were < 20) from the raw data. Phred = −10log10 (e), e is defined as the error probability of sequencing for every single base. Q20, Q30 and GC content of the clean data were calculated. Subsequently, Bowtie (v2.0.6)  and Tophat2  (v2.0.9) was used to align paired-end clean reads to the reference genome (version GCA_000003055.5_Bos_taurus_UMD_3.1.1). The default parameters for Tophat2 were set as ‘-read-mismatches=2 (≤2 mismatches are allowed) and -read-gap-length=2’ (≤2 gaps are allowed). The mapped reads were assembled using Cufflinks (v2.1.1) in a reference-based approach . Cufflinks was run with ‘min-frags-per-transfrag = 0’ and ‘-library-type’, other parameters were set as default. Fragments per kb for a million reads (FPKM) of both coding genes and lncRNAs were analyzed using Cuffdiff (v2.1.1) . Differential expression analysis between two groups was performed using the DESeq R package (1.8.3). The P-values were adjusted using the Benjamini and Hochberg method . P-adjust (q-value) < 0.05 and the absolute value of fold change≥2 were set as the threshold for significantly differential expression.
Identification the annotated and novel lncRNAs
NONCODE database was used to characterize the annotated lncRNAs in bovine from the assembled transcripts. To identify bovine novel lncRNAs, the steps were followed as Wang et al. (2017) described  with little modification. (1) transcripts with length < 200 bp were removed; (2) transcripts with predicted ORF > 300 were removed; (3) transcripts were compared with mRNA, rRNA, tRNA, snRNA, snoRNA and pre-miRNA (https://www.ncbi.nlm.nih.gov/) using Cuffcompare v2.1.1 to remove the same or similar transcripts . (4) transcripts with FPKM < 1 were removed; (5) transcripts that did not pass the protein-coding-score test were removed using the Coding Potential Calculator (CPC)  and PFAM database .
Prediction and functional enrichment analysis of lncRNA target genes
To reveal the potential function of lncRNAs, their target genes were predicted in trans and cis. For cis role, it refers to lncRNA’s action on neighboring target genes. In this study, the coding genes from 100 kb upstream and downstream of an lncRNA were searched. The trans role refers to the influence of lncRNAs on other genes at expression level. RNAplex bioinformatics software (http://www.bioinf.uni-leipzig.de/~htafer/index.html) was used to predict lncRNA target genes in trans. Genome distribution of differentially expressed lncRNAs and mRNAs was illustrated with Circos (http://circos.ca/). GO enrichment analysis was performed with DAVID database (https://david.ncifcrf.gov/). As to KEGG analysis, differentially expressed genes were analyzed with KEGG online website (http://www.genome.jp/kegg/). Protein-protein interaction network between differentially expressed genes were analyzed by STRING database (https://string-db.org/), and further visualized with Cytoscape (http://www.cytoscape.org/).
Prediction of potential miRNAs interacted with lncRNAs
To obtain potential miRNAs interacted with lncRNAs, PicTar (https://pictar.mdc-berlin.de/), PITA (https://genie.weizmann.ac.il/pubs/mir07/mir07_prediction.html), and RNAhybrid (https://bibiserv.cebitec.uni-bielefeld.de/rnahybrid/) were used to predict the candidate miRNAs interacted with lncRNAs. The miRNAs shared by the above three tools were selected as the candidate miRNAs to assume the mechanism of lncRNAs interacted with lactation. The potential target genes of miR-221 were predicted by Targetscan (http://www.targetscan.org/vert_71/), PITA and miRanda. The genes shared by the above three tools were selected as the candidate target genes.
Verification of sequencing data using qRT-PCR
To confirm the sequencing results, quantitative real time PCR (qRT-PCR) method was used to measure the relative expression of DEGs and DELs. The total RNAs from the cows used for lncRNA-seq were also used for qRT-PCR. The first-strand cDNA was obtained using a PrimeScript RT reagent Kit (TaKaRa, Dalian, China) following the manufacturer’s instructions. The qRT-PCR was performed in triplicate using SYBR® Premix Ex Taq™ II (TaKaRa) on the Bio-Rad CFX96 Touch™ Real Time PCR Detection System (Bio-Rad, Hercules, CA, USA). The qRT-PCR protocol was initiated with an incubation of 3 min at 95 °C, followed by 40 cycles of 95 °C for 12 s and optimized annealing temperature for 40 s. The relative expression of DEGs and DELs were analyzed using the 2−ΔΔCt method and normalized by GAPDH gene . The qRT-PCR primers for DELs and DEGs were designed with Primer Premier 6.0 (Premier, British Columbia, Canada) spanning at least one intron and shown in Additional file 1.
BMECs transfection and expressions of four lncRNAs and ErBb3
BMECs were cultured in Dubecco’s modified Eagle medium (DMEM) /F12 containing 10% fetal bovine serum and 1% penicillin streptomycin (All from Gibco, Grand Island, NY, USA) at 37 °C with 5% CO2. The cells were seeded in 24-well plates, then transfected with miR-221 mimic, mimic-NC, miR-221 inhibitor and inhibitor-NC (GenePharma, Shanghai, China) at approximately 50% confluence using Lipofectamine 2000 Transfection Reagent (Invitrogen, USA) according to the manufacturer’s instructions, respectively. Briefly, Lipofectamine 2000 Transfection Reagent (1 μL) was diluted in 25 μL of Opti-MEM (Invitrogen) with miR-221 mimic, mimic-NC, inhibitor, or inhibitor-NC to yield final concentrations of 20, 20, 40, and 40 nM, respectively. The mixture was incubated at room temperature for 20 min and added to the BMECs. The transfection efficiency was assessed by fluorescence microscopy after 48 h.
Total RNA was extracted from the transfected BMECs using Trizol reagent (Invitrogen) 72 h after transfection. Then, the cDNA was obtained as above mentioned method. qRT-PCR protocol was also performed as above mentioned. Primers for the four lncRNAs and ErBb3 are shown in Additional file 1.
Overview of sequencing data in cow mammary gland
A total of 58,411,766 and 69,114,038 raw reads were produced from cow mammary glands using the Illumina Hiseq3000 platform in dry and lactation periods, respectively. After discarding adaptor sequences and low-quality sequences, 54,805,136 and 65,069,892 corresponding clean reads were obtained, and the percentage of clean reads was 93.83 and 94.15%, respectively (Table 1). The whole expression feature of transcripts was shown in Fig. 2a. The expression level of the transcripts in lactation was slightly higher than that in dry period. Similarly, 64,316 and 61,791 known transcripts, 44,651 and 43,094 known mRNAs were also obtained in dry and lactation period of cow mammary gland, respectively. Of those transcripts, 928 were novel in dry stage and 841 were in lactation, respectively.
The genome distribution of clean reads was shown in Fig. 2b. The results revealed that most clean reads in the two periods were located in exonic region, only few reads were in intergenic region. And about a quarter of the total clean reads were in the intronic region of the bovine genome.
GO and KEGG enrichment analysis of DEGs
Totally 2890 differentially expressed genes (DEGs) were obtained according the fold change≥2 and q < 0.01 (Additional file 2). And 2300 genes were down-regulated, whereas 590 were up-regulated in lactation compared with dry period. Then functional enrichment analysis of the 590 significantly up-regulated genes was performed based on their relative expression and fold changes. As illustrated in Fig. 3a, 69 significantly enriched GO terms (P < 0.05) were identified, such as negative regulation of cell growth (GO: 0030308), growth factor activity (GO: 0008083), positive regulations of ERK1 and ERK2 cascade (GO: 0070374), and sodium channel complex (GO: 0034706), etc. ERK1 and ERK2 had been suggested to play an important role in regulating cell invasion, cell proliferation, and colony formation in triple-negative breast cancer cell lines . Growth factors, especially the epidermal growth factors (EGFs) and insulin-like growth factors (IGFs), were involved in development of normal mammary gland and pathogenesis of breast cancer .
Two thousand two hundred eight significantly down-regulated genes were also selected to carry out the functional enrichment analysis. Forty five terms, including cell division (GO: 0051301), cell adhesion (GO: 0007155), cell communication (GO: 0007154), cation transmembrane transport (GO: 0098655), and cell surface receptor signaling pathway (GO: 0035556), were shown in Fig. 3b. It was worthy to notice that the process of cell adhesion and division were closely associated with mammary gland architecture construction, maintenance, development, and lactation [41, 42].
Meanwhile, KEGG results from the significantly up- and down-regulated genes indicated that 53 significantly signaling pathways were enriched (Fig. 3c, P < 0.05), such as cell adhesion molecules (CAMs), PI3K-Akt, PPAR, TNF, cytokine-cytokine, cell cycle, and Wnt signaling pathways. Previous studies had determined that Wnt, PPAR, CAMs, PI3K-Akt, and TNF signaling pathways could regulate mammary gland development, milk fat synthesis, mastitis, BMECs proliferation and apoptosis [43,44,45,46,47]. The expression pattern of genes involved in PI3K-Akt and PPAR signaling was displayed in Fig. 3d. The heat map showed that PIK3R3, CSF3, TNC, TLR2, GNG11, DOIT43, NOS3, THBS3, IL3RA, FN1, SORBS1, FAPPS, SLC27A, ACP7, GK, ACSL4 and ANGPTL4 genes expressed higher in lactation than that in dry period at mRNA level.
DELs and their potential interacted miRNAs involved in lactation
Totally 23,495 expressed lncRNAs were found in the two different periods, of which 5893 were novel (Additional files 3 and 4) and 17,602 were annotated in NONCODE (Additional file 5). A total of 3746 significantly differentially expressed lncRNA transcripts were found in lactation compared with that in dry period, including 2732 down- and 1014 up-regulated lncRNA transcripts (fold change≥2, P < 0.05) (Additional files 6 and 7). The genome wide distribution of DELs and DEGs was shown in Fig. 4a. It could be seen that the DELs and DEGs was harmoniously located on autosomes and X chromosome. The stage-specific expressed lncRNAs and genes were shown in Additional files 8 and 9. One hundred forty five lncRNAs were identified and subjected to further analysis according their expression fold changes and genomic locations. The potential interacted miRNAs with the selected 145 lncRNAs were predicted by miRanda, PITA and RNAhybrid, and partial result was shown in Table 2. It could be seen that miR-103, miR-21-3p, miR-27a-5p, miR-107, and miR-24-3p were predicted to interact with lncRNA TCONS_00120917. The interaction network between candidate lncRNAs and their potential miRNAs was constructed (Fig. 4b). It showed that several miRNAs might interact with multiple lncRNAs, such as miR-221 interacting with lncRNAs TCONS_00032404, TCONS_00032444, TCONS_00040268, TCONS_00137654, TCONS_00071659, TCONS_00143274 and TCONS_00000352. Coincidentally, it had demonstrated that the expression of miR-221 was found to be higher in peak lactation than in early lactation, suggesting a role in the control of endothelial cell proliferation or angiogenesis which was closely related with lactation , so those above mentioned lncRNAs could be considered as important candidates for lactation.
The interaction between lncRNAs and miR-221 as well as miR-221 and ErBb3
To verify the interaction between lncRNAs and miR-221 and the interaction between miR-221 and ErBb3, overexpression and inhibition of miR-221 was employed in BMECs. Our data showed that the expressions of the four lncRNAs (TCONS_00040268, TCONS_00137654, TCONS_00071659, and TCONS_00000352) decreased with the overexpression of miR-221. On the contrary, their expressions were increased with the inhibition of miR-221 in BMECs (Fig. 5a), which indicated that miR-221 could interact with these four lncRNAs; Furthermore, the overexpression of miR-221 could reduce the expression of ErBb3 gene in miR-221 mimic group compared with mimic-NC group. While the inhibition of miR-221 increased ErBb3 gene expression, which indicated that ErBb3 gene was the target of miR-221 (Fig. 5b).
Interaction network between the candidate lncRNAs and their potential target genes
To reveal the potential function of these selected lncRNAs, their target genes were predicted in cis and trans roles, and we found that the mammary gland biology related genes, such as PRLR, FAS, and MAP3K7 genes, could be targeted by several lncRNAs in cis or trans (Table 3). In addition, the network between lncRNAs and target genes was constructed. The potential target genes of those lncRNAs included S100A1, UPK1B, BARD1, FGF2, IGFBP1, ALOX15, KLH34, SYT12, WNT10A, SULF1, SLC10A6 and CLE7A (Fig. 6). Moreover, the expression level of those genes in dry period were significantly higher than that in lactation period in our study while the corresponding lncRNAs were significantly lower, suggesting an opposite expression trend exists between those genes and corresponding lncRNAs. Previous data had demonstrated that IL1B might function as a retinoic acid induced gene and could inhibit the proliferation of normal human mammary epithelial cells . Therefore, it could be deduced that TCONS_00075230 might involve in regulating the proliferation of mammary epithelial cells by interacting with IL1B gene in cow mammary gland.
GO and KEGG analysis of lncRNAs potential target genes
A total of 47 potential target genes in cis of lncRNAs (Additional file 10) were selected to carry out functional enrichment analysis. Our results showed that nine significantly enriched GO terms focused on the negative regulation of JAK-STAT cascade (GO: 0046426), positive regulation of transcription from RNA polymerase II promoter (GO: 0045944), cytokine-mediated signaling pathway (GO: 0019221), and positive regulation of phosphorylation (GO: 0042327), etc. (Fig. 7a, P < 0.05). Previous study had suggested that JAK-STAT pathway made a marked contribution to dairy milk production traits . In addition, phosphorylation regulation of mammary S6 K1 and eIF2 genes could control milk protein yield .
Similarly, a total of 459 potential target genes in trans of lncRNAs (Additional file 11) were selected to perform functional enrichment analysis (Fig. 7a, P < 0.05). These terms were protein glycosylation (GO: 0006486), positive regulation of cell division (GO: 0051781), DNA-directed RNA polymerase III complex (GO: 0005666), and growth factor binding signaling pathway (GO: 0019838), etc. It was worth mentioning that cell division signaling pathway had been demonstrated to be closely associated with mammary gland architecture construction and maintenance . Additionally, protein glycosylation in mammal membrane played important roles in milk quality and biomodulator properties .
Meanwhile, the KEGG results from the target genes of candidate lncRNAs in trans and cis illustrated 15 pathways (Fig. 6b). Among these pathways, PI3K-Akt and MAPK pathways had been confirmed to be associated with cow lactation [43, 52].
Interaction network of protein-protein
Combine the bioinformatics analysis of DEGs and target genes of DELs, the interaction network between proteins was produced using String software (Fig. 8). The results revealed that protein-protein interaction focused on CCN (CCND1, CCNA1, CCNB2, CCNA2, CCNB1, CCNE1, and CCNE2) and CDC protein family (CDC6, CDC20, CDC25B, CDC25C, and CDCA8). CCN and CDC protein families had been demonstrated to be involved in cyclin and cell division cycle, respectively [53, 54].
Seven lncRNAs and five genes, which were significantly differentially expressed between the two periods, were selected to verify the transcriptome sequencing data using qRT-PCR. The seven lncRNAs were selected based on maximal fold change (and P < 0.05) and location of the genome. For the five genes, they were selected because of their fold change (and P < 0.05) and potential function in lactation. The results showed that the relative expression of the selected lncRNAs and genes was consistent with their sequencing data (Fig. 9).
Increasing data had shown that lncRNAs could regulate gene expression both at transcriptional and post-transcriptional levels, including the regulation of splicing, mRNA processing, and translation . It had been demonstrated that some lncRNAs implicated in breast normal development and cancer based on their expression pattern, function and localization in human and mouse [56, 57]. However, limited researches had been reported the regulation mechanism of lncRNAs on bovine lactation. Koufariotis et al. (2015) identified and annotated novel lncRNA transcripts in the bovine genome across 18 tissues (including mammary gland) via RNA sequencing. In addition, they found most lncRNA in bovine mammary gland were downregulated compared with other tissues . Tong et al. (2017) identified 36 lincRNAs located in milk related quantitative trait loci (QTL) from five RNA-seq datasets of bovine mammary glands whereas one lincRNA was within clinical mastitis QTL region, which indicated that these lncRNAs were involved in many biological functions including susceptibility to clinical mastitis as well as milk quality and production . In our study, lactation-related lncRNAs were assessed by directly lncRNA-seq from cow mammary glands at the dry and lactation period, respectively. Given the functional enrichment analysis of lncRNAs target genes, 24 lncRNAs were identified to have potential role in lactation (Table 3). Previous study had suggested that lncRNAs could regulate biological processes by sponging miRNAs , so 15 lncRNAs interacted with lactation-associated miRNAs were also considered as potential regulators for lactation (Table 2). Among those lncRNAs, five lncRNAs were shared by the above two prediction methods.
LncRNAs might regulate expression of lactation-associated genes in trans or cis
It was well recognized that a series of genes involved in lactation initiation, maintenance as well as mammary gland growth, development and breast cancer by direct or indirect regulation [8, 9, 11, 60]. In this study, the DEGs functional enrichment results showed that these genes were related to some biological processes, including cell division, adhesion, cycle, ERK1 and ERK2 cascade, PI3K-Akt, PPAR, and TNF pathway, which were closely associated with lactation [43,44,45,46].
LncRNAs can regulate neighboring gene expression, therefore high expression correlations exert between lncRNAs and mRNAs (known as in cis) . LncRNAs also can change the expression of distant mRNAs through the pairing of lncRNAs-mRNA (known as in trans). Bioinformatics analysis of lncRNAs target genes in trans showed that these genes played an important role in some pathways, such as MAPK, PI3K-Akt, prolactin, NF-kappa B, and Toll-like receptor signaling pathways, which played important roles during mammary gland development and lactation [4, 8, 62]. Consequently, many lncRNAs might function through targeting mRNA which played important roles in mammary gland from non-lactation to lactation period. For example, lncRNA TCONS_00075230 and IL1B had an opposite expression trend between dry and lactation period, and IL1B had reported to inhibit the proliferation of normal human mammary epithelial cells function as a retinoic acid induced gene . So TCONS_00075230 might involve in lactation process through altering the expression of IL1B.
Regulating role of lncRNA-miRNA-mRNA network
It was known that the regulatory networks of lncRNAs, miRNA, and ceRNAs communicated with each other to regulate gene expression . LncRNA SNHG7 could accelerate prostate cancer proliferation and cell cycle progression through cyclin D1 by sponging miR-503 . The axis of lncRNA H19-miR-675-TGFBI had possible diagnostic and therapeutic potential for advanced prostate cancer . LncRNA APF could control the expression of miR-188-3p, which suppressed myocardial infarction and autophagy by targeting ATG7 . A series of miRNAs had been demonstrated to regulate lactation-associated processes, including lactation cycle, milk fat accumulation, hormone receptor activity, mammary gland involution and development, and lactation activity of mammary epithelial cells [14, 66, 67]. In this study, some candidate lncRNAs, including TCONS_00040268, TCONS_00137654, TCONS_00071659, and TCONS_00000352, could interact with lactation-related miR-221, and miR-221 could target ErBb3 gene. It had known that ErBb3 gene could drive mammary epithelial survival and differentiation during pregnancy and lactation . Therefore, these lncRNAs might play an important role in regulating lactation, which would provide a new insight into lactation process of cattle.
In this study, 3746 significantly dysregulated lncRNA transcripts were found in lactation period compared with that in dry period, including 2732 down- and 1014 up-regulated lncRNA transcripts (fold change≥2, P < 0.05). Functional enrichment analysis of target genes and interacted miRNAs prediction of 34 lncRNAs indicated these lncRNAs might be important regulators for lactation in dairy cattle. This study would provide a resource for bovine lncRNA study involving in lactation biology.
Bovine mammary epithelial cells
Competing endogenous RNAs
Differentially expressed genes
Differentially expressed lncRNAs
Dubecco’s modified Eagle medium
Long non-coding RNAs
Quantitative real time PCR
Watson CJ, Khaled WT. Mammary development in the embryo and adult: a journey of morphogenesis and commitment. Development. 2008;135(6):995–1003.
Topper YJ, Freeman CS. Multiple hormone interactions in the developmental biology of the mammary gland. Physiol Rev. 1980;60(4):1049–106.
Liu X, Robinson GW, Wagner KU, Garrett L, Wynshawboris A, Hennighausen L. Stat5a is mandatory for adult mammary gland development and lactogenesis. Genes Dev. 1997;11(2):179–86.
Rädler PD, Wehde BL, Wagner KU. Crosstalk between STAT5 activation and PI3K/AKT functions in normal and transformed mammary epithelial cells. Mol Cell Endocrinol. 2017;451(8):31–9.
Furth PA, Nakles RE, Millman S, Diaz-Cruz ES, Cabrera MC. Signal transducer and activator of transcription 5 as a key signalling pathway in normal mammary gland developmental biology and breast cancer. Breast Cancer Res. 2011;13(5):220.
Ersahin T, Tuncbag N, Cetin-Atalay R. The PI3K/AKT/mTOR interactive pathway. Mol BioSyst. 2015;11(7):1946.
Sobolewska A, Gajewska M, Zarzyńska J, Gajkowska B, Motyl T. IGF-I, EGF, and sex steroids regulate autophagy in bovine mammary epithelial cells via the mTOR pathway. Eur J Cell Biol. 2009;88(2):117–30.
Schmidt JW, Wehde BL, Sakamoto K, Triplett AA, Anderson SM, Tsichilis PN, Leone G, Wagner KU. Stat5 regulates the phosphatidylinositol 3-kinase/Akt1 pathway during mammary gland development and tumorigenesis. Mol Cell Biol. 2004;34(7):1363–77.
Dai WT, Zou YX, White RR, Liu JX, Liu HY. Transcriptomic profiles of the bovine mammary gland during lactation and the dry period. Funct Integr Genomics. 2018;18(2):125–40.
Li Z, Liu HY, Jin XL, Lo LJ, Liu JX. Expression profiles of microRNAs from lactating and non-lactating bovine mammary glands and identification of miRNA related to lactation. BMC Genomics. 2012;13(1):731.
Williams MM, Vaught DB, Joly MM, Hicks DJ, Sanchez V, Owens P, Rahman B, Elion DL, Balko JM, Cook RS. ErbB3 drives mammary epithelial survival and differentiation during pregnancy and lactation. Breast Cancer Res. 2017;19(1):105.
Bionaz M, Loor JJ. Gene networks driving bovine milk fat synthesis during the lactation cycle. BMC Genomics. 2008;9(1):366.
Bionaz M, Loor JJ. Gene networks driving bovine milk protein synthesis during the lactation cycle. Bioinform Biol Insights. 2011;5:83–98.
Do DN, Li R, Dudemaine PL, Ibeagha-Awemu EM. MicroRNA roles in signaling during lactation: an insight form differential expression, time course and pathway analyses of deep sequence data. Sci Rep. 2017;7:44605.
Tang KQ, Wang YN, Zan LS, Yang WC. miR-27a controls triacylglycerol synthesis in bovine mammary epithelial cells by targeting peroxisome proliferator-activated receptor gamma. J Dairy Sci. 2017;100(5):4102–12.
Lee MJ, Yoon KS, Cho KW, Kim KS, Jung HS. Expression of miR-206 during the initiation of mammary gland development. Cell Tissue Res. 2013;353(3):425–33.
Li D, Xie XJ, Wang J, Bian YJ, Li QZ, Gao XJ, Wang CM. miR-486 regulates lactation and targets the PTEN gene in cow mammary glands. PLoS One. 2015;10(3):e0118284.
Wang M, Moisá S, Khan MJ, Wang J, Bu D, Loor JJ. MicroRNA expression patterns in the bovine mammary gland are affected by stage of lactation. J Dairy Sci. 2012;95(11):6529–35.
Chu MQ, Zhao Y, Yu S, Hao YN, Zhang PF, Feng YN, Zhang HF, Ma DX, Liu J, Cheng M, et al. MicroRNA-221 may be imvolved in lipid metabolism in mammary epithelial cell. Intl J Biochem Cell Biol. 2018;97:118–27.
Ucar A, Vafaizadeh V, Jarry H, Fiedler J, Klemmt PAB, Thum T, Groner B, Chowdhury K. miR-212 and miR-132 are required for epithelial stromal interactions necessary for mouse mammary gland development. Nat Genet. 2010;42(12):1101–8.
Guttman M, Rinn JL. Modular regulatory principles of large non-coding RNAs. Nature. 2012;482(7385):339–46.
Shore AN, Rosen JM. Regulation of mammary epithelial cell homeostasis by lncRNAs. Int J Biochem Cell Biol. 2014;54:318–30.
Venkatraman A, He XC, Thorvaldsen JL, Sugimura R, Perry JM, Tao F, Zhao M, Christenson MK, Sanchez R, Yu JY, et al. Maternal imprinting at the H19-Igf2 locus maintains adult haematopoietic stem cell quiescence. Nature. 2013;500(7462):345–9.
Ginger MR, Shore AN, Contreras A, Rijnkels M, Miller J, Gonzalez-Rimbau MF, Rosen JM. A noncoding RNA is a potential marker of cell fate during mammary gland development. Proc Natl Acad Sci U S A. 2006;103(15):5781–6.
Askarian-Amiri ME, Crawford J, French JD, Smart CE, Smith MA, Clark MB, Ru K, Mercer TR, Thompson ER, Lakhani SR, et al. SNORD-host RNA Zfas1 is a regulator of mammary development and a potential marker for breast cancer. RNA. 2011;17(5):878–91.
Standaert L, Adriaens C, Radaelli E, Van Keymeulen A, Blanpain C, Hirose T, Nakagawa S, Marine JC. The long noncoding RNA Neat1 is required for mammary gland development and lactation. RNA. 2014;20(12):1844–9.
Cesana M, Cacchiarelli D, Legnini I, Santini T, Sthandier O, Chinappi M, Tramontano A, Bozzoni I. A long noncoding RNA controls muscle differentiation by functioning as a competing endogenous RNA. Cell. 2011;147:358–69.
Yu S, Zhao Y, Lai F, Chu M, Hao Y, Feng Y, Zhang H, Liu J, Cheng M, Li L, Shen W, Min L. LncRNA as ceRNAs may be involved in lactation process. Oncotarget. 2017;8(58):98014–28.
Shore AN, Kabotyanski EB, Roarty K, Smith MA, Zhang Y, Creighton CJ, Dinger ME, Rosen JM. Pregnancy-induced noncoding RNA (PINC) associates with polycomb repressive complex 2 and regulates mammary epithelial differentiation. PLoS Genet. 2012;8(7):e1002840.
Langmead RBB, Trapnell C, Pop M, Salzberg LS. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol. 2009;10(3):R25.
Trapnell C, Pachter L, Salzberg SL. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics. 2009;25(9):1105–11.
Trapnell C, Roberts A, Goff L, Pertea G, Kim D, Kelley DR, Pimentel H, Salzberg SL, Rinn JL, Pachter L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat Protoc. 2012;7(3):562–78.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, van Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation. Nat Biotechnol. 2010;28:511–5.
Storey JD. The positive false discovery rate: a Bayesian interpretation and the q-value. Ann Stat. 2003;31(6):2013–35.
Wang SH, Ge W, Luo ZX, Guo Y, Jiao BL, Qu L, Zhang ZY, Wang X. Integrated analysis of coding genes and non-coding RNAs during hair follicle cycle of cashmere goat (Capra hircus). BMC Genomics. 2017;18(1):767.
Kong L, Zhang Y, Ye ZQ, Liu XQ, Zhao SQ, Wei L, Gao G. CPC: assess the protein-coding potential of transcripts using sequence features and support vector machine. Nucleic Acids Res. 2007;35:W345–9.
Punta M, Coggill PC, Eberhardt RY, Mistry J, Tate J, Boursnell C, Pang N, Forslund K, Ceric G, Clements J, et al. The Pfam protein families database. Nucleic Acids Res. 2012;40(D1):D290–301.
Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(T)(−Delta Delta C) method. Methods. 2001;25(4):402–8.
Park S, Sung Y, Jeong J, Choi M, Lee J, Kwon W, Jang S, Park SJ, Kim HS, Lee MH, et al. hMAGEA2 promotes progression of breast cancer by regulating Akt and Erk1/2 pathways. Oncotarget. 2017;8(23):37115–27.
Hynes NE, Watson CJ. Mammary gland growth factors: roles in normal development and in cancer. Cold Spring Harb Perspect Biol. 2010;2(8):a003186.
Streuli CH. Cell adhesion in mammary gland biology and neoplasia - preface. J Mammary Gland Biol. 2003;8(4):375–81.
Knight CH. The importance of cell division in udder development and lactation. Livest Prod Sci. 2000;66(2):169–76.
Samant GV, Sylvester PW. Gamma-Tocotrienol inhibits ErbB3-dependent PI3K/Akt mitogenic signalling in neoplastic mammary epithelial cells. Cell Proli. 2006;39(6):563–30.
Yee LD, Guo Y, Bradbury J, Suster S, Clinton SK, Seewaldt VL. The antiproliferative effects of PPAR gamma ligands in normal human mammary epithelial cells. Breast Cancer Res Treat. 2003;78:179–92.
Frey RS, Singletary KW. Genistein activates p38 mitogen-activated protein kinase, inactivates ERK1/ERK2 and decreases Cdc25C expression in immortalized human mammary epithelial cells. J Nutr. 2003;133(1):226–31.
Hojilla CV, Jackson HW, Khokha R. TIMP3 regulates mammary epithelial apoptosis with immune cell recruitment through differential TNF dependence. PLoS One. 2011;6(10):e26718.
Kessenbrock K, Smith P, Steenbeek SC, Pervolarakis N, Kumar R, Minami Y, Goga A, Hinck L, Werb Z. Diverse regulation of mammary epithelial growth and branching morphogenesis through noncanonical Wnt signaling. Proc Natl Acad Sci U S A. 2017;114(12):3121–6.
Liu LM, Gudas LJ. Retinoic acid induces expression of the interleukin-1b gene in cultured normal human mammary epithelial cells and in human breast carcinoma lines. J Cell Physiol. 2002;193:244–52.
Arun SJ, Thomson PC, Sheehy PA, Khatkar MS, Raadsma HW, Williamson P. Targeted analysis reveals an important role of JAK-STAT-SOCS genes for milk production traits in Australian dairy cattle. Front Genet. 2015;6(187):342.
Toerien CA, Trout DR, Cant JP. Nutritional stimulation of milk protein yield of cows is associated with changes in phosphorylation of mammary eukaryotic initiation factor 2 and ribosomal s6 kinase 1. J Nutr. 2010;140(2):285–92.
O'Riordan N, Kane M, Joshi L, Hickey RM. Structural and functional characteristics of bovine milk protein glycosylation. Glycobiology. 2014;24(3):220–36.
Jurek B, Slattery DA, Maloumby R, Hillerer K, Koszinowski S, Neumann ID, van den Burg EH. Differential contribution of hypothalamic MAPK activity to anxiety-like behaviour in virgin and lactating rats. PLoS One. 2012;7(5):e37060.
Yang Y, Gu CY, Luo C, Li F, Wang M. BUB1B promotes multiple myeloma cell proliferation through CDC20/CCNB axis. Med Oncol. 2015;32(3):81.
Cooper CR, Szaniszlo PJ. Evidence for 2 cell-division cycle (CDC) genes that govern yeast bud emergence in the pathogenic fungus wangiella-dermatitidis. Infect Immun. 1993;61(5):2069–81.
Hansji H, Leung EY, Baguley BC, Finlay GJ, Askarian-Amiri ME. Keeping abreast with long non-coding RNAs in mammary gland development and breast cancer. Front Genet. 2014;5:379.
Huang J, Zhou N, Watabe K, Lu Z, Wu F, Xu M, Mo YY. Long non-coding RNA UCA1 promotes breast tumor growth by suppression of p27 (Kip1). Cell Death Dis. 2014;5:e1008.
Zhang YF, Wagner EK, Guo XY, May I, Cai QY, Zheng W, He CY, Long JR. Long intergenic non-coding RNA expression signature in human breast cancer. Sci Rep. 2016;6:37821.
Koufariotis LT, Chen YPP, Chamberlain A, Jagt CV, Hayes BJ. A catalogue of novel bovine long noncoding RNA across 18 tissues. PLoS One. 2015;10(10):e0141225.
Tong C, Chen QL, Zhao LL, Ma JF, Ibeagha-Awemu EM, Zhao X. Identification and characterization of long intergenic noncoding RNAs in bovine mammary glands. BMC Genomics. 2017;18:468.
Bach K, Pensa S, Grzelak M, Hadfield J, Adams DJ, Marioni JC, Khaled WT. Differentiation dynamic of mammary epithelial cells revealed by single-cell RNA sequencing. Nat Commun. 2017;8(1):2128.
Guttman M, Amit I, Garber M, French C, Lin MF, Feldser D, Huarte M, Zuk O, Carey BW, Cassady JP, et al. Chromatin signature reveals over a thousand highly conserved large non-coding RNAs in mammals. Nature. 2009;458(7235):223–7.
Garcia MC, Lopez M, Gualillo O, Seoane LM, Dieguez C, Senaris RM. Hypothalamic levels of NPY, MCH, and prepro-orexin mRNA during pregnancy and lactation in the rat: role of prolactin. FASEB J. 2003;17(11):1392–400.
Qi H, Wen B, Wu Q, Cheng W, Lou J, Wei J, Huang J, Yao X, Weng G. Long noncoding RNA SNHG7 accelerates prostate cancer proliferation and cycle progression through cyclin D1 by sponging miR-503. Biomed Phamacother. 2018;102:326–32.
Zhu MJ, Chen Q, Liu X, Sun Q, Zhao X, Deng R, Wang YL, Huang J, Xu M, Yan JS, et al. lncRNA H19/miR-675 axis represses prostate cancer metastasis by targeting TGFBI. FEBS J. 2014;281(16):3766–75.
Wang K, Liu CY, Zhou LY, Wang JX, Wang M, Zhao B, Zhao WK, Xu SJ, Fan LH, Zhang XJ, et al. APF lncRNA regulates autophagy and myocardial infarction by targeting miR-188-3p. Nat Commun. 2015;6:6779.
Lin XZ, Luo J, Zhang LP, Wang W, Gou DM. MiR-103 controls milk fat accumulation in goat (Capra hircus) mammary gland during lactation. PLoS One. 2013;8(11):e79258.
Feuermann Y, Kang K, Shamay A, Robinson GW, Hennighausen L. MiR-21 is under control of STAT5 but is dispensable for mammary development and lactation. PLoS One. 2014;9(1):e85123.
This study was supported by the Program for New Excellent Talents in University of China (NCET-13-0485), National natural Science Foundation of China (31772573) and Agricultural Scientific and Technological Innovation Project of Shandong Academy of Agricultural Science (CXGC2016B03).
Availability of data and materials
The RNA-seq datasets supporting the conclusions of this article are available in the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/Traces/sra_sub/). BioProject asscession: PRJNA482783. BioSamples: SAMN09714469 and SAMN09714468.
The datasets supporting the conclusions of this article are included within the article and its additional files.
All the experimental procedures and animal care with cows involved in this study were approved by the Experimental Animal Manage Committee of Northwest A&F University (China).
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
qRT-PCR primers for lncRNAs and coding genes. (DOCX 21 kb)
The expression profile of DEGs in the dry and lactation periods. (XLSX 387 kb)
The expression profile of novel lncRNAs in the dry and lactation periods. (XLSX 980 kb)
The sequence information of novel lncRNAs. (XLSX 2011 kb)
The expression profile of annotated lncRNAs in the dry and lactation periods. (XLSX 1732 kb)
The expression file of novel DELs in the dry and lactation periods. (XLSX 230 kb)
The expression file of annotated DELs in the dry and lactation periods. (XLSX 224 kb)
Stage specific expressed lncRNAs in the dry or lactation periods. (XLSX 41 kb)
Stage specific expressed genes in the dry or lactation periods. (XLSX 21 kb)
Predicted target genes and mRNAs of DELs (in cis). (XLSX 24 kb)
Predicted target genes and mRNAs of DELs (in trans). (XLSX 84 kb)
About this article
Cite this article
Yang, B., Jiao, B., Ge, W. et al. Transcriptome sequencing to detect the potential role of long non-coding RNAs in bovine mammary gland during the dry and lactation period. BMC Genomics 19, 605 (2018). https://doi.org/10.1186/s12864-018-4974-5
- Lactation cycle
- Cow mammary gland