Next Article in Journal
Presence of Helicobacter pylori and H. suis DNA in Free-Range Wild Boars
Next Article in Special Issue
Detection of Genomic Regions with Pleiotropic Effects for Growth and Carcass Quality Traits in the Rubia Gallega Cattle Breed
Previous Article in Journal
A Peculiar Distribution of the Emerging Nematode Angiostrongylus cantonensis in the Canary Islands (Spain): Recent Introduction or Isolation Effect?
Previous Article in Special Issue
Impact of Censored or Penalized Data in the Genetic Evaluation of Two Longevity Indicator Traits Using Random Regression Models in North American Angus Cattle
 
 
Search for Articles:
Title / Keyword
Author / Affiliation / Email
Journal
Article Type
 
 
Advanced Search
 
Section
Special Issue
Volume
Issue
Number
Page
 
Journals
Animals
Volume 11
Issue 5
10.3390/ani11051268
animals-logo
Submit to this Journal Review for this Journal Propose a Special Issue
► Article Menu

Article Menu

share Share announcement Help format_quote Cite question_answer Discuss in SciProfiles thumb_up
...
Endorse textsms
...
Comment
first_page
settings
Order Article Reprints
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Long Intergenic Non-Coding RNAs in the Mammary Parenchyma and Fat Pad of Pre-Weaning Heifer Calves: Identification and Functional Analysis

by
Shengchao Zhang
1,
Sibtain Ahmad
2,
Yuxia Zhang
1,
Guohua Hua
1,3,* and
Jianming Yi
1,3,*
1
Key Laboratory of Agricultural Animal Genetic, Breeding, and Reproduction for Ministry of Education, College of Animal Science, Huazhong Agricultural University, Wuhan 430070, China
2
Faisalabad-Sub Campus Depalpur, University of Agriculture, Okara 56150, Pakistan
3
National Center for International Research on Animal Genetics, Breeding and Reproduction (NCIRAGBR), Huazhong Agricultural University, Wuhan 430070, China
*
Authors to whom correspondence should be addressed.
Animals 2021, 11(5), 1268; https://doi.org/10.3390/ani11051268
Submission received: 17 March 2021 / Revised: 10 April 2021 / Accepted: 24 April 2021 / Published: 28 April 2021
(This article belongs to the Collection Advances in Cattle Breeding, Genetics and Genomics)
Browse Figures
Versions Notes

Abstract

:

Simple Summary

The development of mammary gland is directly related to the productivity of dairy animals. Some studies showed that feeding the enhanced plane of nutrition at pre-weaning stage are advantageous to the development of mammary gland. However, regulators which are involved in this biological process remain largely unknown. In this work, we have identified some long intergenic non-coding RNAs (lincRNAs) in mammary parenchyma (PAR) and mammary fat pad (MFP) of heifer calves under different levels of nutrition at pre-weaning stage by using the published RNA-seq database. Furthermore, those putative lincRNAs, which were highly correlated with these key protein-coding genes in mammary gland development, were highlighted. Our results not only confirmed the advantages of feeding calves with enhanced feeding plane in pre-weaning stage, but also provided fundamental base for further research on the biological processes of mammary gland development.

Abstract

Enhanced plane of nutrition at pre-weaning stage can promote the development of mammary gland especially heifer calves. Although several genes are involved in this process, long intergenic non-coding RNAs (lincRNAs) are regarded as key regulators in the regulated network and are still largely unknown. We identified and characterized 534 putative lincRNAs based on the published RNA-seq data, including heifer calves in two groups: fed enhanced milk replacer (EH, 1.13 kg/day, including 28% crude protein, 25% fat) group and fed restricted milk replacer (R, 0.45 kg/day, including 20% crude protein, 20% fat) group. Sub-samples from the mammary parenchyma (PAR) and mammary fat pad (MFP) were harvested from heifer calves. According to the information of these lincRNAs’ quantitative trait loci (QTLs), the neighboring and co-expression genes were used to predict their function. By comparing EH vs R, 79 lincRNAs (61 upregulated, 18 downregulated) and 86 lincRNAs (54 upregulated, 32 downregulated) were differentially expressed in MFP and PAR, respectively. In MFP, some differentially expressed lincRNAs (DELs) are involved in lipid metabolism pathways, while, in PAR, among of DELs are involved in cell proliferation pathways. Taken together, this study explored the potential regulatory mechanism of lincRNAs in the mammary gland development of calves under different planes of nutrition.

1. Introduction

The mammary gland is a complex organ, distinguishing mammals from other animal species, undergoing through a series of developmental changes at different physiological stages starting from embryo to pubertal stage following reproductive stage, also known as mammogenesis. It is composed of different types of cells like parenchyma cells, mammary fat pad (MFP), fibroblasts and vascular endothelial cells [1]. Parenchyma (PAR) is the key tissue in synthesis and secretion of milk, while, the mammary fat pad (MFP) tissue is used to provide the protection and support for PAR [2]. MFP is necessary for development of the secretory epithelium and provides signals that mediate ductal morphogenesis and potentially alveolar differentiation. Moreover, these MFPs are susceptible to dietary changes [3]. Although, mammary gland ducts’ elongation and branching mainly occur in the pubertal period. But, the pre-weaning stage of mammary gland played an important role in heifer’s future milk yield [4,5]. Previous studies showed that the EH and R nutritional levels have different impacts on the mammary gland mass and composition of PAR and MFP in the pre-weaning stage. A number of genes have been identified that are involved in mammary gland development under different dietary planes, like, EGF, FGF2, IGF-1, etc. [6]. Still, there is a dire need to identify novel genes and their interaction in different tissues of the mammary gland.
Long non-coding RNAs (lncRNAs) have been defined as transcripts of length ≥200 nt that lacks protein-coding potential [7]. According to the genomic location and context, lncRNAs are divided into five classes, including intergenic, sense, antisense, intronic and bidirectional non-coding RNAs and vast majority is long intergenic non-coding RNAs (lincRNAs) [8,9]. LincRNAs have diverse different features form messenger RNA (mRNA) and exercise functions such as chromatin modifications and transcriptional regulation in nucleus, and also implied in post-transcriptional regulation in cytoplasm [10]. Besides, some of the lincRNAs have been confirmed to be the key regulators and biomarkers in the development of mammary gland, like, nuclear-enriched abundant transcript 1 (NEAT1), pregnancy induced noncoding RNA (PINC) and zinc finger NFX1-Type containing 1 (znfx1) antisense RNA 1(ZFAS1), etc. [11,12,13]. However, those lincRNAs, which might be involved in mammary development under different nutrient supply in pre-weaning period of Holstein calves are yet to be known [14,15,16,17].
In this study, we identified the lincRNAs in MFP and PAR of pre-weaning Holstein heifers, under enhanced (EH) and restricted (R) plane of nutrition [18]. A total of 534 transcripts originating from 434 gene loci were identified as putative lincRNAs. The function of these putative lincRNAs’ was predicted by analyzing neighboring genes and significantly correlated differentially expressed genes (DEGs) in MFP or PAR under EH and R nutrition supply. The present study broadens the knowledge of lincRNA annotation in bovine as well as facilitates future research about the mammary gland development.

