Callum MacPhillamy
Yan Ren
Tong Chen1
Stefan Hiendleder
Wai Yee Low- 1The Davies Research Centre, School of Animal and Veterinary Sciences, University of Adelaide, Roseworthy, SA, Australia
- 2Robinson Research Institute, The University of Adelaide, North Adelaide, SA, Australia
Introduction: MicroRNAs (miRNAs) play a crucial role in regulating gene expression during key developmental processes, including fetal development. Brahman (Bos taurus indicus) and Angus (Bos taurus taurus) cattle breeds represent two major cattle subspecies with strikingly different phenotypes.
Methods: We analyzed miRNA expression in liver samples of purebred and reciprocal crosses of Angus and Brahman to investigate breed and parent-of-origin effects at the onset of accelerated fetal growth.
Results: We identified eight novel miRNAs in fetal liver samples and 14 differentially expressed miRNAs (DEMs) between purebred samples. Correlation of gene expression modules and miRNAs by breed and parent-of-origin effects revealed an enrichment of genes associated with breed-specific differences in traits such as heat tolerance (Brahman) and fat deposition (Angus). We demonstrate that genes predicted to be targets of DEMs were more likely to be differentially expressed than non-targets (p-value < 0.05). We identified several miRNAs (bta-miR-187, bta-miR-216b, bta-miR-2284c, bta-miR-2285c, bta-miR-2285cp, bta-miR-2419-3p, bta-miR-2419-5p, and bta-miR-11984) that showed similar correlation patterns as bta-miR-2355-3p, which has been associated with the glutamatergic synapse pathway, a key facilitator of heat tolerance. Furthermore, we report Angus-breed-specific miRNAs (bta-miR-2313-5p, btamiR-490, bta-miR-2316, and bta-miR-11990) that may be involved in fat deposition. Finally, we showed that the DEMs identified in fetal liver are involved in Rap1, MAPK, and Ras signalling pathways, which are important for fetal development, muscle development and metabolic traits such as fat metabolism.
Conclusion: Our work sheds light on the miRNA expression patterns that contribute to gene expression differences driving phenotypic differences in indicine and taurine cattle.
Introduction
The two main lineages of modern cattle breeds are derived from two separate domestications of the wild auroch (Bos primigenius) (McTavish et al., 2013). The first domestication event occurred in the Fertile Crescent ∼10,000 years ago and gave rise to Bos taurus taurus (Bruford et al., 2003; Ajmone-Marsan et al., 2010; MacHugh et al., 2017). A second domestication from a different auroch subspecies occurred ∼1,500 years later in the Indus Valley (Loftus et al., 1994) and gave rise to Bos taurus indicus. The domesticated subspecies are referred to as taurine and indicine cattle; the Angus breed is a representative of taurine cattle, and the Brahman breed of indicine cattle (McTavish et al., 2013). The two breeds have starkly contrasting phenotypes; e.g., Angus cattle have been selected for growth traits and meat production (Elzo et al., 2012) and hence have shorter gestation length and mature earlier than Brahman cattle with a muscle fibre composition conducive to meat tenderness (Casas et al., 2011; Xiang et al., 2013). Brahman cattle, in contrast, have superior heat and disease tolerance traits and are thus adapted to tropical environments, but they mature more slowly and have coarser muscle fibres and tougher meat (Chase et al., 1997; Elzo et al., 2012; Dikmen et al., 2018; Goszczynski et al., 2018).
Although Angus and Brahman cattle differ by ∼1% genetically at the nucleotide level and have a different Y chromosome structure (Goldammer et al., 1997; Low et al., 2020), they produce fertile inter-subspecific hybrids (Hiendleder et al., 2008) that have enabled the investigation of genetic and epigenetic factors involved in phenotypic differences of quantitative traits (Hiendleder et al., 2008; Mao et al., 2008; Xiang et al., 2013; Xiang et al., 2014; Akanno et al., 2018; Gobena et al., 2018; Andrade et al., 2022). Genetic and epigenetic regulatory changes that underlie the phenotypic differences between taurine and indicine cattle are of substantial interest from both scientific and economic perspectives, i.e., understanding complex trait biology and using this knowledge to combine desirable traits from both subspecies (Li et al., 2023).
Differences in the expression of microRNA (miRNA), small non-coding RNAs of ∼22 bp essential for epigenetic regulation of gene expression (Eulalio et al., 2008), could be one of the main drivers of extensive transcriptomic and phenotypic differences observed between taurine and indicine cattle (Liu et al., 2021). However, little is known about miRNAs that differentiate the subspecies (Hanif et al., 2018; Dong et al., 2023). MiRNAs have been implicated in a wide range of biological processes, including growth and development (Awamleh et al., 2019; Morales-Roselló et al., 2022), cell differentiation (Galagali and Kim, 2020), metabolism (Rottiers and Näär, 2012) and response to environmental stimuli (Vrijens et al., 2015; Liu et al., 2020). In general, miRNAs are highly conserved across species (Macfarlane and Murphy, 2010) but frequently display species- and tissue-specific expression patterns (Jopling, 2012; Sun et al., 2014). Regulatory changes by miRNA affect gene expression and can contribute to the emergence of new traits and phenotypic diversity (Yan et al., 2013; Li et al., 2020; Hao et al., 2023).
Postnatal phenotype is, to a large extent, determined during fetal growth and development. For example, myogenesis and adipogenesis are crucial developmental processes that begin in utero and have far-reaching effects on the manifestation of postnatal traits such as muscle mass, as the number of muscle fibres is determined before birth (Picard et al., 2002) and fat deposition, which plays a prominent role in the tenderness and taste of meat (Du et al., 2015; Louveau et al., 2016). The liver is a key metabolic organ, and its role in energy metabolism (Lotto et al., 2023), tightly links it with many of the economically important traits like adipogenesis and myogenesis. Moreover, the fetal liver is responsible for allocating incoming energy into different developmental processes (Rui, 2014), like muscle growth and thus, any transcriptional differences observed may shed light on possible causes of phenotypic differences between the two breeds. Our previous studies in purebred and crossbred taurine and indicine fetuses have demonstrated extensive parent-of-origin effects on the fetal musculoskeletal system (Xiang et al., 2013; Xiang et al., 2014) and revealed non-mendelian parent-of-origin effects on gene expression profiles (Liu et al., 2021) that may in part be orchestrated by differences in miRNA expression. Here, we report differential miRNA expression by breed and parent-of-origin effects in the same fetal resource of purebred taurine Angus, indicine Brahman and reciprocal cross fetuses, explore correlations with mRNA expression, and identify potential biological pathways targeted by these miRNAs to shed light on molecular differences that contribute to the diverse taurine and indicine cattle phenotypes.
Materials and methods
Study animals and sample collection
Liver samples of concepti used in the present study were the same as those described in our previous work (Liu et al., 2021). All animal experiments and procedures described in this study complied with Australian guidelines, approved by the University of Adelaide’s Animal Ethics Committee and followed the ARRIVE Guidelines (https://arriveguidelines.org/) (Approval No. S-094-2005). Experimental concepti and samples were obtained as described in Xiang et al. (2013) and Liu et al. (2021). Briefly, the parents were purebred Angus (Bos taurus taurus) and purebred Brahman (Bos taurus indicus), denoted as BT and BI, respectively. Primiparous dams of ∼16–20 months of age were grazed on pasture supplemented with silage. The animals were inseminated with purebred Angus (BTBT) or purebred Brahman (BIBI) semen and pregnancy tested via ultrasound. Pregnant dams and fetuses were ethically sacrificed at day 153 ± 1 of gestation and fetuses were dissected to obtain liver tissue samples from the Lobus hepatis sinister. Samples were snap-frozen in liquid nitrogen and stored at −80°C until further use. The liver samples represented three female and three male individuals from each of the four genetic combinations: BT × BT, BT × BI, BI × BT, and BI × BI. The four genetic groups were denoted with the paternal breed listed first. For example, BIBT represents a fetus whose sire was Brahman (BI) and whose dam was Angus (BT).
microRNA extraction and sequencing
Total RNA was extracted from fetal liver tissue using Qiagen RNeasy® Plus Universal kit according to the manufacturer’s instruction. The microRNA libraries were made by using the Bioo Scientific® NEXTflex™ Small RNA-Seq kit v3 and sequenced at the Australian Cancer Research Foundation, Adelaide, Australia, using an Illumina® NextSeq 500. Individual sample names and their corresponding genetic group are outlined in Supplementary Table S1.
Sequence alignment
MicroRNA sequences were aligned using a custom Nextflow pipeline (Di Tommaso et al., 2017). Sequencing reads were first checked for quality using FastQC (v. 0.12.1) (Andrews, 2010). Adapters from the 5′ and 3′ ends of the reads were then trimmed using CutAdapt (v. 4.3) (Martin, 2011). We then filtered reads to retain those within the 17-28 bp range with a mean sequence quality of 25 using Prinseq (v. 0.20.4) (Schmieder and Edwards, 2011). Reads passing the filtering were reassessed with FastQC. Reads were then filtered for various bovine small RNA species that included ribosomal RNA (rRNA), transfer RNA (tRNA), small nuclear RNA (snRNA), and small nucleolar RNA (snoRNA). These small RNAs were downloaded from the Rfam database (Kalvari et al., 2018; Kalvari et al., 2020) via RNA Central and accessed on the 28th of April 2023. Next, we filtered reads against non-coding RNA (ncRNA) and coding DNA (cDNA) from Ensembl release 109 for ARS-UCD1.2. Reads that did not map to these RNA species were considered potential miRNAs and used in the subsequent analyses. Code relating to all analyses is available at: https://github.com/DaviesCentreInformatics/MicroRNA_BiVsBt. MiRNA sequencing reads are available from BioProject: PRJNA626458.
Messenger RNA from RNA-seq reads of the same samples (BioProject PRJNA626458) were aligned to ARS-UCD1.2 (Rosen et al., 2020) using the parameters described in Liu et al. (2021). Reads were quantified using featureCounts from the Rsubread package (v. 2.14.2) (Liao et al., 2019) and the ARS-UCD1.2 Ensembl gene annotation version 109.
Quantification of known and discovery of novel miRNAs
Potential miRNA reads were used as input for the miRDeep2 pipeline (Friedländer et al., 2011). We first used “mapper.pl” from miRDeep2 with default parameters, except “-l,” which we set to 17. This step produced collapsed reads and alignments in the miRDeep2 “arf” format, which were then used in the miRDeep2 quantification and discovery steps. As input to “miRDeep2.pl,” we used the collapsed reads from the “mapper.pl” step, the ARS-UCD1.2 reference genome, the “.arf” file from the “mapper.pl” step, mature and hairpin miRNAs belonging to Bos taurus, which is denoted as bta from miRbase (Kozomara and Griffiths-Jones, 2011). We also used mature miRNAs from Capra hircus (chi) and Ovis ares (oar) in the miRDeep2 pipeline. We used a mirDeep2 score of ≥4, an estimated probability of being a true positive ≥70%, a significant Randfold p-value and the precursor location as an ID to identify novel miRNAs (Mukiibi et al., 2020).
Identification of differentially expressed miRNAs and mRNAs
We used the output from miRDeep2 to generate counts for differentially expressed miRNA (DEM) analysis. We followed a standard differential expression workflow using DESeq2 (Love et al., 2014). All miRNA samples were sequenced in a single batch. The model parameters were breed and sex, with the contrasts being the six pairwise comparisons of each breed and the comparison between males and females.
The mRNA counts were quantified using featureCounts, and the output was used as input to DESeq2 for differential expression analyses (Love et al., 2014). Here, the model parameters were breed, sex and batch. The contrasts were the same as those used in the miRNA analyses. Only DEGs and DEMs with an adjusted p-value less than 0.05 were considered significant.
miRNA target prediction
The miRanda (v. 3.3a) (Enright et al., 2003) miRNA target prediction software was used to identify possible gene targets of miRNAs. We extracted the three prime untranslated regions (3′UTRs) from the ARS-UCD1.2 Ensembl release 109. We then extracted the sequences of all mature miRNAs of interest, e.g., DEMs and group-specific miRNAs. We then aligned mature miRNA sequences to the 3′UTRs of ARS-UCD1.2 using miRanda with default parameters (Enright et al., 2003).
miRNA and mRNA co-expression
The weighted gene co-expression network analysis (WGCNA) package (v. 1.72) was used to identify miRNA and mRNA co-expression networks (Langfelder and Horvath, 2008). WGCNA enables users to identify which genes have similar expression profiles (likely co-expressed) and to correlate that gene expression with other data like phenotype data or miRNA expression.
WGCNA removes genes with insufficient counts across samples, outlier genes, and outlier samples from the normalized count matrix. We used a variance stabilising transformation to normalise the matrix, as recommended by Langfelder and Horvath (2008). The input matrix is of the form m × n where m is the number of samples and n is the number of genes.
We then constructed the gene co-expression network by calculating an adjacency matrix from the filtered and normalized count matrix. The adjacency matrix is an n × n matrix, where n is the number of genes in the count matrix. The adjacency matrix is populated with values between 0 and 1, such that aij gives the connection strength between gene i and gene j. We used a “softPower” of five and seven in calculating the adjacency matrix for the parent of origin and breed comparisons, respectively. We then transformed the adjacency matrix into a topological overlap matrix (TOM) to calculate which genes have high topological overlap, i.e., are connected to roughly the same genes as one another. We next identified the dissimilarity TOM by subtracting the TOM from 1. Following this, we performed hierarchical clustering to identify genes that grouped into modules of co-expressed genes. Each module needed to contain at least 90 genes to be considered separate from another module, and modules needed a correlation of at least 0.8 to be merged. We then correlated the expression of these gene modules with each of the genetic groups and the corresponding miRNA expression data. We considered any correlation between miRNA, mRNA and genetic group with a p-value below 0.05 as significant. Significant modules were labelled first by a prefix denoting the comparison they were identified in and then with a letter. For example, module A, identified between the breeds, was labelled module br-A, module B, identified in the maternal comparison, was labelled module m-B, and module C, identified in the paternal comparison, was labelled module p-C.
DEGs in miRNA targets
To determine if the predicted targets of a given miRNA were more likely to be DE within the present study, we tested how likely the targets of each miRNA were to be DE over all DE genes. Here, we defined the number of trials (n) as the number of predicted targets for a given miRNA, I.e., there were n targets that could be DE or not. The successes (k) were the number of predicted targets that were also DE. Finally, the probability of success (p) was the number of DEGs in a given comparison divided by the total number of genes in the count matrix. We paired the DEM and DEG comparisons such that if we were comparing DEMs identified between BIBI and BTBT, we only considered DEGs that were also identified between BIBI and BTBT.
KEGG pathway analyses
The clusterProfiler R package (v. 4.8.1) (Yu et al., 2012) was used to perform KEGG pathway analysis. The “enrichKEGG” function was used to identify KEGG pathways significantly enriched in target genes. We used a cut-off of 0.05 for p-values, and the Benjamini-Hochberg method was used to adjust p-values to give q-values, where we also used a cut-off of 0.05 for significance testing. We performed KEGG analyses for the predicted targets of each miRNA that was either DE or specific to a particular group.
Determination of most frequently targeted pathways
To determine which pathways were most frequently targeted by either DE or group-specific miRNAs, we first performed target prediction for each miRNA. Next, we performed KEGG pathway enrichment with clusterProfiler (v. 4.8.1) (Yu et al., 2012) for the genes targeted by each miRNA. We then identified all unique pathways and counted the number of times they were targeted by a different miRNA. Word clouds to represent frequency of pathways were made using wordcloud (v. 1.9.2. http://amueller.github.io/word_cloud/) against Angus and Brahman cattle outline as background.
Results
miRNA and mRNA sequence quality
There was a sample average of ∼16.9 million reads before filtering (Supplementary Table S2). No reads were shorter than our minimum length threshold, but around 588,000 reads per sample were longer than our maximum length threshold, and around 779,000 reads were below our quality threshold and thus removed. This initial filtering left ∼15.5 million reads to be assessed as potential miRNAs. Just over 50% of reads were found to be other RNA species, not miRNAs (Figure 1A). Most non-miRNA RNA species were non-coding RNA (ncRNA), with an average of 7 million reads per sample mapped to known bovine ncRNA sequences (Figure 1A; Supplementary Table S2). After filtering, ∼8.2 million candidate miRNA reads remained for downstream analyses. Our previous work has reported the RNA-seq alignment quality in detail (Liu et al., 2021). Briefly, after read trimming, an average of 49.9 million reads remained for alignment, and between 72% and 77% of them could be assigned to genes.
FIGURE 1. (A) Bar chart showing mean proportions of sequence reads belonging to the filtered, cDNA, tRNA, ncRNA, rRNA, snoRNA, and snRNA, too short, too long and “low qual” categories. Filtered refers to reads that are not contaminants and passed length and quality thresholds. cDNA, tRNA, ncRNA, snoRNA, and snRNA refer to any reads that mapped to them, i.e., contaminants. (B) Upset plot outlining the intersections of miRNAs with a count of at least one in all samples within a group. The x-axis of the bar plot represents each combination of groups e.g., the first point refers to the group of miRNAs found in all groups. The second point refers to the miRNAs found in BIBI, BTBT and BTBI but not BIBT. BT denotes Angus, and BI denotes Brahman. The y-axis is the number of miRNAs found in each point on the x-axis. (C) PCA plot showing how all 24 samples cluster when using all expressed miRNAs. Samples are coloured by breed, where blue denotes BTBT, orange is BTBI, green is BIBT and red is BIBI. (D) Same as C but only using the top ten most variable miRNAs determined by variance. The colours are the same as in (C). (E) PCA plot showing clustering of samples using all expressed miRNAs. Dots represent an individual, and colours represent genetics. Blue dots represent BTBT, orange dots represent BTBI, green dots represent BIBT, and red dots represent BIBI. The x-axis represents PC1, and the y-axis represents PC2. (F) PCA plot showing clustering of samples using only the top ten most variable miRNAs based on variance. Colours and axes are the same as in (E).
Ten miRNAs can differentiate taurine and indicine breeds
We performed principal component analyses to determine if the samples clustered by genetic group using the miRNA data as they did with RNA-seq data (Liu et al., 2021). We observed no discernible clustering pattern when using all miRNAs (Figure 1E). We then filtered the miRNA data and used the ten most variable miRNAs, as measured by variance. We subsequently observed a clear demarcation between BTBT and BIBI samples (Figure 1F). The hybrid groups (BIBT; BTBI) showed considerable overlap with one another and the purebred samples. The BTBT samples tended to cluster more closely to other samples within their group than any other group with BIBT, BTBI, and BIBI samples spanning the full range of PC1 (Figure 1F).
Known miRNA expression and novel miRNA profiles
As the hybrid samples showed high similarity with both BIBI and BTBT genetic groups, we focused known miRNA expression on the purebred genetic groups. There were 414 and 416 known miRNAs with a normalized count of at least one in all BIBI and BTBT samples, respectively (Figure 1B; Supplementary Table S3). Most (388) known miRNAs were shared between BIBI and BTBT samples (Figure 1B). In the BTBT group, we found 28 miRNAs not expressed in the BIBI group. Conversely, we found 26 miRNAs in the BIBI group not expressed in the BTBT group (Figure 1B). In 21 of the 24 liver samples analysed, bta-miR-122 was the most highly expressed miRNA, with an average of 39% of all mapped reads. This liver-specific miRNA is conserved across species with pleiotropic functions, including in cholesterol, glucose and iron homeostasis and lipid metabolism (Thakral and Ghoshal, 2015).
When we separated fetal liver samples based on the breed of their maternal and paternal parents, there were between 376 and 398 miRNAs expressed in all 12 samples for each parental group, i.e., maternal BI, maternal BT, paternal BI and paternal BT (Supplementary Table S3; Figures 1C, D). We observed 361 miRNAs shared between the maternal BT and BI groups; a similar number was observed when we compared paternal BI and BT groups. The maternal BI group samples contained 37 miRNAs not expressed in the maternal BT group, while 15 miRNAs were unique to the maternal BT group when compared to the maternal BI group (Figure 1B). The inverse of this pattern was observed between paternal groups, where samples from the paternal BT group had 33 unique miRNAs compared to 16 in those from the paternal BI group (Figure 1C).
We discovered eight potentially novel miRNAs not listed in miRBase in BTBT and BIBI samples (Supplementary Table S4). We considered these high-confidence novel miRNAs as they were found in all samples within a particular group. Six of these novel miRNAs were shared between BIBI and BTBT samples, with one being exclusive to each group, i.e., one novel miRNA was found in all six BTBT samples but not BIBI samples and vice versa (Supplementary Table S4). We have included the sequences of these novel miRNAs in Supplementary Table S4. Given all predicted novel miRNAs were lowly abundant (less than 0.2% of filtered reads) and these are not validated, we did not analyse them further.
Gene expression modules correlated with breed
Using the WGCNA method (Langfelder and Horvath, 2008) and a significant correlation threshold of p < 0.05, we identified seven gene modules with significant correlations and one approaching (p = 0.09) our significance threshold (p < 0.05) among BIBI and BTBT groups (Figure 2A). These modules, A-H, comprised a total of 7,707 genes, with 239 of these being differentially expressed between BIBI and BTBT samples. We identified a total of 2,525 DEGs between BIBI and BTBT (Table 1), and together, these eight modules contained 239 of 2,525 (∼9%) DEGs identified between the two groups. Modules A, B, D and E displayed significant positive correlations with BTBT samples and significant negative correlations with BIBI samples. Modules F, G and H displayed significant negative correlations with BTBT samples and significant positive correlations with BIBI samples. Module C had a positive trend with BTBT samples and a negative trend with BIBI samples but was not statistically significant (p = 0.09).
FIGURE 2. (A) Correlation heat map displaying the correlations between traits and the mRNA expression data. The x-axis refers to each phenotypic trait, i.e., breed. The y-axis refers to each of the gene modules identified by WGCNA. Each cell in the heatmap is coloured by the strength of the correlation, where a darker shade of red denotes an increasingly positive correlation, and darker shades of blue denote an increasingly negative correlation. The values in the heatmap denote the p-value associated with the correlation. (B) Same as A but comparing maternal genetics. (C) Same as (A,B) but comparing paternal genetics. (D) Correlation heat map displaying the correlations between mRNA and miRNA. The x-axis denotes each of the miRNAs deemed of interest, i.e., was differentially expressed or only found in a particular group. The y-axis and correlation colour scheme are the same as (A). The matrix underneath the heatmap denotes which comparison or group within which a given miRNA was identified. A yellow box denotes a miRNA that was upregulated in BIBI, and a blue box denotes a miRNA that was upregulated in BTBT. BT denotes Angus, and BI denotes Brahman. A black square denotes the presence of that miRNA in that group. (E) Same as (D) but comparing miRNAs expressed in samples with different maternal genetics. A yellow box denotes miRNA that was upregulated in maternal BI. (F) Same as (D,E) but comparing samples with different paternal genetics. Yellow boxes denote miRNAs upregulated in paternal BI, and blue boxes denote miRNAs upregulated in paternal BT.
TABLE 1. Number of differentially expressed miRNAs and genes of breed and parent-of-origin effects.
We then performed KEGG pathway analysis on the genes within each module to determine if they were enriched in pathways that shed light on phenotypic differences between Brahman and Angus. Among the eight gene modules correlated with breed, four modules were enriched with pathways. Module D from the breed comparison (br-D) was enriched in eight pathways: amino sugar and nucleotide sugar metabolism, biosynthesis of nucleotide sugars, apoptosis, MAPK signalling, phosphatidylinositol signalling, one carbon pool by folate and lipid and atherosclerosis pathways. Module br-F displayed a single enrichment in the glutamatergic synapse pathway. Module br-G had enrichments in relaxin signalling, chemical carcinogenesis—reactive oxygen species, nicotinate and nicotinamide metabolism, thermogenesis, retrograde endocannabinoid signalling, GnRH signalling and purine metabolism pathways. Module br-C, while not significant, did approach a significant correlation with breed (p = 0.09) and displayed an enrichment of 32 pathways, including thermogenesis, TGF-beta signalling, mTOR signalling, neurotrophin signalling and Ras signalling (Supplementary Table S5). Modules br-A, br-B, br-E, and br-H were not enriched in any pathways.
Gene expression modules correlated with maternal genome
Three gene modules identified in the maternal comparison (m-A, m-B, and m-C) significantly correlated with maternal BI or BT groups (Figure 2B). These modules contained 214 DEGs, which was about 34% of the total differentially expressed genes (635) identified between maternal breed groups (Table 1). Modules m-A and m-C were positively correlated with the maternal BI group and negatively correlated with the maternal BT group. Module m-B was negatively correlated with the maternal BI group and positively correlated with the maternal BT group (Figure 2B). All three modules correlated with the maternal breed were enriched in pathways (Supplementary Table S5). Module m-A showed enrichment for the glutamatergic synapse and several other pathways (Supplementary Table S5). Pathways enriched in module m-B included ECM-receptor interaction and PI3K-Akt signalling pathways, among others (Supplementary Table S5). Module m-C displayed enrichments in cell cycle and DNA replication pathways, as well as several others (Supplementary Table S5).
Gene expression modules correlated with paternal genome
Using WGCNA, we identified four gene modules with significant correlations and one module approaching significance (p = 0.07) with paternal groups (Figure 2C). There were 13,805 genes in these modules, with 1,903 differentially expressed, i.e., ∼65% of the 2,097 DEGs identified between paternal groups were contained within these five gene modules. Modules A and B identified in the paternal group comparison (p-A, p-B) were positively correlated with paternal BI samples and negatively correlated with paternal BT samples. The remaining modules (pC-E) were negatively correlated with paternal BI and positively correlated with paternal BT samples (Figure 2C).
All five modules correlated with paternal breeds were enriched in pathways (Supplementary Table S5). Purine metabolism, bile secretion, glutamatergic synapse and pancreatic secretion were all enriched in module p-A. Module p-B was enriched in 84 pathways, too many to list; refer to (Supplementary Table S5). Module p-C saw, among others, an enrichment in MAPK signalling and CoA biosynthesis pathways (Supplementary Table S5). Module p-D was enriched for the spliceosome pathway. Module p-E was enriched for 20 pathways, including thermogenesis and purine metabolism.
miRNA-module correlations
There were 78 miRNAs identified as being either differentially expressed or only expressed in a breed, maternal breed or paternal breed group. Of these, 35 had a significant positive correlation with at least one module, and 27 had a significant negative correlation with at least one module (Supplementary Table S7). The greatest number of modules that were significantly correlated with a single miRNA was four. Bta-miR-2284c was positively correlated with modules br-F, br-G, br-H, m-A, and p-A and negatively correlated with br-B, br-D (Figures 2D–F; Table 2). Bta-miR-2285cp was positively correlated with br-F, br-G, br-H, m-A, and p-A, and negatively correlated with br-B, br-D and p-C (Figures 2D–F; Table 2).
TABLE 2. Positive and negative correlations between miRNAs and gene expression modules.
Bta-miR-2415-3p was positively correlated with br-D and br-E, while it also exhibited negative correlations with br-F, br-G, br-H, m-A, and p-A (Figures 2D–F; Table 2). Bta-miR-2440 was positively correlated with br-B, br-D and br-E and was negatively correlated with br-F, br-G, br-H, m-A, and p-A (Figures 2D–F; Table 2). Interestingly, bta-miR-2284c and bta-miR-2285cp had very similar module correlation patterns, the only difference being that bta-miR-2285cp was also negatively correlated with module p-C. Similarly, bta-miR-2415-3p and bta-miR-2440 had similar module correlation patterns, the only difference being that bta-miR-2440 was also positively correlated with module br-B. Furthermore, bta-miR-2284c and bta-miR-2285cp were positively correlated with the BIBI, maternal BIBI and paternal BIBI groups and displayed expression below our significance threshold in BTBT samples (Figures 2D–F).
Predicted targets of DEMs more likely to be DEGs
We compared the predicted targets of each DEM identified between BIBI and BTBT samples and determined whether they were more likely to be DE than the background genes, i.e., genes not predicted to be targeted by a given miRNA. We found that predicted targets of 12 of the 14 DEMs identified between BIBI and BTBT samples were significantly more likely to be DE than non-target genes (binomial test, p < 0.05) (Supplementary Table S8). Furthermore, the targets of the single DEM identified between maternal BI and maternal BT were significantly more likely to be DE (binomial test, p < 0.05). The targets of four of the five DEMs identified between paternal groups were significantly more likely to be DE than non-target genes (binomial test, p < 0.05).
Pathways predicted to be targeted by DE and group-specific miRNAs
There were 23, 15, 33, 29, 37, and 17 miRNAs upregulated or only expressed in breed BTBT, maternal BT, paternal BT, breed BIBI, maternal BI and paternal BI, respectively (Figure 3A; Supplementary Table S9). We used the miRanda target prediction software (v. 3.3a) (Enright et al., 2003) to identify potential targets of each miRNA (Figure 3). Between 300 and 5,115 targets were predicted among these miRNAs (Supplementary Table S10), which were then used to perform KEGG pathway analysis. Genes involved in the MAPK and Rap1 signalling pathways were consistently targeted by miRNAs, with between 37% and 56% of miRNAs in each group targeting these pathways (Supplementary Table S9). All pathways targeted by miRNAs expressed in each group can be seen in Figure 3; Supplementary Table S9.
FIGURE 3. Schematic overview of how the most frequently targeted pathways were identified. For each differentially expressed and group-specific miRNA, we predicted their mRNA targets. We then performed pathway enrichment, where the input was a list of genes predicted to be targeted by a given miRNA. We then identified all the unique pathways predicted to be targeted by all the upregulated and expressed miRNAs for a particular group (e.g., all upregulated and expressed miRNAs in the BTBT group) and counted how many miRNAs targeted each pathway. BT denotes Angus, and BI denotes Brahman. The number of miRNAs upregulated and expressed in each group was BIBI = 29, BTBT = 28, maternal BI = 37, maternal BT = 15, paternal BI = 17 and paternal BT = 33.
No DEMs were identified between the sexes
To determine potential differences in miRNA expression between the sexes, we compared all males against all females. We observed no DEMs, nor did we find significant correlations between gene modules and sex. Consequently, there was no pathway analysis performed. The lack of DEMs between males and females was consistent with the limited number of DEGs observed between male and female samples.
Discussion
To our knowledge, this study is the first to report on miRNA expression differences between taurine and indicine cattle and correlate the results with mRNA expression. We utilized miRNA expression data from fetal liver samples representing male and female individuals of Brahman, Angus and their reciprocal crosses. The inclusion of reciprocal crosses enabled us to disentangle not only miRNA expression differences that could be related to breed but also whether the maternal or paternal genomes impacted miRNA expression. Previous efforts to understand potential miRNA drivers of differences between the two cattle subspecies have focused on adult tissue. For example, Deb and Sengar (2021) investigated the miRNA expression profiles between Sahiwal (indicine) and Frieswal (indicine x taurine) cattle from blood samples in response to summer heat stress, identifying DEMs that interact with heat shock protein 70 (Hsp70). Similarly, Dong et al. (2023) investigated miRNA expression differences in the mature testes of Mongolian (taurine) and Hainan (indicine) cattle, demonstrating breed differences in spermatogenesis-related miRNAs. However, it is reasonable to assume that molecular drivers of physiological and phenotypic differences may not be detectable by the time the trait is fully developed in adults.
Our previous study on mRNA expression data in the same fetal resource clearly distinguished genetic groups BTBT, BIBT, BTBI, and BIBI (Liu et al., 2021). In contrast, miRNA data only separated purebred taurine and indicine fetal liver samples. Furthermore, despite observing thousands of differentially expressed genes between BTBT and BIBI samples, we identified only 14 DEMs between the two breeds. This result is similar to a study that observed only 23 DEMs when comparing miRNA expression in Hereford x Limousine beef cattle and Holstein-Friesian dairy cattle muscle cells during myogenic differentiation (Sadkowski et al., 2018). As a single miRNA targets multiple mRNAs (Peterson et al., 2014), the relatively low number of DEMs is not unexpected as these miRNAs may profoundly impact genes and gene networks differentially expressed between breeds.
Using the mRNA expression data, we identified gene modules correlated with the two breeds, Angus and Brahman, that represent taurine and indicine cattle. Moreover, by correlating the same gene expression data with miRNA expression, we identified miRNAs that may contribute to the gene expression differences that distinguish the two types of cattle. For example, module br-C was positively correlated with BTBT samples and was enriched for several pathways, including thermogenesis and Ras signalling. Ras signalling plays an essential role in adipogenesis (Murholm et al., 2010); adipocyte hyperplasia is evident in cattle at the mid-gestation developmental stage (Zhao et al., 2019). Interestingly, bta-miR-11990 was predicted to target genes involved in the Ras signalling pathway. This miRNA shares a similar correlation pattern with several BTBT-specific miRNAs (bta-miR-2313-5p, bta-miR-490, and bta-miR-2316), which were predicted to target pathways also involved in adipogenesis, such as NF-kappa B signalling (Peng et al., 2022) and VEGF signalling pathways (Park et al., 2017). Moreover, bta-miR-2316 has previously been reported as being expressed in the adipose tissue of cattle, suggesting a role in adipogenesis (Romao et al., 2014). However, further work is needed to confirm its role in adipogenesis in cattle. As Angus cattle are known to have superior fat deposition performance in cold climates compared to Brahman (Boyles and Riley, 1991), these miRNAs may be contributing to post-transcriptional regulation that conveys this trait to Angus cattle. However, the potential correlation of these miRNAs with fat deposition needs to be tested experimentally.
An advantage of including reciprocal crosses in breed comparisons is that it enables the investigation of maternal and paternal genome effects. We observed a single DEM between different maternal groups. In addition to this single DEM (bta-miR-187), we observed several miRNAs that were only expressed in maternal BI samples (bta-miR-2284c, bta-miR-2285c, bta-miR-2285cp, bta-miR-2285cz, bta-miR-2355-3p, and bta-miR-11984). These miRNAs were positively correlated with module m-A, a module that was significantly correlated with maternal BI samples and was enriched for the glutamatergic synapse pathway. Moreover, bta-miR-2355-3p was also predicted to target genes involved in this pathway. We observed a similar pattern of module-pathway associations in the sire of origin comparison, with bta-miR-2355-3p being positively correlated to module p-A, which was also enriched for the glutamatergic synapse pathway. Several miRNAs (bta-miR-187, bta-miR-216b, bta-miR-2284c, bta-miR-2285c, bta-miR-2285cp, bta-miR-2419-3p, bta-miR-2419-5p, and bta-miR-11984) shared a similar correlation pattern to bta-miR-2355-3p. The glutamatergic synapse has a known role in heat tolerance of an individual as glutamatergic neurons transmit peripheral and central heat signals to the hypothalamic preoptic area of the brain (Sun et al., 2022), which then begins a coordinated response to lower the temperature. The liver is known to play an integral role in coordinating this heat stress response via increased production of heat shock proteins, increasing metabolic rate and increased vasodilation (Thorne et al., 2020).
Further evidence to support the possible role of glutamatergic synapses in conveying heat tolerance in cattle can be found in another recent study of Holstein dairy cattle (Cheruiyot et al., 2021). Given the developmental timeline of the liver and the observed pathway enrichments (Tiniakos et al., 1996; Giancotti et al., 2019), it is possible that bta-miR-2355-3p and those miRNAs with similar correlation patterns play a role in modulating neuron development in the liver, priming it for hotter temperatures later in life and that this can be conveyed by either a Brahman sire or dam.
While we observed several DEMs that were correlated to gene modules, many more miRNAs were only expressed in one of the two groups in each comparison, e.g., bta-miR-490 in the breed comparison, bta-miR-11998 in the maternal breed comparison and bta-miR-187 in the paternal breed comparison. There was some overlap between DEMs and group-specific miRNAs, e.g., bta-miR-184 in the paternal-breed comparison. However, this was due to differences in the detection of expression threshold used to identify miRNAs as being expressed and how DESeq2 determines if a gene or miRNA should be retained. In any case, more miRNAs presented as only being expressed in one group, e.g., BIBI or BTBT, than displaying differential expression between groups. This pattern leads us to posit that group-specific miRNAs may be more important in driving gene regulatory differences than DEMs between Brahman and Angus cattle.
We identified a range of predicted targets for each of the differentially expressed and group-specific miRNAs and performed biological pathway analyses on these targets to gain insights into the pathways that may be affected. Notably, a substantial proportion of the targeted pathways identified in each of the three comparisons were signalling pathways, such as the Rap1, MAPK, and Ras signalling pathways. Each of these pathways plays an important role in fetal development. Rap1 signalling is important for vascular morphogenesis (Chrzanowska-Wodnicka, 2013). Additionally, the Rap1 signalling pathway has been implicated as a possible cause of the differences in meat production traits between high and low-performing meat goat breeds (Shen et al., 2022). Furthermore, Rap1 signalling ablation in the brain of mice has been shown to protect mice from high-fat diet-induced obesity (Kaneko et al., 2016), suggesting a critical role of Rap1 signalling in fat deposition. While this study investigated Rap1 signalling in the brain, given the close interactions between the brain and liver to monitor glucose and lipid homeostasis, it is possible that Rap1 signalling in the liver also has a role in fat storage. Ras and MAPK signalling pathways are critical to cell proliferation and differentiation (Zhang and Liu, 2002). Moreover, the Ras signalling pathway is purported to have a role in adipogenesis as ectopic expression of the pathway can induce preadipocyte formation in the absence of insulin and insulin-like growth factor 1 (Igf-1) (MacDougald and Lane, 1995), suggesting a possible role in modulating adipogenesis in the fetus. In addition, MAPK signalling can act as a negative regulator of muscle development (Xie et al., 2018), suggesting that differences in how this pathway is regulated may influence how the fetus develops muscle. The differences in miRNA expression between taurine and indicine cattle appear to result in differential modulation of signalling cascades involved in fetal liver development and growth regulation, which may contribute to phenotypic differences in traits such as heat tolerance, fat deposition and growth rate. In conclusion, we have identified differentially expressed miRNAs in the rapidly growing fetus that are involved in the genetic and epigenetic architecture of a range of complex traits, including fat deposition and heat tolerance, that differentiate taurine and indicine cattle, which are vital for global food supply.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.
Ethics statement
The animal study was approved by the University of Adelaide’s Animal Ethics Committee and followed the ARRIVE guidelines. The study was conducted in accordance with the local legislation and institutional requirements.
Author contributions
CM: Conceptualization, Data curation, Formal Analysis, Investigation, Methodology, Writing–original draft, Writing–review and editing. YR: Data curation, Methodology, Writing–original draft. TC: Investigation, Methodology, Resources, Writing–original draft. SH: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Writing–review and editing. WL: Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing–original draft, Writing–review and editing.
Funding
The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The study was funded by the JS Davies bequest through a Professorial Fellowship to SH and the Davies Livestock Research Centre. This work was supported with supercomputing resources provided by the Phoenix HPC service at the University of Adelaide.
Acknowledgments
We thank Prof John Williams for being a part of the team who initiated the miRNA work.
Conflict of interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Publisher’s note
All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.
Supplementary material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fgene.2023.1329939/full#supplementary-material
References
Ajmone-Marsan, P., Garcia, J. F., and Lenstra, J. A. (2010). On the origin of cattle: how aurochs became cattle and colonized the world. Evol. Anthropol. Issues, News, Rev. 19, 148–157. doi:10.1002/evan.20267
Akanno, E. C., Chen, L., Abo-Ismail, M. K., Crowley, J. J., Wang, Z., Li, C., et al. (2018). Genome-wide association scan for heterotic quantitative trait loci in multi-breed and crossbred beef cattle. Genet. Sel. Evol. 50, 48. doi:10.1186/s12711-018-0405-y
Andrade, A. N. N., Hernandez, A., Rodriguez, E. E., Davila, K. M. S., Mateescu, R., and Rodríguez, E. (2022). 38 effect of breed composition and genome-wide association study on epidermis thickness in a multibreed angus-brahman population. J. Animal Sci. 100, 9. doi:10.1093/jas/skac247.015
Andrews, S. (2010). FastQC. Available at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
Awamleh, Z., Gloor, G. B., and Han, V. K. M. (2019). Placental microRNAs in pregnancies with early onset intrauterine growth restriction and preeclampsia: potential impact on gene expression and pathophysiology. BMC Med. Genomics 12, 91. doi:10.1186/s12920-019-0548-x
Boyles, S. L., and Riley, J. G. (1991). Feedlot performance of Brahman x Angus versus Angus steers during cold weather. J. Anim. Sci. 69, 2677–2684. doi:10.2527/1991.6972677x
Bruford, M. W., Bradley, D. G., and Luikart, G. (2003). DNA markers reveal the complexity of livestock domestication. Nat. Rev. Genet. 4, 900–910. doi:10.1038/nrg1203
Casas, E., Thallman, R., and Cundiff, L. (2011). Birth and weaning traits in crossbred cattle from Hereford, Angus, brahman, boran, tuli, and Belgian blue sires. J. Animal Sci. 89, 979–987. doi:10.2527/jas.2010-3142
Chase, C. C., Chenoweth, P. J., Larsen, R. E., Olson, T. A., Hammond, A. C., Menchaca, M. A., et al. (1997). Growth and reproductive development from weaning through 20 months of age among breeds of bulls in subtropical Florida. Theriogenology 47, 723–745. doi:10.1016/s0093-691x(97)00030-7
Cheruiyot, E. K., Haile-Mariam, M., Cocks, B. G., Macleod, I. M., Xiang, R., and Pryce, J. E. (2021). New loci and neuronal pathways for resilience to heat stress in cattle. Sci. Rep. 11, 16619. doi:10.1038/s41598-021-95816-8
Chrzanowska-Wodnicka, M. (2013). Distinct functions for Rap1 signaling in vascular morphogenesis and dysfunction. Exp. Cell Res. 319, 2350–2359. doi:10.1016/j.yexcr.2013.07.022
Deb, R., and Sengar, G. S. (2021). Comparative miRNA signatures among Sahiwal and Frieswal cattle breeds during summer stress. 3 Biotech. 11, 79. doi:10.1007/s13205-020-02608-4
Dikmen, S., Mateescu, R. G., Elzo, M. A., and Hansen, P. J. (2018). Determination of the optimum contribution of Brahman genetics in an Angus-Brahman multibreed herd for regulation of body temperature during hot weather. J. Animal Sci. 96, 2175–2183. doi:10.1093/jas/sky133
Di Tommaso, P., Chatzou, M., Floden, E. W., Barja, P. P., Palumbo, E., and Notredame, C. (2017). Nextflow enables reproducible computational workflows. Nat. Biotechnol. 35, 316–319. doi:10.1038/nbt.3820
Dong, Z., Ning, Q., Liu, Y., Wang, S., Wang, F., Luo, X., et al. (2023). Comparative transcriptomics analysis of testicular miRNA from indicine and taurine cattle. Anim. Biotechnol. 34, 1436–1446. doi:10.1080/10495398.2022.2029466
Du, M., Wang, B., Fu, X., Yang, Q., and Zhu, M.-J. (2015). Fetal programming in meat production. Meat Sci. 109, 40–47. doi:10.1016/j.meatsci.2015.04.010
Elzo, M. A., Johnson, D. D., Wasdin, J. G., and Driver, J. D. (2012). Carcass and meat palatability breed differences and heterosis effects in an Angus–Brahman multibreed population. Meat Sci. 90, 87–92. doi:10.1016/j.meatsci.2011.06.010
Enright, A. J., John, B., Gaul, U., Tuschl, T., Sander, C., and Marks, D. S. (2003). MicroRNA targets in Drosophila. Genome Biol. 5, R1. doi:10.1186/gb-2003-5-1-r1
Eulalio, A., Huntzinger, E., and Izaurralde, E. (2008). Getting to the root of miRNA-mediated gene silencing. Cell 132, 9–14. doi:10.1016/j.cell.2007.12.024
Friedländer, M. R., Mackowiak, S. D., Li, N., Chen, W., and Rajewsky, N. (2011). miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades. Nucleic Acids Res. 40, 37–52. doi:10.1093/nar/gkr688
Galagali, H., and Kim, J. K. (2020). The multifaceted roles of microRNAs in differentiation. Curr. Opin. Cell Biol. 67, 118–140. doi:10.1016/j.ceb.2020.08.015
Giancotti, A., Monti, M., Nevi, L., Safarikia, S., D’Ambrosio, V., Brunelli, R., et al. (2019). Functions and the emerging role of the foetal liver into regenerative medicine. Cells 8, 914. doi:10.3390/cells8080914
Gobena, M., Elzo, M. A., and Mateescu, R. G. (2018). Population structure and genomic breed composition in an angus–brahman crossbred cattle population. Front. Genet. 9, 90. doi:10.3389/fgene.2018.00090
Goldammer, T., Brunner, R. M., and Schwerin, M. (1997). Comparative analysis of Y chromosome structure in Bos taurus and B. indicus by FISH using region-specific, microdissected, and locus-specific DNA probes. Cytogenet Cell Genet. 77, 238–241. doi:10.1159/000134584
Goszczynski, D. E., Corbi-Botto, C. M., Durand, H. M., Rogberg-Muñoz, A., Munilla, S., Peral-Garcia, P., et al. (2018). Evidence of positive selection towards Zebuine haplotypes in the BoLA region of Brangus cattle. Animal 12, 215–223. doi:10.1017/S1751731117001380
Hanif, Q., Farooq, M., Amin, I., Mansoor, S., Zhang, Y., and Khan, Q. M. (2018). In silico identification of conserved miRNAs and their selective target gene prediction in indicine (Bos indicus) cattle. PLoS One 13, e0206154. doi:10.1371/journal.pone.0206154
Hao, D., Wang, X., Yang, Y., Chen, H., Thomsen, B., and Holm, L.-E. (2023). MicroRNA sequence variation can impact interactions with target mRNA in cattle. Gene 868, 147373. doi:10.1016/j.gene.2023.147373
Hiendleder, S., Lewalski, H., and Janke, A. (2008). Complete mitochondrial genomes of Bos taurus and Bos indicus provide new insights into intra-species variation, taxonomy and domestication. Cytogenet. Genome Res. 120, 150–156. doi:10.1159/000118756
Jopling, C. (2012). Liver-specific microRNA-122: biogenesis and function. RNA Biol. 9, 137–142. doi:10.4161/rna.18827
Kalvari, I., Nawrocki, E. P., Argasinska, J., Quinones-Olvera, N., Finn, R. D., Bateman, A., et al. (2018). Non-coding RNA analysis using the Rfam database. Curr. Protoc. Bioinforma. 62, e51. doi:10.1002/cpbi.51
Kalvari, I., Nawrocki, E. P., Ontiveros-Palacios, N., Argasinska, J., Lamkiewicz, K., Marz, M., et al. (2020). Rfam 14: expanded coverage of metagenomic, viral and microRNA families. Nucleic Acids Res. 49, D192–D200. doi:10.1093/nar/gkaa1047
Kaneko, K., Xu, P., Cordonier, E. L., Chen, S. S., Ng, A., Xu, Y., et al. (2016). Neuronal Rap1 regulates energy balance, glucose homeostasis, and leptin actions. Cell Rep. 16, 3003–3015. doi:10.1016/j.celrep.2016.08.039
Kozomara, A., and Griffiths-Jones, S. (2011). miRBase: integrating microRNA annotation and deep-sequencing data. Nucleic acids Res. 39, D152–D157. doi:10.1093/nar/gkq1027
Langfelder, P., and Horvath, S. (2008). WGCNA: an R package for weighted correlation network analysis. BMC Bioinforma. 9, 559. doi:10.1186/1471-2105-9-559
Liao, Y., Smyth, G. K., and Shi, W. (2019). The R package Rsubread is easier, faster, cheaper and better for alignment and quantification of RNA sequencing reads. Nucleic Acids Res. 47, e47. doi:10.1093/nar/gkz114
Li, B., Dong, J., Yu, J., Fan, Y., Shang, L., Zhou, X., et al. (2020). Pinpointing miRNA and genes enrichment over trait-relevant tissue network in Genome-Wide Association Studies. BMC Med. Genomics 13, 191. doi:10.1186/s12920-020-00830-w
Liu, G., Liao, Y., Sun, B., Guo, Y., Deng, M., Li, Y., et al. (2020). Effects of chronic heat stress on mRNA and miRNA expressions in dairy cows. Gene 742, 144550. doi:10.1016/j.gene.2020.144550
Liu, R., Tearle, R., Low, W. Y., Chen, T., Thomsen, D., Smith, T. P. L., et al. (2021). Distinctive gene expression patterns and imprinting signatures revealed in reciprocal crosses between cattle sub-species. BMC Genomics 22, 410. doi:10.1186/s12864-021-07667-2
Li, Z., He, J., Yang, F., Yin, S., Gao, Z., Chen, W., et al. (2023). A look under the hood of genomic-estimated breed compositions for brangus cattle: what have we learned? Front. Genet. 14, 1080279. doi:10.3389/fgene.2023.1080279
Loftus, R. T., Machugh, D. E., Bradley, D. G., Sharp, P. M., and Cunningham, P. (1994). Evidence for two independent domestications of cattle. Proc. Natl. Acad. Sci. 91, 2757–2761. doi:10.1073/pnas.91.7.2757
Lotto, J., Stephan, T. L., and Hoodless, P. A. (2023). Fetal liver development and implications for liver disease pathogenesis. Nat. Rev. Gastroenterology Hepatology 20, 561–581. doi:10.1038/s41575-023-00775-2
Louveau, I., Perruchot, M. H., Bonnet, M., and Gondret, F. (2016). Invited review: pre- and postnatal adipose tissue development in farm animals: from stem cells to adipocyte physiology. Animal 10, 1839–1847. doi:10.1017/S1751731116000872
Love, M. I., Huber, W., and Anders, S. (2014). Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. doi:10.1186/s13059-014-0550-8
Low, W. Y., Tearle, R., Liu, R., Koren, S., Rhie, A., Bickhart, D. M., et al. (2020). Haplotype-resolved genomes provide insights into structural variation and gene content in Angus and Brahman cattle. Nat. Commun. 11, 1–14. doi:10.1038/s41467-020-15848-y
Macdougald, O. A., and Lane, M. D. (1995). Transcriptional regulation of gene expression during adipocyte differentiation. Annu. Rev. Biochem. 64, 345–373. doi:10.1146/annurev.bi.64.070195.002021
Macfarlane, L. A., and Murphy, P. R. (2010). MicroRNA: biogenesis, function and role in cancer. Curr. Genomics 11, 537–561. doi:10.2174/138920210793175895
Machugh, D. E., Larson, G., and Orlando, L. (2017). Taming the past: ancient DNA and the study of animal domestication. Annu. Rev. Animal Biosci. 5, 329–351. doi:10.1146/annurev-animal-022516-022747
Mao, W., Albrecht, E., Teuscher, F., Yang, Q., Zhao, R., and Wegner, J. (2008). Growth-and breed-related changes of fetal development in cattle. Asian-Australasian J. Animal Sci. 21, 640–647. doi:10.5713/ajas.2008.70293
Martin, M. (2011). Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.J. 17, 10. doi:10.14806/ej.17.1.200
Mctavish, E. J., Decker, J. E., Schnabel, R. D., Taylor, J. F., and Hillis, D. M. (2013). New World cattle show ancestry from multiple independent domestication events. Proc. Natl. Acad. Sci. 110, E1398–E1406. doi:10.1073/pnas.1303367110
Morales-Roselló, J., Loscalzo, G., García-Lopez, E. M., García-Gimenez, J. L., and Perales-Marín, A. (2022). MicroRNA-132 is overexpressed in fetuses with late-onset fetal growth restriction. Health Sci. Rep. 5, e558. doi:10.1002/hsr2.558
Mukiibi, R., Johnston, D., Vinsky, M., Fitzsimmons, C., Stothard, P., Waters, S. M., et al. (2020). Bovine hepatic miRNAome profiling and differential miRNA expression analyses between beef steers with divergent feed efficiency phenotypes. Sci. Rep. 10, 19309. doi:10.1038/s41598-020-73885-5
Murholm, M., Dixen, K., and Hansen, J. B. (2010). Ras signalling regulates differentiation and UCP1 expression in models of brown adipogenesis. Biochim. Biophys. Acta 1800, 619–627. doi:10.1016/j.bbagen.2010.03.008
Park, J., Kim, M., Sun, K., An, Y. A., Gu, X., and Scherer, P. E. (2017). VEGF-A-Expressing adipose tissue shows rapid beiging and enhanced survival after transplantation and confers IL-4-independent metabolic improvements. Diabetes 66, 1479–1490. doi:10.2337/db16-1081
Peng, S., Zhang, X., Yu, L., Xu, Y., Zhou, Y., Qian, S., et al. (2022). NF-κB regulates brown adipocyte function through suppression of ANT2. Acta Pharm. Sin. B 12, 1186–1197. doi:10.1016/j.apsb.2021.10.023
Peterson, S., Thompson, J., Ufkin, M., Sathyanarayana, P., Liaw, L., and Congdon, C. B. (2014). Common features of microRNA target prediction tools. Front. Genet. 5, 23. doi:10.3389/fgene.2014.00023
Picard, B., Lefaucheur, L., Berri, C., and Duclos, M. J. (2002). Muscle fibre ontogenesis in farm animal species. Reprod. Nutr. Dev. 42, 415–431. doi:10.1051/rnd:2002035
Romao, J. M., Jin, W., He, M., Mcallister, T., and Guan Le, L. (2014). MicroRNAs in bovine adipogenesis: genomic context, expression and function. BMC Genomics 15, 137. doi:10.1186/1471-2164-15-137
Rosen, B. D., Bickhart, D. M., Schnabel, R. D., Koren, S., Elsik, C. G., Tseng, E., et al. (2020). De novo assembly of the cattle reference genome with single-molecule sequencing. Gigascience 9, giaa021. doi:10.1093/gigascience/giaa021
Rottiers, V., and Näär, A. M. (2012). MicroRNAs in metabolism and metabolic disorders. Nat. Rev. Mol. Cell Biol. 13, 239–250. doi:10.1038/nrm3313
Rui, L. (2014). Energy metabolism in the liver. Compr. Physiol. 4, 177–197. doi:10.1002/cphy.c130024
Sadkowski, T., Ciecierska, A., Oprządek, J., and Balcerek, E. (2018). Breed-dependent microRNA expression in the primary culture of skeletal muscle cells subjected to myogenic differentiation. BMC Genomics 19, 109. doi:10.1186/s12864-018-4492-5
Schmieder, R., and Edwards, R. (2011). Quality control and preprocessing of metagenomic datasets. Bioinformatics 27, 863–864. doi:10.1093/bioinformatics/btr026
Shen, J., Hao, Z., Luo, Y., Zhen, H., Liu, Y., Wang, J., et al. (2022). Deep small RNA sequencing reveals important miRNAs related to muscle development and intramuscular fat deposition in longissimus dorsi muscle from different goat breeds. Front. Veterinary Sci. 9, 911166. doi:10.3389/fvets.2022.911166
Sun, J., Zhou, Y., Cai, H., Lan, X., Lei, C., Zhao, X., et al. (2014). Discovery of novel and differentially expressed MicroRNAs between fetal and adult backfat in cattle. PLOS ONE 9, e90244. doi:10.1371/journal.pone.0090244
Sun, Y., Wang, H., Wang, W., Lu, J., Zhang, J., Luo, X., et al. (2022). Glutamatergic and GABAergic neurons in the preoptic area of the hypothalamus play key roles in menopausal hot flashes. Front. Aging Neurosci. 14, 993955. doi:10.3389/fnagi.2022.993955
Thakral, S., and Ghoshal, K. (2015). miR-122 is a unique molecule with great potential in diagnosis, prognosis of liver disease, and therapy both as miRNA mimic and antimir. Curr. Gene Ther. 15, 142–150. doi:10.2174/1566523214666141224095610
Thorne, A. M., Ubbink, R., Brüggenwirth, I. M. A., Nijsten, M. W., Porte, R. J., and Meijer, V. E. D. (2020). Hyperthermia-induced changes in liver physiology and metabolism: a rationale for hyperthermic machine perfusion. Am. J. Physiology-Gastrointestinal Liver Physiology 319, G43–G50. doi:10.1152/ajpgi.00101.2020
Tiniakos, D. G., Lee, J. A., and Burt, A. D. (1996). Innervation of the liver: morphology and function. Liver 16, 151–160. doi:10.1111/j.1600-0676.1996.tb00721.x
Vrijens, K., Bollati, V., and Nawrot, T. S. (2015). MicroRNAs as potential signatures of environmental exposure or effect: a systematic review. Environ. Health Perspect. 123, 399–411. doi:10.1289/ehp.1408459
Xiang, R., Ghanipoor-Samami, M., Johns, W. H., Eindorf, T., Rutley, D. L., Kruk, Z. A., et al. (2013). Maternal and paternal genomes differentially affect myofibre characteristics and muscle weights of bovine fetuses at midgestation. PLOS ONE 8, e53402. doi:10.1371/journal.pone.0053402
Xiang, R., Lee, A. M., Eindorf, T., Javadmanesh, A., Ghanipoor-Samami, M., Gugger, M., et al. (2014). Widespread differential maternal and paternal genome effects on fetal bone phenotype at mid-gestation. J. Bone Mineral Res. 29, 2392–2404. doi:10.1002/jbmr.2263
Xie, S.-J., Li, J.-H., Chen, H.-F., Tan, Y.-Y., Liu, S.-R., Zhang, Y., et al. (2018). Inhibition of the JNK/MAPK signaling pathway by myogenesis-associated miRNAs is required for skeletal muscle development. Cell Death Differ. 25, 1581–1597. doi:10.1038/s41418-018-0063-1
Yan, X., Huang, Y., Zhao, J. X., Rogers, C. J., Zhu, M. J., Ford, S. P., et al. (2013). Maternal obesity downregulates microRNA let-7g expression, a possible mechanism for enhanced adipogenesis during ovine fetal skeletal muscle development. Int. J. Obes. (Lond) 37, 568–575. doi:10.1038/ijo.2012.69
Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. Omics 16, 284–287. doi:10.1089/omi.2011.0118
Zhang, W., and Liu, H. T. (2002). MAPK signal pathways in the regulation of cell proliferation in mammalian cells. Cell Res. 12, 9–18. doi:10.1038/sj.cr.7290105
Keywords: mRNA-miRNA, differentially expressed micrornas (DE-miRNAs), differentially expressed genes (DEGs), Brahman, Angus, fetal liver
Citation: MacPhillamy C, Ren Y, Chen T, Hiendleder S and Low WY (2023) MicroRNA breed and parent-of-origin effects provide insights into biological pathways differentiating cattle subspecies in fetal liver. Front. Genet. 14:1329939. doi: 10.3389/fgene.2023.1329939
Received: 30 October 2023; Accepted: 01 December 2023;
Published: 13 December 2023.
Edited by:
Penny Riggs, Texas A and M University, United StatesReviewed by:
Tara G. McDaneld, Agricultural Research Service (USDA), United StatesShudong Liu, Tarim University, China
Copyright © 2023 MacPhillamy, Ren, Chen, Hiendleder and Low. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Wai Yee Low, wai.low@adelaide.edu.au