Identification of long non-coding RNAs in the early growth stage of Holstein mammary gland

Abstract

In Holsteins, emerging evidence suggests that long non-coding RNAs (lncRNAs) are actively involved in the regulation of mammary gland development. However, the lncRNA expression pattern in the early growth stage of Holstein mammary gland remains little known. In this study, we comprehensively analyzed 11 mammary gland tissue RNA-Seq datasets derived from Holstein heifer calves (ca. 8 weeks old) retrieved from the Sequence Read Archive database. We identified a total of 588 lncRNAs in the intergenic region. By comparing the sequences with NONCODEv5 database, 241 novel lncRNAs were detected in the mammary gland. Prediction of cis-regulated target gene revealed that 364 out of 588 lncRNAs potentially targeted 417 protein-coding genes. Furthermore, Gene Ontology enrichment determined that potential target genes were mainly involved in the terms associated with energy metabolism and cell differentiation, such as cellular respiration, respiratory electron transport chain, electron transport chain, and regulation of cell differentiation, which were closely related to mammary gland development. This study provides potential candidate lncRNAs in regulating mammary gland development in the early growth stage of Holsteins, as well as a better general understanding of Holstein mammary gland development.

Keywords:

• Bioinformatics
• Holsteins
• LncRNAs
• mammary gland

Introduction

Holsteins are one of the most popular cow breeds in the world mainly because of their high milk yield. Their mammary glands are an important organ for Holsteins as they affect both the yield and quality of the produced milk intended for human consumption. Therefore, improving the development of the mammary gland has long become an interest in the dairy industry. Mammary gland development is a continuous process in the entire life of Holsteins starting from the embryonic, pre-weaning, pubertal, and until adult development (Macias & Hinck 2012; Vailati-Riboni et al. 2018). In particular, a well development of the mammary gland in pre-weaning may contribute to greater future production (Soberon et al. 2012). The association mechanisms and underlying genetic factor for the regulation of mammary gland development in the early growth stage have recently garnered increased attention.

The vast majority of genome transcripts have been found to be non-coding RNAs (ncRNAs) in eukaryotes (Stein 2004). NcRNAs are functional RNAs that do not encode proteins but play important roles in various cellular processes (Li et al. 2014; Chen & Yang 2015; He et al. 2018). Specifically, long non-coding RNAs (lncRNAs) are a class of ncRNAs that have more than 200 nucleotides and lack functional open reading frames (ORFs) (Derrien et al. 2012; Rinn & Chang 2012; Batista & Chang 2013). They can be classified into four types based on their location in the genome, namely antisense lncRNAs, intronic lncRNAs, sense lncRNAs, and intergenic lncRNAs (lincRNAs).

Interestingly, large-scale genome sequencing has revealed that lncRNAs are positively involved in mammary gland metabolism and development. For example, 28,562 lncRNAs have been found to be expressed in the bovine mammary gland using RNA-Seq, of which 174 were differentially expressed during heat stress, suggesting that lncRNAs were involved in the physiological adaptation of mammary glands (Li et al. 2018). Ibeagha-Awem et al. also found 4955 lncRNAs expressed in the mammary gland of treated cows with a dietary response to 5% linseed oil or 5% safflower oil, demonstrating that lncRNAs expression pattern in the mammary gland could be modulated by nutritional signal (Ibeagha-Awemu et al. 2018). Yang et al., on the other hand, identified 3746 differentially expressed lncRNAs from the dry and lactation mammary glands of Holstein cows, suggesting that the lncRNAs were important regulators for lactation cycle (Yang et al. 2018). Meanwhile, Tong et al. reported 36 lincRNA in the Holstein mammary gland, which were related to 172 milk quantitative trait loci (QTLs) with one lincRNA found within clinical mastitis QTL region (Tong et al. 2017).

However, the lncRNA expression patterns in the early growth stage of Holstein mammary gland remain largely unknown. To investigate the lncRNA profiles and identify lncRNAs that were involved in mammary gland development in the early growth stage of Holstein, we comprehensively analyzed 11 RNA-Seq datasets that were generated from mammary gland parenchyma of Holsteins heifer calves (ca. 8 weeks old) and identified candidate lncRNAs including novel ones in the mammary gland using bioinformatics analysis. Gene Ontology (GO) enrichment revealed these lncRNAs might modulate mammary gland development via cis-regulating. Our findings provide a list of potential lncRNAs involved in regulating mammary gland development in the early growth stage of Holsteins, thus may contribute to improving yield and quality of the milk production, as well as expand our understanding of Holstein mammary gland.

Materials and methods

Ethics statement

No ethical approval is required for this study.

Data collection and reads quality check

In this study, 11 single-end RNA-Seq datasets (30–101 bp) generated by sequencing using an Illumina HiSeq 2500 platform were retrieved from the Sequence Read Archive (SRA) database (Table S1). As described in the original study by Loor and colleagues (Vailati-Riboni et al. 2018), those samples were subsets of Gene Expression Omnibus database GSE102435 and collected from mammary parenchyma of Holstein heifer calves (ca. 8 weeks old). The cattle (Bos taurus) reference genome (ARS-UCD1.2) and genome annotation file (GTF format file, release 96) were downloaded from the Ensembl website (http://asia.ensembl.org/). Known splice sites in the cattle genome were extracted from the genome annotation file using HISAT2 program (Kim et al. 2015). The quality of these retrieved reads in FASTQ format was evaluated using Fastqc program (Babraham Institute, http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). QC check results yielded from Fastqc showed that all reads were of high quality and did not contain adapters (Table S2). Those retrieved reads were then used for other downstream analyses.

Transcript reconstruction and lncRNA identification

An inhouse lncRNA analysis pipeline was designed to identify lncRNAs in the mammary gland (Figure 1). All clean reads for each of the samples were mapped back onto the cattle reference genome using read aligner HISAT2 (2.1.0) (Kim et al. 2015) with parameters “–known-splicesite-infile” and “–dta”. The transcripts were reconstructed for each sample by Stringtie program (v1.3.4d) (Pertea et al. 2016) with the parameters “-f 0.1“, “-a 10“, “-j 1“, and “-c 2.5“. Then, the 11 samples were merged using Stringtie with “-merge” parameter to generate a global, unified set of transcripts across multiple samples. The fragment per kilobase of transcript per million mapped reads (FPKM) were estimated for the merged transcripts using the Stringtie program (Pertea et al. 2016) with the parameters “-e” and “-B” (Pertea et al. 2016). To obtain new transcripts, merged transcripts were compared with the reference genome annotation file using Gffcompare software (v0.10.6, available online http://ccb.jhu.edu/software/stringtie/gffcompare.shtml), and only the unknown intergenic transcripts (marked with class code “u”) were subjected to further analyses (Trapnell et al. 2012). For the unknown intergenic transcripts, stringent conditions were performed to identify true lncRNAs by removing the transcripts that (1) had FPKM <0.5, (2) single exons, and (3) length size shorter than 200 bp. Then, (4) for the remaining transcripts, five softwares of CPC2 (Kong et al. 2007), CNCI (Sun et al. 2013), CPAT-H (model trained by human lncRNA) (Wang et al. 2013), CPAT_M (model trained by mouse lncRNA) (Wang et al. 2013), and Pfam (Finn et al. 2014) were employed to discriminate non-coding transcripts from the coding ones by calculating for coding ability. Further, (5) we selected the intersection of the non-coding transcripts yield from those five softwares (CPC2, CNCI, CPAT-H, CPAT_M, and Pfam) as the lncRNA set. Lastly, BLASTN (v2.8.1) was used to compare those predicted lncRNAs with NONCODEv5 database (Fang et al. 2018), and the standard for sequence similarity was defined with e-value of <0.000001. In addition, the sequence similarity of the identified lncRNAs was detected using software CD-HIT (Fu et al. 2012), and sequence identity threshold was set with default parameter.

Figure 1. Pipeline for systematic identification of lncRNAs in the mammary gland of Holstein

Potential target prediction and functional enrichment

Most studies suggest that lncRNA expression could act on neighboring protein-coding RNAs by cis-regulating way (Ren et al. 2016; Wang et al. 2016). In this study, the coding genes located at the lncRNA within the distance of 10 kb were considered potential cis-regulated target genes, which were extracted using BEDtools windowBed program (Quinlan 2014). Cis-regulating networks that constructed based on lncRNA–target pairs were screened out using Cytoscape software (Shannon et al. 2003). The prediction of lncRNA functions was based on functional annotation of the potential cis-regulated target genes. GO enrichment for the potential target genes was performed using the David database (Huang et al. 2009).

Validation of lncRNAs by RT-PCR

The RT-PCR primers in the validation of lncRNAs were designed by online software Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/) (Table S3). Total RNA was extracted from bovine mammary epithelial cells (MAC-T) by Trizol reagent (Invitrogen, Hong Kong, China) and was converted into cDNA using PrimeScripts RT Reagent Kit containing gDNA Eraser (TAKARA, Dalian, China). PCR was performed using 2× Taq PCR Master Mix (TIANGEN, Beijing, China) and amplification reaction program used 94°C for 10 s, followed by 35 cycles of 94°C for 30 s, annealing temperature for 30 s and 72°C for 30 s, and final extension for 10 min at 72°C. PCR product was analyzed using 1.5% agarose gel.

Results

In this study, a total of 292,855,791 reads were retained from those originally retrieved from the SRA datasets after a quality check. About 96.84% of reads were successfully aligned back to the reference genome for each sample. An average of 35,975 transcripts (34,124–39,211) were reconstructed from all 11 samples (Table S1). After comparing with the reference annotation using Gffcompare program, a total of 1151 unknown intergenic transcripts were obtained. Most of the unknown intergenic transcripts contained 2–4 exons and were shorter than 1000 bp. Comparing with other chromosomes, chromosome 21 contained more unknown intergenic transcripts and was observed to be slightly different in unknown intergenic transcript length distribution and exon content (Figure 2).

Figure 2. Unknown intergenic transcripts in the mammary gland of Holstein heifer calves

LncRNA-sets contained 767, 759, 714, 739, and 723 lncRNA transcripts that were predicted by CPC2, CNCI, CPAT-H, CPAT-M, and Pfam, respectively. The intersection of these lncRNA-sets yielded 588 lncRNA transcripts (Figure 3(a) and Table S4). Those lncRNAs were further classified into lincRNA, which were located on the intergenic region of the cattle genome. The distribution of 588 lncRNAs on each chromosome is shown in Figure 3(b). The ORF length size of most lncRNA transcripts was less than 100 (Figure 3(c)). The length size of the identified lncRNAs ranged from 203 to 12,421 bp, most of which ranged from 300 bp to 900 bp (Figure 3(d)). In addition, the distance between identified autosomal lncRNAs ranged from 114 bp to 56 M bp, showing a large span in distance. Sequence similarity analysis showed that the 558 lncRNAs were classified into 557 clusters and were very different in sequence structure. After comparing lncRNA sequences with NONCODEv5 database using BLASTN (v2.8.1), 347 out of the 588 lncRNAs were aligned to the database significantly (e-value <0.000001) (Table S5), and 241 novel lncRNAs were detected.

Figure 3. Identification and characterization of lncRNAs in the mammary gland of Holstein heifer calves. (a) A Venn diagram showing the lncRNA transcripts generated by CPC2, CPAT-H, CPAT-M, CNCI, and Pfam. The lncRNA identified by all five softwares was used for downstream analyses. (b) Distribution of lncRNAs on each chromosome. (c) ORF length size distribution of lncRNA transcripts. (d) Length distribution of lncRNA transcripts

To obtain potential cis-regulated targets of identified lncRNAs, we explored the cis-regulated targets based on the different thresholds of ±100 kb, ±50 kb, and ±10kb from protein-coding gene to lncRNA. Figure 4(a) shows that approximately 2–6, 1–4, and 1–2 potential genes were targeted by each of the lncRNAs in the three respective different threshold conditions. Figure 4(b) shows that approximately 1–2, 1, and 1 lncRNAs were targeting each of the protein-coding genes in the three respective different threshold conditions. To increase the accuracy in the prediction of cis-regulated target gene, the ±10 kb threshold was selected as the filtering condition. With this threshold, for 364 out of 588 lncRNAs, we obtained 417 potential protein-coding RNAs, and two major lncRNA–target networks are shown in Figure 4(c,d). In addition, for the remaining 224 lncRNAs, we have not obtained their cis-regulated target on their ±10 kb window.

Figure 4. Prediction of potential cis-regulated target genes of lncRNAs. (a) Distribution of potential cis-regulated target gene number for each of the lncRNAs in different threshold conditions. (b) Distribution of the number of lncRNAs that targeting each of the coding genes in different threshold conditions. (c) and (d) Two major cis-regulating networks. The color of orange represents protein-coding genes, and the color of green represents lncRNAs

GO enrichment analysis further showed that the 417 potential cis-regulated targets were significantly enriched in 157 GO terms (p < 0.05), of which cellular respiration (GO:0045333), respiratory electron transport chain (GO:0022904), and electron transport chain (GO:0022900) were the top three listed GO terms in the biological process (BP) category; respiratory chain (GO:0070469), respiratory chain complex (GO:0098803), and mitochondrial respiratory chain (GO:0005746) were the top three listed GO terms in the cellular component (CC) category; and NADH dehydrogenase (ubiquinone) activity (GO:0008137), NADH dehydrogenase (quinone) activity (GO:0050136), and NADH dehydrogenase activity (GO:0003954) were the top three listed GO terms in the molecular function (MF) category (Figure 5). In addition, cell development (GO:0048468), mesenchymal cell development (GO:0014031), and regulation of cell differentiation (GO:0060284) were also significantly enriched (Table S6).

Figure 5. GO enrichment analysis for potential cis-regulated target genes of lncRNAs. Histogram shows the top ten significantly enriched terms in BP, CC, and MF (p < 0.05)

To confirm whether the lncRNAs identified by our identification pipeline were indeed expressed, we randomly selected 7 lncRNAs from overall 588 identified lncRNAs and examined them by RT-PCR. The RT-PCR result showed that the RT-PCR products of all seven lncRNAs (MSTRG.1772.2, MSTRG.8240.1, MSTRG.3052.1, MSTRG.229.1, MSTRG.8927.2, MSTRG.1214.1, and MSTRG.65.1) were amplified successfully and showed the expected size (Figure 6). Therefore, the lncRNAs we identified in this study were reliable.

Figure 6. Validation of seven randomly selected lncRNAs by RT-PCR. Seven randomly selected lncRNAs numbered from 1 to 7 were described in detail in Table S3

Discussion

Transcripts including protein-coding RNAs and ncRNAs have been widely discovered or further studied by RNA-Seq (Ibeagha-Awemu et al. 2018; Li et al. 2018; Yun et al. 2018). The ability of discovering new transcripts is a major advantage of RNA-Seq technology because of the global and systematic summary of whole transcripts for biological samples (Chen et al. 2017). In this study, 1151 unknown intergenic transcripts were assembled by reconstruction from a transcriptome. Credible lncRNAs, not including false-positive ones, were obtained by calculating for the protein-coding ability of the unknown intergenic transcripts using different models from five different softwares including CPC2, CNCI, CPAT-H, CPAT-M, and Pfam. The intersection of non-coding results generated from these five softwares was considered as the lncRNAs. Our data showed that most ORFs of Holstein lncRNAs had a length size less than 100, consistent with previous studies that have reported short length size of ORFs for lncRNAs in some species of rabbit and mouse (Ruizorera et al. 2014; Wang et al. 2018). The length size of the most lncRNAs ranged from 300 bp to 900 bp, corroborating those reported by Tong et al. in adult Holstein (Tong et al. 2017).

In the past, lncRNAs were considered as a class of by-products of genomic transcription and did not possess any biological function (Doolittle 2013). However, in recent years, more and more studies reported that lncRNAs were widely involved in biological processes such as DNA methylation, histone modification, and chromatin remodeling (Rinn & Chang 2012) and regulated the expression of target genes at the transcriptional, post-transcriptional, and epigenetic levels (Chen & Carmichael 2010; Wierzbicki 2012; Bai et al. 2015). In this study, the results of GO enrichment showed that the lncRNAs were associated with energy metabolism and cell differentiation. Energy metabolism is the process of generating energy (ATP) from diet and closely related to organ development. In Holsteins, energy metabolism mainly depends on the composition of the diet of an individual. The development of the mammary gland in Holstein heifers was reported to be affected by different dietary metabolizable proteins, metabolizable energy ratios (Albino et al. 2015), and enhanced by high plane of nutrition (Vailati-Riboni et al. 2018). The enriched GO terms of cellular respiration, respiratory electron transport chain, and electron transport chain in the present study indicate that lncRNAs might regulate energy metabolism in the early growth stages of Holstein mammary gland. Cell proliferation and differentiation are also required for organ development, and recent studies have shown that lncRNAs such as H19 (Li et al. 2019), Neat1 (Standaert et al. 2014), and XIST (Ma et al. 2019) were positively involved in growth and proliferation of mammary epithelial cells. Functional prediction of lncRNAs in the present study showed that some biological processes such as cell development, mesenchymal cell development, and regulation of cell differentiation were significantly enriched, suggesting that lncRNAs might modulate mammary gland development by regulating cell differentiation in the early growth stage of Holstein. However, the exact functions of these lncRNAs further need more focused investigation to verify our observations.

Conclusion

In conclusion, lncRNAs in Holstein mammary gland were analyzed using bioinformatics analysis. A list of credible lncRNAs including novel ones was obtained via bioinformatics analysis and RT-PCR validation. Prediction of potential cis-regulated target gene revealed that lncRNAs might regulate mammary gland development by regulating energy metabolism and cell differentiation. Since mammary gland development was closely related to the milk production, lncRNAs identified in this study might also be the molecular markers underlying the trait of milk production. This study provides valuable insights and could be the basis for further investigations to determine the functions of lncRNAs in the mammary gland development in Holstein. It also generally contributes to a better understanding of Holstein mammary gland development.

Disclosure statement

No potential conflict of interest was reported by the authors.

Supplemental data

Supplemental data for this article can be accessed here.

Additional information

Funding

This work was supported by key technologies and applications of breeding management for producing high-quality milk sources [2018NZ0003].