2. Materials and Methods

2.1. Experiment Design and Library Construction

The experiment was previously published by Vailati-Riboni et al. [18]. Briefly, 12 Holstein heifer calves (6.0 ± 2 d old, 39.0 ± 4.4 kg) were divided into two groups under the same forage and feeding management conditions. Feeding was started at the end of 4th week and both of treatments were reduced to 50% at 8th week, to induce weaning. During the trials, both milk replacers were fed in two equal portions twice daily at 06:00 and 17:00 h and calves were provided with drinking water supply. Total RNAs was extracted from the MFP and PAR after the calves were euthanized and their whole mammary glands were removed and dissected. The RNA-seq cDNA libraries were constructed using the Illumina TruSeq Stranded mRNA Sample Prep kit. The single-end read library construction following the manufacturer’s instructions with mRNA enrichment.

2.2. Databases

A total of 22 single-read RNA-seq data were downloaded from NCBI Sequence Read Archive (SRA) database. The Bos taurus UMD3.1.1 reference genome FASTA file and the Gene Transfer Format (GTF) file were downloaded from the ensembl website (ftp://ftp.ensembl.org/pub/release-98/fasta/bos_taurus/dna/ (accessed on 1 June 2014)) and (ftp://ftp.ensembl.org/pub/release-98/gtf/bos_taurus/ (accessed on 1 June 2014)). The UniRef90 (UniProt Reference Clusters) database was downloaded from the UniProt website (http://www.ebi.ac.uk/uniprot/database/download.html (accessed on 1 January 2019)). Moreover, non-redundant reference sequence (RefSeq) NR data was downloaded from (ftp://ftp.ncbi.nih.gov/blast/db/ (accessed on 1 January 2016)). The Homo sapiens and Mus musculus over.chain file (converting genome coordinates intermediate files) downloaded from (https://hgdownload.soe.ucsc.edu/goldenPath/bosTau8/liftOver/ (accessed on 1 October 2014)). The bed file downloaded from (http://asia.ensembl.org/biomart/martview/ (accessed on 1 June 2014)).

2.3. Alignment and Assembly of RNA-Seq Data

Reads were aligned to Bos taurus reference genome (UMD3.1.1) by Hisat2 (version2.1.0, Iowa State University, Ames, IA, USA) with default parameters [19]. Mapped reads were assembled and 22 assembled transcripts files (GTF format) of four groups were then merged into a non-redundant transcriptome by StringTie (version 1.3.5, Johns Hopkins University, Baltimore, MD, USA) [20].

2.4. Identification and Characterization of Putative LincRNAs

LincRNAs are the intergenic transcript which have been defined as transcribed non-coding RNAs ≥ 200 nucleotides. Based on this, our pipeline to identify lincRNAs has been shown in (Figure 1a).
We used non-redundant transcriptome to identify lincRNAs: (1) only the “u’’ class as candidate linRNAs by gffcompare (version 0.10.6, Johns Hopkins University, Baltimore, MD, USA), which represented the unknown, intergenic transcripts, were retained [21]; (2) transcripts were selected to do further analysis (having exon ≥ 2 and length ≥ 200 bp); (3) We used coding potential calculation (CPC) tool (version -0.9-r2, Tsinghua University, Beijing, China) to calculate the coding potential of transcripts in both strands, CPC < 0 were retained [22]; (4) we used command “transeq” and “hmmerscan” of HMMER (version 3.2.1, HHMI Janelia Fam Research Campus, Ashburn, VA, USA) tool to translated the retained transcript sequence into six possible protein sequences, which had a significant hit in the Pfam database (E-value < 1 × 10−5) were removed [23]; (5) we compared retained transcript sequences with NCBINR and UniRef90 database by BLAST(version v2.6.0 +, National Center for Biotechnology Information, Bethesda, MA, USA), with similarity to known proteins (E-value < 1 × 10−5) were removed [24]; (6) those transcripts, having fragments per kilobase of transcript per million mapped reads (FPKM) score more than 0.1 at least one sample were retained. In addition, we found that 234 out of 534 putative lincRNA sequences were highly similar to the known lncRNA by comparing the NONCODE database (http://www.noncode.org (accessed on 1 January 2018)) [25]. Kolomogorv–Smirnov (KS) test was used for the characterization of putative lincRNAs and protein-coding genes. Furthermore, UCSC website liftOver tools (http://genome.ucsc.edu/cgi-bin/hgLiftOver (accessed on 2 July 2015)) and the over.chain file were used to compare evolutionary conservation of putative lincRNAs and protein-coding genes.

2.5. Express Analysis mRNA and Putative LincRNAs in MFP and PAR

Gene expression levels were estimated based on read counts through feature count (version 1.6.4) software [26]. R package DEseq2 was used to conduct differential expression tests between the two groups. And it was considered differentially expressed, if the value of |log2(Fold-change)| ≥ 1 and padj ≤ 0.05 [27]. In addition, unpaired t-test was used to identify differentially expressed lincRNAs (DELs).

2.6. Neighboring Gene Identification and Correlation Analysis of DELs of MFP and PAR

We identified the neighboring gene (<100 kb) of putative lincRNAs by Bedtools (version 2.29.0, University of Virginia School of Medicine, Charlottesville, VA, USA) [28]. Pearson correlation coefficient as applied to calculate the correlation between DELs and neighboring gene or DEGs. In addition, only putative lincRNAs’ function was inferred according to the pathway analysis result of the genes, which were adjacent to these DELs or significantly correlated with DELs.

2.7. GO Ontology and Pathway Analysis

In order to predict these putative lincRNAs’ function, the KOBAS (version 3.0) (http://kobas.cbi.pku.edu.cn/kobas3/?t=1 (accessed on 5 July 2011)) was used in order to conduct the pathway analysis [29]. p-value of pathway and enrichGO less than 0.05 were considered statistically significant.

2.8. Quantification of LincRNAs through QRT-PCR

PCR primers for three randomly selected lincRNAs were designed by the oligo7 program (Supplement Table S1). Total RNA was extracted from each target tissue of bovine by using TRIzol reagent (Invitrogen, California, CA, USA) as per the instruction of manufacturer. The Glyceraldehyde-3-Phosphate Dehydrogenase (GAPDH) served as the endogenous control gene. QRT-PCR was performed with SYBR® premix Ex Taq II (Tli RNaseH Plus) (2×) (TAKARA, Kyoto, Japan) to assess the expression level of these three lincRNAs and the results were calculated using the 2−ΔΔCt [30].

3. Results

3.1. Identification and Characterization of the Putative LincRNAs

A total of 22 RNA-seq data were collected from Gene Expression Omnibus (GEO) database. These data represent all the transcripts of MFP and PAR tissues from pre-weaning Holstein heifer calves, which were treated with restricted milk replacer (R, 0.45 kg/day, including 20% crude protein, 20% fat), and enhanced milk replacer (EH, 1.13 kg/day, including 28% crude protein, 25% fat) to identify lincRNAs in these tissues related to energy supply [18]. By using HISAT2, approximately 579.8 of 600.1 million clean reads were mapped on the Bos taurus reference genome (UMD3.1.1) (Supplement Table S2). Then 97,541 non-redundant transcripts were reconstructed for each sample by merged command of the stringtie software. The results of gffcompare showed that 9131 transcripts were intergenic transcripts. Finally, 534 transcripts originated from 434 gene loci were identified as the putative lincRNA based on their protein coding abilities (Figure 1a) (Supplement Table S3). By comparing with the NONCODEv5_cow database (http://www.noncode.org (accessed on 6 January 2018)), we found that 234 out of 534 transcripts have high similarity with known lincRNAs (Supplement Table S4). The other 300 transcripts were considered as novel lincRNAs (Figure 1b). More than 30 lincRNAs were enriched in BTA 3 (Bos taurus autosome 3), followed by BTA 5, 2, 10 and 11 (Figure 1c).
There were many differences between lincRNAs and protein-coding genes, including different exon number, length of transcript, expression level and evolutionary conservation, etc. [31] to characterize putative lincRNAs, these were compared with those protein-coding transcripts, which have been annotated in Ensembl database. The characteristics of protein-coding genes and putative lincRNAs were revealed by Kolomogorv–Smirnov (KS) test. These results showed that the average transcript length of these putative lincRNAs was significantly shorter when compared with the protein-coding genes (mean value 967 bp vs. 2058 bp, respectively; p-value < 2.2 × 10−16). (Figure 1d). Meanwhile, we found that the average exon number of the putative lincRNAs was less than protein-coding genes (mean value 2.4 vs. 9, respectively; p-value < 2.2 × 10−16) (Figure 1e). Most of the putative lincRNAs had two exons. Fragments per kilobase of transcript per million (FPKM) analysis results showed that the expression level of putative lincRNAs was lower than protein-coding gene (mean value 0.77 vs. 31.57, respectively; p-value = 9.508 × 10−13) (Figure 1f). After converting genome coordinates, we found that 3.7% and 0.56% of exonic regions of putative lincRNAs in Bos taurus have orthologous regions in Homo sapiens and Mus musculus. And, 10.6% and 5.6% of exonic regions of protein-coding genes in Bos taurus have orthologous regions in Homo sapiens and Mus musculus. It can be implied that lincRNAs showed lower evolutionary conservation than protein-coding transcripts both in Homo sapiens and Mus musculus (Figure 1g). Taken together, these putative lincRNAs showed shorter length, less exon number, lower expression and evolutionary conservation than protein-coding transcripts.
Previous studies showed that lincRNAs were remarkably tissue-specific than mRNA [32]. To confirm these identified lincRNAs, we randomly picked four novel lincRNAs to compare the expression level in mammary gland and other tissues by using the quantitative reverse transcription polymerase chain reaction (QRT-PCR). The results showed that these lincRNAs were highly expressed in mammary gland (Figure 2a–d).

3.2. Expression Analysis and Functional Prediction of DELs in MFP

A total of 79 DELs (61 upregulated, and 18 downregulated; EH vs R) were identified in response to different energy levels in MFP by using the unpaired t-test (Figure 3a). LincRNAs have been found to act as cis-element to participate in the transcriptional regulation of neighboring genes (<100 KB) [33,34]. According to the position of the genome, 261 genes (including 217 protein-coding genes) were found near those DELs (Supplement Table S4). These genes were involved in glucose metabolism and DNA repair signaling pathway (Figure 3b, upper yellow panel). Three of these nearby genes were also identified as DEGs in MFP (Table 1).
Besides, lncRNA could also act as trans-elements to regulate distant genes [35,36]. Those genes having similar co-expression patterns could be used to predict the putative lincRNA function. Therefore, we calculated the correlation between DELs and DEGs in MFP by using the Pearson correlation coefficient method (Supplement Table S5). A total of 488 DEGs (including 15 transcription factors) were found to be negatively or positively correlated with these 79 differentially expressed lincRNAs (p-value < 0.05). These DEGs were enriched in metabolism and signal transmission-related signaling pathways (Figure 3b, lower blue panel). Among them, peroxisome proliferators-activated receptors (PPARs) signaling pathway having the highest log2(p-value). Further analysis matched these DEGs with DELs to a total of 11,536 pairs. For example, TCONS_00016005 was positively correlated with apolipoprotein C3 (APOC3) and angiopoietin-related protein 4 (ANGPTL4), which were both enriched in PPAR signaling pathway and cholesterol metabolism pathway (Figure 3c).

3.3. Expression Analysis and Functional Prediction of DELs in PAR

By using the same method, 86 differentially expressed lincRNAs (54 upregulated, 32 downregulated; EH vs R) in mammary parenchyma (PAR) were identified (Figure 4a). There were 281 genes (including 256 protein-coding genes) near these DELs (Supplement Table S4), including 19 DEGs (Table 1). These neighboring genes were mainly enriched in various metabolism pathways (Figure 4b, upper yellow panel). A total of 1453 DEGs (including 92 transcription factors) were found to be negatively or positively correlated with these 86 differentially expressed lincRNAs (Supplement Table S5). These DEGs were enriched in immune, cell proliferation, cell cycle and signal transduction pathways (Figure 4b, lower blue panel). Further analysis matched these DEGs with DELs to a total of 47,099 pairs. For example, TCONS_0062567 and TCONS_0091815 were positively and TCONS_0066624 was negatively correlated with Insulin Like Growth Factor I (IGF-I), which were enriched in MAPK signaling pathway, PI3K-Akt signaling pathway, Focal adhesion, Ras signaling pathway and epidermal growth factor receptor (EGFR) tyrosine kinase inhibitor resistance pathway and involved in the mammary gland development.

4. Discussion

Some studies showed that feeding heifer claves with restricted milk replacer was the right strategy to save rearing cost and increase milk production [37,38]. With the development of research, it has been confirmed that increased nutrient supply would benefit the development of mammary gland ultimately leading towards increased milk yield [39,40,41]. Furthermore, nutrient supply can affect the total protein, total DNA, total fat, protein and fat concentration changes in MFP and PAR. In order to reveal this molecular mechanism, a number of transcriptomic studies have been reported [42,43]. However, it is still largely unknown. Over the last decade, more and more studies showed that lncRNAs were the key regulators in various biological processes, and most of them were lincRNAs [10]. Although, PINC, NEAT1 and ZFAS1 were reported to be involved in mammary epithelial cell proliferation. Which lincRNAs are involved in the mammary gland development under different feeding regimes remained largely unknown.
In this study, we have identified 534 putative lincRNAs from 22 samples from two tissues (MFP and PAR) and two treatments (EH, 1.13 kg/day, including 28% crude protein, 25% fat; R, 0.45 kg/day, including 20% crude protein, 20% fat), by using published high throughput RNA-seq data [18]. Results of characterization for these putative lincRNAs in our study are consistent with previous reports [44,45] which were also confirmed by QRT-PCR and it was concluded that most of the putative lincRNAs were highly expressed in mammary tissues.
In MFP tissue, those DELs neighboring genes and significantly correlated DEGs were enriched in various metabolism and signal transmission processes. Such as, PPAR signaling pathway, Calcium signaling pathway and Cytokine-cytokine receptor interaction, etc. Especially, the PPARs signaling pathway played a key function in preadipocyte proliferation [38,39]. In addition, we also found that these neighboring genes and significantly correlated DEGs were associated with lipid metabolism, which played a functional role in adipocyte tissue like, angiopoietin-related protein 4 (ANGPTL4), apolipoprotein C3 (APOC3), etc. [46,47]. And, these DEGs were downregulated in MFP under the restricted supply. Therefore, we speculated that these DELs might regulate neighboring genes and are significantly correlated with DEGs to promote MFP development under the enhanced nutrients supply.
In PAR tissue, those DELs’ neighboring genes and significantly correlated DEGs were enriched in metabolism-related and cell proliferation-related pathways. Such as, PI3K-Akt signaling pathway, JAK-STAT signaling pathway and Ras signaling pathway, etc. Moreover, TCONS_0062567 was negatively correlated with Insulin Like Growth Factor I (IGF-I), while, TCONS_0091815 and TCONS_0066624 were positively correlated with IGF-I, which were considered as the important regulators in the lactation process [48]. Some other significantly correlated DEGs play indispensable role in mammary gland development like, Epidermal Growth Factor Receptor (EGFR) and Cyclin D2 (CCND2), etc. [49,50]. And, these DEGs were downregulated in PAR under the R supply suggesting that these DELs may promote the PAR development under the EH supply.
Mammary gland development is a complex process, including various tissues and cell types [1,51,52]. Both MFP and PAR play important roles in various stages of mammary gland development. In this study, only lincRNAs with ploy A tail were identified for this biological process. Functional prediction was accomplished for DELs highly correlated DEGs. However, further research can be conducted to explore specific functions and regulated mechanisms for these lincRNAs.

5. Conclusions

In summary, we identified 79 and 86 DELs in MFP and PAR of pre-weaning heifer calves under the enhanced and restricted nutrient supply. These putative lincRNAs may influence metabolism, cell proliferation and tissue interaction in MFP and PAR by positive or negative regulation to promote the mammary gland development and tissue communication under the enhanced nutrition supply. In this study, only those putative lincRNAs were studied which were highly correlated with these key protein-coding genes in mammary gland development. Moreover, it is suggested to assess the functional analysis of these putative lincRNAs by experiment.

Supplementary Materials

The following are available online at https://www.mdpi.com/article/10.3390/ani11051268/s1, Table S1: Primers for real-time PCR assay used in this study, Table S2: The number of reads and mapping ratio of each sample, Table S3: The information of putative lincRNAs, Table S4: The alignment result of putative lincRNAs, Table S5: The correlation between DELs and DEGs.

Author Contributions

Conceptualization, methodology, data analyzes and manuscript draft by S.Z., S.A. and Y.Z.; experimental design and manuscript revision by J.Y. and G.H. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by National Key Research and Development Program of China (No. 2016YFD0500507; 2017YFD0501903), and Fundamental Research Funds for the Central Universities (No.2662018PY091).

Institutional Review Board Statement

All experiments in our study were approved by the Ethics Committee of Huazhong Agricultural University. And, all methods in our study were also approved by the committee mentioned above and performed in accordance with the relevant guidelines and regulations (Approval No: YJM-201001).

Data Availability Statement

Data analyzed in this study were a re-analysis of published RNA-seq data, which are openly available at locations cited in the reference section [18]. Further documentation about data processing is available at [NCBI GEO DataSets] at [GSE102435].

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Macias, H.; Hinck, L. Mammary gland development. Wiley Interdiscip. Rev. Dev. Biol. 2012, 1, 533–557. [Google Scholar] [CrossRef] [Green Version]
  2. Hovey, R.C.; Aimo, L. Diverse and active roles for adipocytes during mammary gland growth and function. J. Mammary Gland. Biol. Neoplasia 2010, 15, 279–290. [Google Scholar] [CrossRef] [Green Version]
  3. Lohakare, J.D.; Sudekum, K.H.; Pattanaik, A.K. Nutrition-induced changes of growth from birth to first calving and its impact on mammary development and first-lactation Milk yield in dairy heifers: A review. Asian Australas. J. Anim. Sci. 2012, 25, 1338–1350. [Google Scholar] [CrossRef] [Green Version]
  4. Piantoni, P.; Bionaz, M.; Graugnard, D.E.; Daniels, K.M.; Everts, R.E.; Rodriguez-Zas, S.L.; Lewin, H.A.; Hurley, H.L.; Akers, M.; Loor, J.J. Functional and gene network analyses of transcriptional signatures characterizing pre-weaned bovine mammary parenchyma or fat pad uncovered novel inter-tissue signaling networks during development. BMC Genom. 2010, 11, 331. [Google Scholar] [CrossRef] [Green Version]
  5. Watson, C.J.; Khaled, W.T. Mammary development in the embryo and adult: A journey of morphogenesis and commitment. Development 2008, 135, 995–1003. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  6. Butner, J.D.; Chuang, Y.L.; Simbawa, E.; Al-Fhaid, A.S.; Mahmoud, S.R.; Cristini, V.; Wang, Z. A hybrid agent-based model of the developing mammary terminal end bud. J. Theor. Biol. 2016, 407, 259–270. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  7. Kertz, A.F.; Hill, T.M.; Quigley, J.D.; Heinrichs, A.J., 3rd; Linn, J.G.; Drackley, J.K. A 100-Year Review: Calf nutrition and management. J. Dairy Sci. 2017, 100, 10151–10172. [Google Scholar] [CrossRef]
  8. Kopp, F.; Mendell, J.T. Functional Classification and Experimental Dissection of Long Noncoding RNAs. Cell 2018, 172, 393–407. [Google Scholar] [CrossRef] [Green Version]
  9. Beermann, J.; Kirste, D.; Iwanov, K.; Lu, D.; Kleemiß, F.; Kumarswamy, R.; Schimmel, K.; Bär, C.; Thum, T. A large shRNA library approach identifies lncRNA Ntep as an essential regulator of cell proliferation. Cell Death Differ. 2017, 25, 307–318. [Google Scholar] [CrossRef]
  10. Guttman, M.; Garber, M.; Levin, J.Z.; Donaghey, J.; Robinson, J.; Adiconis, X.; Fan, L.; Koziol, M.J.; Gnirke, A.; Nusbaum, C.; et al. Ab initio reconstruction of cell type-specific transcriptomes in mouse reveals the conserved multi-exonic structure of lincRNAs. Nat. Biotechnol. 2010, 28, 503–510. [Google Scholar] [CrossRef]
  11. Liu, S.; Wang, Z.; Chen, D.; Zhang, B.; Tian, R.R.; Wu, J.; Zhang, Y.; Xu, K.; Yang, L.-M.; Cheng, C.; et al. Annotation and cluster analysis of spatiotemporal- and sex-related lncRNA expression in rhesus macaque brain. Genome Res. 2017, 27, 1608–1620. [Google Scholar] [CrossRef] [Green Version]
  12. Standaert, L.; Adriaens, C.; Radaelli, E.; Van Keymeulen, A.; Blanpain, C.; Hirose, T.; Nakagawa, S.; Marina, J.-C. The long noncoding RNA Neat1 is required for mammary gland development and lactation. RNA 2014, 20, 1844–1849. [Google Scholar] [CrossRef] [Green Version]
  13. Hansji, H.; Leung, E.Y.; Baguley, B.C.; Finlay, G.J.; Cameron-Smith, D.; Figueiredo, V.C.; Askarian-Amiri, M.E. ZFAS1: A long noncoding RNA associated with ribosomes in breast cancer cells. Biol. Direct. 2016, 11, 62. [Google Scholar] [CrossRef] [Green Version]
  14. Shore, A.N.; Kabotyanski, E.B.; Roarty, K.; Smith, M.A.; Zhang, Y.; Creighton, C.J.; Dinger, M.E.; Rosen, J.M. Pregnancy-induced noncoding RNA (PINC) associates with polycomb repressive complex 2 and regulates mammary epithelial differentiation. PLoS Genet. 2012, 8, e1002840. [Google Scholar] [CrossRef] [Green Version]
  15. Tong, C.; Chen, Q.; Zhao, L.; Ma, J.; Ibeagha-Awemu, E.M.; Zhao, X. Identification and characterization of long intergenic noncoding RNAs in bovine mammary glands. BMC Genom. 2017, 18, 468. [Google Scholar] [CrossRef]
  16. Ibeagha-Awemu, E.M.; Li, R.; Dudemaine, P.L.; Do, D.N.; Bissonnette, N. Transcriptome Analysis of Long Non-Coding RNA in the Bovine Mammary Gland Following Dietary Supplementation with Linseed Oil and Safflower Oil. Int. J. Mol. Sci. 2018, 19, 3610. [Google Scholar] [CrossRef] [Green Version]
  17. Yang, B.; Jiao, B.; Ge, W.; Zhang, X.; Wang, S.; Zhao, H.; Wang, X. Transcriptome sequencing to detect the potential role of long non-coding RNAs in bovine mammary gland during the dry and lactation period. BMC Genom. 2018, 19, 605. [Google Scholar] [CrossRef] [Green Version]
  18. Vailati-Riboni, M.; Bucktrout, R.E.; Zhan, S.; Geiger, A.; McCann, J.C.; Akers, R.M.; Loor, J.J. Higher plane of nutrition pre-weaning enhances Holstein calf mammary gland development through alterations in the parenchyma and fat pad transcriptome. BMC Genom. 2018, 19, 900. [Google Scholar] [CrossRef] [Green Version]
  19. Kim, D.; Langmead, B.; Salzberg, S.L. HISAT: A fast spliced aligner with low memory requirements. Nat. Methods 2015, 12, 357–360. [Google Scholar] [CrossRef] [Green Version]
  20. Pertea, M.; Pertea, G.M.; Antonescu, C.M.; Chang, T.C.; Mendell, J.T.; Salzberg, S.L. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat. Biotechnol. 2015, 33, 290–295. [Google Scholar] [CrossRef] [Green Version]
  21. Pertea, G.; Pertea, M. GFF Utilities: GffRead and GffCompare. F1000Research 2020, 9. [Google Scholar] [CrossRef]
  22. Kang, Y.J.; Yang, D.C.; Kong, L.; Hou, M.; Meng, Y.Q.; Wei, L.; Gao, G. CPC2: A fast and accurate coding potential calculator based on sequence intrinsic features. Nucleic Acids Res. 2017, 45, W12–W16. [Google Scholar] [CrossRef] [Green Version]
  23. El-Gebali, S.; Mistry, J.; Bateman, A.; Eddy, S.R.; Luciani, A.; Potter, S.C.; Qureshi, M.; Richardson, L.J.; Salazar, G.A.; Smart, A.; et al. The Pfam protein families database in 2019. Nucleic Acids Res. 2019, 47, D427–D432. [Google Scholar] [CrossRef]
  24. Johnson, M.; Zaretskaya, I.; Raytselis, Y.; Merezhuk, Y.; McGinnis, S.; Madden, T.L. NCBI BLAST: A better web interface. Nucleic Acids Res. 2008, 36, W5–W9. [Google Scholar] [CrossRef]
  25. Fang, S.; Zhang, L.; Guo, J.; Niu, Y.; Wu, Y.; Li, H.; Zhao, L.; Li, X.; Teng, X.; Sun, X.; et al. NONCODEV5: A comprehensive annotation database for long non-coding RNAs. Nucleic Acids Res. 2018, 46, D308–D314. [Google Scholar] [CrossRef]
  26. Liao, Y.; Smyth, G.K.; Shi, W. featureCounts: An efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 2014, 30, 923–930. [Google Scholar] [CrossRef] [Green Version]
  27. Love, M.I.; Huber, W.; Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014, 15, 550. [Google Scholar] [CrossRef] [Green Version]
  28. Quinlan, A.R.; Hall, I.M. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics 2010, 26, 841–842. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  29. Xie, C.; Mao, X.; Huang, J.; Ding, Y.; Wu, J.; Dong, S.; Kong, L.; Gao, G.; Li, C.-Y.; Wei, L. KOBAS 2.0: A web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 2011, 39, W316–W322. [Google Scholar] [CrossRef] [Green Version]
  30. Livak, K.J.; Schmittgen, T.D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 2001, 25, 402–408. [Google Scholar] [CrossRef]
  31. Ransohoff, J.D.; Wei, Y.; Khavari, P.A. The functions and unique features of long intergenic non-coding RNA. Nat. Rev. Mol. Cell Biol. 2018, 19, 143–157. [Google Scholar] [CrossRef] [PubMed]
  32. Statello, L.; Guo, C.J.; Chen, L.L.; Huarte, M. Gene regulation by long non-coding RNAs and its biological functions. Nat. Rev. Mol. Cell Biol. 2021, 22, 96–118. [Google Scholar] [CrossRef]
  33. Zhao, W.; Mu, Y.; Ma, L.; Wang, C.; Tang, Z.; Yang, S.; Zhou, R.; Hu, X.; Li, M.; Li, K. Systematic identification and characterization of long intergenic non-coding RNAs in fetal porcine skeletal muscle development. Sci. Rep. 2015, 5, 8957. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  34. Tan, J.Y.; Smith, A.A.T.; Ferreira da Silva, M.; Matthey-Doret, C.; Rueedi, R.; Sonmez, R.; Ding, D.; Kutalik, Z.; Bergmann, S.; Marques, A.C. cis-Acting Complex-Trait-Associated lincRNA Expression Correlates with Modulation of Chromosomal Architecture. Cell Rep. 2017, 18, 2280–2288. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  35. Xie, X.; Fu, Y.; Liu, J. Chemical reprogramming and transdifferentiation. Curr. Opin. Genet. Dev. 2017, 46, 104–113. [Google Scholar] [CrossRef] [PubMed]
  36. Casero, D.; Sandoval, S.; Seet, C.S.; Scholes, J.; Zhu, Y.; Ha, V.L.; Luong, A.; Parekh, C.; Crooks, G.M. Long non-coding RNA profiling of human lymphoid progenitor cells reveals transcriptional divergence of B cell and T cell lineages. Nat. Immunol. 2015, 16, 1282–1291. [Google Scholar] [CrossRef]
  37. Sejrsen, K.; Huber, J.T.; Tucker, H.A.; Akers, R.M. Influence of nutrition of mammary development in pre- and postpubertal heifers. J. Dairy Sci. 1982, 65, 793–800. [Google Scholar] [CrossRef]
  38. Anderson, K.L.; Nagaraja, T.G.; Morrill, J.L. Ruminal metabolic development in calves weaned conventionally or early. J. Dairy Sci. 1987, 70, 1000–1005. [Google Scholar] [CrossRef]
  39. Van Amburgh, M.E.; Galton, D.M.; Bauman, D.E.; Everett, R.W.; Fox, D.G.; Chase, L.E.; Erb, H.N. Effects of three prepubertal body growth rates on performance of Holstein heifers during first lactation. J. Dairy Sci. 1998, 81, 527–538. [Google Scholar] [CrossRef]
  40. Moallem, U.; Werner, D.; Lehrer, H.; Zachut, M.; Livshitz, L.; Yakoby, S.; Shamay, A. Long-term effects of ad libitum whole milk prior to weaning and prepubertal protein supplementation on skeletal growth rate and first-lactation milk production. J. Dairy Sci. 2010, 93, 2639–2650. [Google Scholar] [CrossRef] [PubMed]
  41. Soberon, F.; Raffrenato, E.; Everett, R.W.; Van Amburgh, M.E. Preweaning milk replacer intake and effects on long-term productivity of dairy calves. J. Dairy Sci. 2012, 95, 783–793. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  42. Hare, K.S.; Leal, L.N.; Romao, J.M.; Hooiveld, G.J.; Soberon, F.; Berends, H.; Van Amburgh, M.E.; Martin-Tereso, J.; Steele, M.A. Preweaning nutrient supply alters mammary gland transcriptome expression relating to morphology, lipid accumulation, DNA synthesis, and RNA expression in Holstein heifer calves. J. Dairy Sci. 2019, 102, 2618–2630. [Google Scholar] [CrossRef] [PubMed]
  43. Holliday, H.; Baker, L.A.; Junankar, S.R.; Clark, S.J.; Swarbrick, A. Epigenomics of mammary gland development. Breast Cancer Res. 2018, 20, 100. [Google Scholar] [CrossRef] [Green Version]
  44. Derrien, T.; Johnson, R.; Bussotti, G.; Tanzer, A.; Djebali, S.; Tilgner, H.; Guernec, G.; Martin, D.; Merkel, A.; Knowles, D.G.; et al. The GENCODE v7 catalog of human long noncoding RNAs: Analysis of their gene structure, evolution, and expression. Genome Res. 2012, 22, 1775–1789. [Google Scholar] [CrossRef] [Green Version]
  45. Fu, Y.; Xu, Z.Z.; Lu, Z.J.; Zhao, S.; Mathews, D.H. Discovery of Novel ncRNA Sequences in Multiple Genome Alignments on the Basis of Conserved and Stable Secondary Structures. PLoS ONE 2015, 10, e0130200. [Google Scholar] [CrossRef]
  46. Kim, H.K.; Youn, B.S.; Shin, M.S.; Namkoong, C.; Park, K.H.; Baik, J.H.; Kim, J.B.; Park, J.-Y.; Lee, K.U.; Kim, Y.-B.; et al. Hypothalamic Angptl4/Fiaf is a novel regulator of food intake and body weight. Diabetes 2010, 59, 2772–2780. [Google Scholar] [CrossRef] [Green Version]
  47. Nordestgaard, B.G.; Nicholls, S.J.; Langsted, A.; Ray, K.K.; Tybjærg-Hansen, A. Advances in lipid-lowering therapy through gene-silencing technologies. Nat. Rev. Cardiol. 2018, 15, 261–272. [Google Scholar] [CrossRef]
  48. Christopoulos, P.F.; Msaouel, P.; Koutsilieris, M. The role of the insulin-like growth factor-1 system in breast cancer. Mol. Cancer 2015. [Google Scholar] [CrossRef] [Green Version]
  49. Sigismund, S.; Avanzato, D.; Lanzetti, L. Emerging functions of the EGFR in cancer. Mol. Oncol. 2018, 12, 3–20. [Google Scholar] [CrossRef]
  50. Zhu, W.; Zhao, M.; Mattapally, S.; Chen, S.; Zhang, J. CCND2 Overexpression Enhances the Regenerative Potency of Human Induced Pluripotent Stem Cell-Derived Cardiomyocytes: Remuscularization of Injured Ventricle. Circ. Res. 2018, 122, 88–96. [Google Scholar] [CrossRef] [PubMed]
  51. Berryhill, G.E.; Trott, J.F.; Hovey, R.C. Mammary gland development--It’s not just about estrogen. J. Dairy Sci. 2016, 99, 875–883. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  52. Inman, J.L.; Robertson, C.; Mott, J.D.; Bissell, M.J. Mammary gland development: Cell fate specification, stem cells and the microenvironment. Development 2015, 142, 1028–1042. [Google Scholar] [CrossRef] [PubMed] [Green Version]
Animals 11 01268 g001 550
Figure 1. Identification and characterization analysis of the putative long intergenic non-coding RNAs (lincRNAs). (a) Pipeline for the identification of putative lincRNAs. (b) Venn diagram of known and novel lincRNAs. (c) Bos taurus chromosome distribution of putative lincRNAs. Blue: novel lincRNAs, gray: known lincRNAs. (UR: unknown region, included NKLS02001272.1, NKLS02002063.1, NKLS02002207.1, NKLS02002210.1, NKLS02001064.1, NKLS02001662.1, NKLS02002183.1). (d) Transcript length of putative lincRNAs and protein-coding gene (X axis represent the range of transcription length (bp); Y axis represent the portion of specific length of transcripts); (e) exon numbers of lincRNAs and protein-coding gene (X axis represent the range of exon number; Y axis represent the portion of specific number of exon); (f) the expression levels of lincRNAs and protein-coding gene; (g) evolutionary conservation between protein-coding genes and putative lincRNAs of Bos taurus (Bos) with Homo sapiens (Homo) and Mus musculus (Mus).
Figure 1. Identification and characterization analysis of the putative long intergenic non-coding RNAs (lincRNAs). (a) Pipeline for the identification of putative lincRNAs. (b) Venn diagram of known and novel lincRNAs. (c) Bos taurus chromosome distribution of putative lincRNAs. Blue: novel lincRNAs, gray: known lincRNAs. (UR: unknown region, included NKLS02001272.1, NKLS02002063.1, NKLS02002207.1, NKLS02002210.1, NKLS02001064.1, NKLS02001662.1, NKLS02002183.1). (d) Transcript length of putative lincRNAs and protein-coding gene (X axis represent the range of transcription length (bp); Y axis represent the portion of specific length of transcripts); (e) exon numbers of lincRNAs and protein-coding gene (X axis represent the range of exon number; Y axis represent the portion of specific number of exon); (f) the expression levels of lincRNAs and protein-coding gene; (g) evolutionary conservation between protein-coding genes and putative lincRNAs of Bos taurus (Bos) with Homo sapiens (Homo) and Mus musculus (Mus).
Animals 11 01268 g001
Animals 11 01268 g002 550
Figure 2. (a–d) Expression profile of four putative long intergenic non-coding RNAs (lincRNAs) in eight different tissues. Y axis represents relative expression level. Results are presented as mean values ± Standard error (SEM).
Figure 2. (a–d) Expression profile of four putative long intergenic non-coding RNAs (lincRNAs) in eight different tissues. Y axis represents relative expression level. Results are presented as mean values ± Standard error (SEM).
Animals 11 01268 g002
Animals 11 01268 g003 550
Figure 3. Functional analysis of differentially expressed lincRNAs (DELs) in mammary fat pad (MFP) under enhanced milk replacer (EH) vs restricted milk replacer (R). (a)The heat map of DELs in MFP under EH vs R. (b) The pathway analysis of DELs’ neighboring genes and significantly correlated genes (yellow, represent DELs’ neighboring genes; blue, represent significantly correlated differentially expressed genes (DEGs)) in MFP. (c) The selected DELs highly correlated with DEGs. (Arrow represent these genes were involved in those pathways)
Figure 3. Functional analysis of differentially expressed lincRNAs (DELs) in mammary fat pad (MFP) under enhanced milk replacer (EH) vs restricted milk replacer (R). (a)The heat map of DELs in MFP under EH vs R. (b) The pathway analysis of DELs’ neighboring genes and significantly correlated genes (yellow, represent DELs’ neighboring genes; blue, represent significantly correlated differentially expressed genes (DEGs)) in MFP. (c) The selected DELs highly correlated with DEGs. (Arrow represent these genes were involved in those pathways)
Animals 11 01268 g003
Animals 11 01268 g004 550
Figure 4. Functional analysis of differentially expressed lincRNAs (DELs) in mammary parenchyma (PAR) under enhanced milk replacer (EH) vs restricted milk replacer (R). (a) Heat map of DELs in PAR under EH vs R. (b) Pathway analysis of DELs’ neighboring genes and significantly correlated genes (yellow, represent DELs’ neighboring genes; blue, represent significantly correlated genes) in PAR. (c) Selected DELs that are highly correlated with DEGs. (Arrow represent the genes were involved in those pathways)
Figure 4. Functional analysis of differentially expressed lincRNAs (DELs) in mammary parenchyma (PAR) under enhanced milk replacer (EH) vs restricted milk replacer (R). (a) Heat map of DELs in PAR under EH vs R. (b) Pathway analysis of DELs’ neighboring genes and significantly correlated genes (yellow, represent DELs’ neighboring genes; blue, represent significantly correlated genes) in PAR. (c) Selected DELs that are highly correlated with DEGs. (Arrow represent the genes were involved in those pathways)
Animals 11 01268 g004
Table
Table 1. Differentially expressed lincRNAs’ (DELs’) neighbored DEGs in PAR and MFP.
Table 1. Differentially expressed lincRNAs’ (DELs’) neighbored DEGs in PAR and MFP.
DELs IDNeighboring Gene NameChange (EH vs R)Tissue
TCONS_00049867ACVR2BupMFP
TCONS_00055163LAMA1downMFP
TCONS_00027412LOC100847119upMFP
TCONS_00063960PYGMupPAR
TCONS_00033404HRCupPAR
TCONS_00043293THEMIS2downPAR
TCONS_00040017LIMS2upPAR
TCONS_00049846STACupPAR
TCONS_00091815PDE10AdownPAR
TCONS_00016752HCKdownPAR
TCONS_00090096NTRK2downPAR
TCONS_00039007MRC2downPAR
TCONS_00033404KCNA7downPAR
TCONS_00015641PTGISupPAR
TCONS_00015642PTGISupPAR
TCONS_00042229SCN7AupPAR
TCONS_00036160ACEupPAR
TCONS_00033404LHBdownPAR
TCONS_00040017MYO7BupPAR
TCONS_00072600GIMAP1downPAR
TCONS_00036008PYYupPAR
TCONS_00033404NTF4downPAR
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zhang, S.; Ahmad, S.; Zhang, Y.; Hua, G.; Yi, J. Long Intergenic Non-Coding RNAs in the Mammary Parenchyma and Fat Pad of Pre-Weaning Heifer Calves: Identification and Functional Analysis. Animals 2021, 11, 1268. https://doi.org/10.3390/ani11051268

AMA Style

Zhang S, Ahmad S, Zhang Y, Hua G, Yi J. Long Intergenic Non-Coding RNAs in the Mammary Parenchyma and Fat Pad of Pre-Weaning Heifer Calves: Identification and Functional Analysis. Animals. 2021; 11(5):1268. https://doi.org/10.3390/ani11051268

Chicago/Turabian Style

Zhang, Shengchao, Sibtain Ahmad, Yuxia Zhang, Guohua Hua, and Jianming Yi. 2021. "Long Intergenic Non-Coding RNAs in the Mammary Parenchyma and Fat Pad of Pre-Weaning Heifer Calves: Identification and Functional Analysis" Animals 11, no. 5: 1268. https://doi.org/10.3390/ani11051268

APA Style

Zhang, S., Ahmad, S., Zhang, Y., Hua, G., & Yi, J. (2021). Long Intergenic Non-Coding RNAs in the Mammary Parenchyma and Fat Pad of Pre-Weaning Heifer Calves: Identification and Functional Analysis. Animals, 11(5), 1268. https://doi.org/10.3390/ani11051268

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop