Introduction

In domestic cattle studies, whole-genome sequencing has been successfully applied to investigate genetic diversity, genomic architecture, history of cattle populations, selection signatures, and economically and physiologically important genomic variations. These are all important investigations for characterizing cattle genetic resources for agriculture and food production1,2. However, thus far, genomic characterization at the functional level is lacking for native cattle breeds, which is crucial to understanding their ability to adapt to various environmental circumstances3. RNA sequencing may allow us to unravel the critical knowledge needed for understanding how environmental circumstances, demographic factors and breeding histories are reflected in gene expression profiles in different cattle breeds3. Moreover, transcriptome profiling can reveal important candidate genes for production, fertility and health traits in domestic animal species and breeds, and the genetic and evolutionary basis of these complex traits3,4,5.

Here, we used high-throughput RNA sequencing technology to investigate and compare gene expression profiles in four adipose tissues (metacarpal, perirenal, tailhead and prescapular adipose tissues), representing visceral, peripheral and bone marrow fat tissue, in three native breeds and one commercial cattle breed. These adipose tissues are important organs for many physiological functions for survival and successful reproduction6. Due to the energy storing and channelling properties of the adipocytes, adipose tissue plays a crucial role in energy metabolism, for example, cold-induced adaptive thermogenesis, as well as in cushioning internal organs and insulating the body7. Adipose tissues help control energy balance and metabolic activity and are also capable of restructuring on the basis of nutritional changes7. There are two main types of adipose tissues in mammals: brown and white. White adipose tissue stores energy, whereas brown adipose tissue (making up to 2% of the bodyweight of calves up to 26 months of age8) dissipates stored energy as heat by burning fatty acids to maintain body temperature9. Fat in bone marrow acts as an energy reservoir contributing to metabolic processes and undergoes dynamic changes, for example, as a result of starvation10. Weldenegodguad et al.9, recently reported that the gene expression profiles of metacarpal tissues were distinct from those of perirenal and prescapular fats in semi-domestic reindeer (Rangifer tarandus), providing interesting insights into the nature of adipose tissues, with candidate genes involved in immune response and energy metabolism.

The native cattle breeds included in this study – Mirandesa cattle from Portugal; Northern Finncattle from Finland; and Yakutian cattle from the Sakha Republic (Yakutia), the Russian Federation – have adapted to various biogeographical and production environments and have different genetic origins and breeding histories. In addition, we investigated the Holstein cattle breed, which is the high-producing and most popular dairy cattle breed worldwide (Fig. 1)6,7,9. Studying the genetic resources of native, locally adapted cattle is essential due to climate change; the current productive commercial breeds, such as the Holstein, may not express the genetic variations relevant for future breeding and adaptation to changing environments11. Yakutian cattle, for example, live in the area of Sakha (Yakutia), which has some of the most extreme low temperatures in the settled world (− 70 °C to + 30 °C); thus, this breed is of great interest when studying adaptation to cold12. Adaptation to cold temperatures simultaneously implies adaptation to a different type of diet consisting natural forage.

Fig. 1
figure 1

Graphical summary of gene expression in adipose tissues. a PCA plot based on the gene expression of the tissue samples representing Yakutian cattle (YAK), Mirandesa cattle (MIR), Holstein cattle (HOL) and Northern Finncattle (NFC) samples based on expression profiles, with dot colours indicating tissue and breed; b shared and unique genes expressed in the tissues: a venn diagram of expressed genes with abundance > 0.1 transcripts per million reads (TPM) in metacarpal adipose tissue (MAT), tailhead adipose tissue (TAT), perirenal adipose tissue (PAT) and prescapular adipose tissue (PSAT).

Full size image

The Northern Finncattle breed is well adapted to harsh northern environments and is used for milk and meat production2,13. Challenging environmental conditions and limited natural resources have had an impact on the traditional diet of Northern Finncattle. Their feed typically included cooked soft stewed mash made from ingredients such as lichen and other less nutritious forage14 while currently the modern feeding system is practiced. Mirandesa is a native cattle breed from the north of Portugal, a region characterized by its moderate-to-warm Mediterranean climate with Continental influence. Previous studies highlighted the genetic differentiation of Mirandesa cattle, confirming a bottleneck and the consequent loss of genetic diversity15. This breed is known for its excellent meat quality, tenderness, longevity, fertility, temperament and adaptation to the rigorous winter and summer environment.

To get insights into the role of adipose tissues in shaping the adaptive capacity of these breeds, we performed an RNA-Seq study on four adipose tissues from three native, and one commercial cattle breed. Here we aim to explore the transcriptome profiles of these tissues and identify differentially expressed genes and their association with key physiological and metabolic activities in cattle breeds included in this study. Moreover, we are interested in identifying differences in gene expression profiles between intensively selected Holstein cattle and locally adapted native breeds. In this study, we have looked for the effect of gender and castration (in Yakutian male samples) in gene expression. However, it should be noted that the results could also be influenced by age, diet and other environmental factors.

Results

RNA-Seq and mapping

A total of 2,830,100,766 read pairs were generated from data of 81 samples of perirenal (n = 26), metacarpal (n = 26), tailhead (n = 26) and prescapular (n = 3) adipose tissues. After sequencing, a total of 1,156 multi-lane Fastq files were merged into 162 forward and reverse files representing 81 samples. The Phred quality scores from the reads of all the samples were > 30 (30.57–40.11). The number of reads per sample ranged between 27.1 million and 63.2 million, with an average of 35.2 million (See Additional file 1, Table S1 and S2).

Gene expression overview

A total of 20,714 genes (of 27,607 total listed bovine genes) were expressed in 81 samples. According to the RNA-Seq Expectation Maximization (RSEM) results, the highest number of genes with an abundance > 0.1 transcripts per million reads (TPM) were expressed in metacarpal adipose tissue (n = 17,401), followed by tailhead adipose tissue (n = 16,992), perirenal adipose tissue (n = 16,802), and prescapular adipose tissue (n = 16,060). Principal component analysis (PCA) of the normalized expression profiles revealed that PC1 and PC2 together explained ~ 50% of the variance. PC1 alone explained 44% of the variance and clearly separated metacarpal adipose tissue from the other types of adipose tissue (Fig. 1a). There were far fewer obvious divisions between the rest of the tissues, with tailhead and perirenal adipose tissues overlapping remarkably independently of breed origin.

As shown in Fig. 1b, 15,431 genes were commonly expressed in all four tissues. The metacarpal adipose tissue (n = 672; See Additional file 2, Table S3) harbored the highest number of uniquely expressed genes, followed by the tailhead (n = 215; See Additional file 2, Table S4), perirenal adipose tissue (n = 133; See Additional file 2, Table S5), and prescapular adipose tissue (n = 125; See Additional file 2, Table S6).

In metacarpal adipose tissue, the HOX family of genes (HOXD13, HOXD12), olfactory receptor genes (e.g., OR10AG82P, OR2AD1, OR2AZ1, OR2J1), and solute carrier family genes (e.g., SLC12A3, SLC25A2, SLC26A9) were uniquely expressed. The top GO terms associated with uniquely expressed genes in metacarpal adipose tissue (See Additional file 3, Table S7) included “GO:0048513 animal organ development”, “GO:0048522 positive regulation of cellular process”, “GO:0051716 cellular response to stimulus”, “GO:0048518 positive regulation of biological process”, “GO:0050896 response to stimulus”, “GO:0016043 cellular component organization”, “GO:0071840 cellular component organization or biogenesis”, “GO:0048856 anatomical structure development”, as well as metabolic process, organic substance metabolic process, cellular metabolic process, primary metabolic process, and macromolecule metabolic process, and “GO:0000003 reproduction”. Moreover, two KEGG pathways (metabolic pathways and neuroactive ligand‒receptor interaction) were associated with uniquely expressed genes in metacarpal adipose tissue (See Additional file 3, Table S8).

The most abundant gene in the tailhead adipose tissue (See Additional file 2, Table S4) was bta-mir-192, which is a microRNA possibly associated with genes targeting lipogenesis and the regulation of adipose deposition and differentiation16. Genes such as HOXC12, HOXC13, SERPINA3-1 and SERPINA3-3 were uniquely expressed in tailhead adipose tissue. Similarly, among the uniquely expressed genes in perirenal adipose tissue (See Additional file 2, Table S5), a set of microRNAs were found, including mir-197 and TEX11. Finally, the uniquely expressed genes in prescapular adipose tissue included PHGR1, KRT7 and DRGX (See Additional file 2, Table S6). We did not find any GO terms or KEGG pathways for the lists of the highest abundance genes in tailhead, perirenal or prescapular adipose tissues.

Castrated vs. uncastrated male samples

We obtained both castrated and noncastrated male samples from Yakutian cattle. Previous adipose tissue studies have shown gene expression differences in the adipose tissues of castrated and noncastrated mammals; for example, in castrated mice, brown adipose tissue will begin to convert to white adipose tissue17. It was therefore important to compare the gene expression profiles of adipose tissue samples between the castrated and noncastrated male samples of Yakutian cattle. However, only a few DEGs were detected between Yakutian cattle castrated and noncastrated male samples. Only one DEG (TDH) was detected in the metacarpal adipose tissue (See Additional file 4, Table S9), and this gene was upregulated in the noncastrated samples. Similarly, four DEGs were present in the tailhead (See Additional file 4, Table S9) adipose tissue, of which two (C11H2orf50 and ENSBTAG00000007075) were upregulated in the castrated animals, and two (SARDH and GSTA2) were upregulated in the noncastrated animals. The expression of all three DEGs (ENSBTAG00000025258, CXCL9, and ENSBTAG00000052522) in perirenal adipose tissue was upregulated in the castrated animals (See Additional file 4, Table S9). Owing to the minimal effect of castration in the present Yakutian cattle samples, all the samples were grouped together in further analyses, taking these DEGs into consideration in our conclusions (See Additional file 4, Table S9).

Gender differences

We further performed an analysis to test for the gender effect within breeds. Our comparison of the gene expression profiles between male and female individuals of the native cattle breeds revealed a total of 345 significant DEGs in the three adipose tissues (Table 1). The highest number of DEGs between the sexes was identified in the tailhead adipose tissue.

Table 1 Number of upregulated genes in males and females.
Full size table

The gender-wise differential gene expression comparisons of metacarpal adipose tissue identified 47, 20 and 26 DEGs in Northern Finncattle, Mirandesa cattle, and Yakutian cattle, respectively (Additional file 5, Tables S10 to S12). Gene expression analysis of perirenal adipose tissue revealed 12, 14, and 47 DEGs between males and females in Yakutian cattle, Mirandesa cattle, and Northern Finncattle respectively (See Additional file 6, Tables S13 to S15). Similar analysis comparing the sexes in tailhead adipose tissue yielded 13, 68, and 98 DEGs in Yakutian cattle, Mirandesa cattle, and Northern Finncattle respectively (See Additional file 7, Tables S16 to S18).

Differential gene expression profiles in the native cattle breeds

Given the fact that we have identified significant within-breed differences deriving from gender, the analysis of DEGs between the three native breeds was conducted correcting for gender effects. The highest number of DEGs (n = 26) was found in tailhead adipose tissue between Mirandesa cattle and Northern Finncattle (Table 2).

Table 2 Number of differentially expressed genes (DEGs) between the native breeds.
Full size table

In metacarpal adipose tissue, the comparison of Yakutian and Mirandesa cattle yielded 17 DEGs, with 7 upregulated for Yakutian and 10 upregulated for Mirandesa (See Additional file 8, Table S19). The gene with the greatest upregulation in Yakutian cattle compared to Mirandesa cattle was NR4A3, which is associated with fat deposition and carbohydrate metabolism18,19. The most highly upregulated gene in Mirandesa was PPP1R14C, which is associated with immune tolerance and longevity in cattle20,21. A comparison of Yakutian cattle and Northern Finncattle yielded 7 genes (See Additional file 8, Table S20), with one (EML6) gene upregulated in Yakutan cattle and 6 genes upregulated in Northern Finncattle. In addition, EML6 is associated with reproductive traits and body composition traits in cattle22,23. Finally, the comparison of Mirandesa cattle and Northern Finncattle cattle yielded only one gene (See Additional file 8, Table S21).

In tailhead adipose tissue, the comparison of Yakutian and Mirandesa cattle yielded four genes in total; one (JPH4) was upregulated in Yakutian, and three (PTHLH, ENSBTAG00000048635, and PPP1R14C) were upregulated in Mirandesa (Additional file 9, Table S22). For PPP1R14C, we observed a pattern of expression similar to that in metacarpal adipose tissue. A comparison of Yakutian cattle and Northern Finncattle yielded two genes, with IGFBP2 upregulated in Yakutian, and SLC2A1 upregulated in Northern Finncattle (See Additional file 9, Table S23). A comparison of tailhead adipose tissue from Mirandesa and Northern Finncattle revealed 26 genes, 19 of which were upregulated in Mirandesa cattle and seven were upregulated in Northern Finncattle (Additional file 9, Table S24). The most upregulated gene in Mirandesa cattle in this comparison was IL6.

In perirenal adipose tissue, a total of nine genes were significantly differentially expressed between Yakutian and Mirandesa cattle, 8 of which were upregulated in Mirandesa cattle (Additional file 10, Table S25). PLEKHG7 was the only gene upregulated in Yakutian cattle. A comparison of Yakutian and Northern Finncattle revealed 17 DEGs, with two genes upregulated in Yakutian cattle, and 15 genes upregulated in Northern Finncattle (Additional file 10, Table S26). The most upregulated gene in Yakutian cattle was PRSS42. A comparison of the Mirandesa and Northern Finncattle genomes yielded 4 genes (Additional file 10, Table S27). Only one gene, MT1E, was upregulated in the Mirandesa cattle and the rest were upregulated in Northern Finncattle.

Differences in gene expression between the commercial and native breeds

The Holstein breed is the most popular dairy cattle worldwide and is selected for milk production traits24. Here, we examined whether there are differences in gene expression between commercial breeds and less intensively selected local native breeds. In general, we identified a higher number of DEGs in pairwise comparisons of native breeds vs Holstein cattle (Table 3). Among the native breeds the number of DEGs is lower (Table 2).

Table 3 Number of differentially expressed genes (DEGs) between female samples of commercial and native breeds.
Full size table

Holstein cattle had the highest number of DEGs when compared with the aboriginal Yakutian cattle in all three adipose tissues (Table 3). In the tailhead adipose tissue, 580 DEGs were identified in the Yakutian cattle vs. Holstein comparison, while in the perirenal adipose tissue, 519 DEGs were found, and in the metacarpal tissue, 78 DEGs were found (See Additional file 11, Table S28 to S36). The expression of genes involved in the immune response and resistance/susceptibility to disease (PRF1, CD8B, IDO1, CCL5, CD3D, CD3E, and GNLY), in milk composition traits (FCGR2B) and in feed intake (GIMAP4) was strongly upregulated in Yakutian cattle (log2-fold change LFC > 4). In Mirandesa, strongly upregulated genes were found in the following categories: fertility (STC1, TCIM), milk production/composition (MAPK15), meat quality (CLK1), feed efficiency (MT1E, CCL8, HI-4), and disease response (HP, IL6). Interestingly, the lowest number of DEGs in all analyses were found between the native Northern Finncattle and Holstein cattle.

The enrichment analysis of GO terms associated with significant DEGs revealed a total of 126 GO terms with a p-value < 0.05 only for the three native breeds vs Holstein comparisons: Yakutian cattle in the tailhead adipose tissue, Yakutian cattle in the perirenal adipose tissue and Mirandesa cattle in the perirenal adipose tissue (See Additional file 11, Table S37 to S42). The DEGs between Yakutian cattle and Holstein cattle in the tailhead adipose tissue had the highest number (n = 86) of associated GO terms. There were 59 GO terms associated with upregulated genes in Yakutian cattle among the top (p-value < 0.05), such as “GO:0031328 positive regulation of biosynthetic process” and “GO:0009891 positive regulation of cellular biosynthetic process”. Other interesting GO terms in this list included “GO:0006955 immune response”, “GO:0006952 defense response”, “GO:0098542 defense response to other organism”, and “GO:0002682 regulation of immune system process”, which indicate differences in disease resistance and general immunity between Holstein cattle and Yakutian cattle. There were 27 GO terms associated with upregulated genes in Holstein cattle compared with Yakutian cattle, the top of which included “GO:1901564 organonitrogen compound metabolic process”, “GO:0008610 lipid biosynthetic process”, “GO:0071702 organic substance transport”, “GO:0046395 carboxylic acid catabolic process”, and “GO:0071704 organic substance metabolic process”. Other interesting GO terms associated with upregulated genes in this category included “GO:0006629 lipid metabolic process”, “GO:0006631 fatty acid metabolic process”, and “GO:0044255 cellular lipid metabolic process”, which suggest a possible difference in lipogenesis between these two breeds.

The KEGG pathways associated with the genes upregulated in the Yakutian cattle vs Holstein cattle in the tailhead adipose tissue comparison included “Th1 and Th2 cell differentiation”, “Th17 cell differentiation”, “antigen processing and presentation”, “chemokine signaling pathway”, and “cytokine‒cytokine receptor interaction” (Additional file 12, Table S43). On the other hand, the KEGG pathways associated with downregulated genes in Holstein cattle in the tailhead adipose tissue were “metabolic pathways”, “carbon metabolism”, “valine, leucine and isoleucine degradation”, “oxidative phosphorylation” and “thermogenesis” (Additional file 12, Table S44).

A total of 30 GO terms were generated for the DEGs between Yakutian cattle and Holstein cattle in the perirenal adipose tissue. There were 5 GO terms associated with upregulated genes in Yakutian cattle, “GO:0002376 immune system process” and “GO:0006955 immune response”, again indicating possible differences in disease resistance between these two breeds. The downregulated genes in Holstein cattle were associated with 25 GO terms, such as “GO:0044281 small molecule metabolic process”, “GO:0044238 primary metabolic process”, “GO:0044237 cellular metabolic process” and “GO:0008152 metabolic process”, as well as “GO:0006629 lipid metabolic process” and “GO:0008610 lipid biosynthetic process”.

Pathway analysis revealed 18 KEGG pathways associated with upregulated genes in Yakutian cattle, for example, “viral protein interaction”, “metabolic pathways” and “thermogenesis” (Additional file 12, Table S45), while no KEGG pathways were found for downregulated genes in the perirenal tissue of Holstein cattle. The significant DEGs between Mirandesa cattle and Holstein cattle in the perirenal adipose tissue generated 10 GO terms in total. There was one upregulated GO term, “immune system process”, and 9 downregulated GO terms, including “lipid metabolic processes”. Overall, top GO terms and KEGG pathways (Tables 4 and 5) were generally unique but had some similarities between tissues and breeds.

Table 4 GO terms of interest based on top differentially expressed genes from comparisons of native cattle breeds versus commercial Holstein.
Full size table
Table 5 Top KEGG pathways associated with DEGs in two adipose tissues (Perirenal and Tailhead) between Yakutian (YAK) and Holstein cattle.
Full size table

YAK Yakutian cattle.

Discussion

Transcriptome studies have shown that gene expression patterns in adipose tissue differ between breeds, sexes, and adipose tissue depots and, more recently, that there are differences in adaptation to extreme environments9,25. Here, we studied the transcriptome profiles of four adipose tissues from three native cattle breeds and one commercial breed, all of which originated from different geographical locations and climates. In the present study, we identified a total of 16,060–17,401 genes in the analysed tissues, which covered approximately 63.8% of the list of genes available for the Bos taurus reference genome (ARS-UCD1.2.). We found that metacarpal adipose tissue displayed a distinct pattern of gene expression compared to the other three tissues, which is in line with the functional role of this tissue, as highlighted in a recent reindeer study9. Similarly, metacarpal adipose tissue had the highest number of uniquely expressed genes (n = 671), and functional annotation of those unique genes indicated the role of this tissue in the development and regulation of cellular processes, reproduction, as well as metabolism. Several genes belonging to the homeobox (HOX) family were uniquely expressed in metacarpal adipose tissue. In particular, HOXD13 had the highest abundance among all uniquely expressed genes in the metacarpal adipose tissue (Additional File 3: S1). HOXD13 is a highly conserved gene belonging to the HOX family of genes and is responsible for morphogenesis, limb development and genital development26,27. In murines, the HOX gene family has been associated with cell differentiation toward adipogenesis28, and in cattle, the HOX gene family is associated with muscularity traits29 and limb development30. The distinct gene expression profiles of metacarpal adipose tissue and the presence of homeobox genes among the uniquely expressed genes in metacarpal adipose tissue are in agreement with the findings of a recent study in reindeer (Rangifer tarandus)9. On the other hand, in tailhead tissue, the uniquely expressed members of the HOX family of genes are HOXC12 and HOXC13, which have been previously associated with thermotolerance in African cattle31. Similarly, SERPINA3-1, one of the uniquely expressed genes in the tailhead adipose tissue may be associated with growth traits and development in Chinese cattle32.

Here we investigated the gene expression profiles of physiologically vital adipose tissues in native cattle breeds adapted to North European boreal (Northern Finncattle), South European Mediterranean (Mirandesa) and North Asian extreme continental (Yakutian cattle) environments. Moreover, the gene expression profiles of these native breeds were compared with those of the commercial Holstein breed that have been intensively selected. In our Yakutian cattle sample set, we also had RNA-Seq data from castrated males, but castration appeared to have a minimal effect on the gene expression profiles; thus, the data from the castrated and uncastrated males were pooled in the subsequent analyses. However, within the native breeds, we identified several DEGs between the female and male samples; therefore, sex was considered as a second factor in the DEG analysis between the three native cattle breeds. It is important to note that there may have been age, diet and other factors related confounding effects on our results. Furthermore, the statistical power of some comparisons is low due to smaller sample sizes. The breedwise comparisons showed that Yakutian cattle exhibited superior disease resistance traits compared to those of Mirandesa and Northern Finncattle, thus promoting adaptation to a challenging environment. Several of the genes downregulated in Yakutian cattle were related to the susceptibility of a host to diseases, including Mycobacterium avium paratuberculosis, bovine tuberculosis, mastitis, brucellosis, and bovine respiratory disease (CD209, EEF1A2, SLC2A1). The CD209 gene plays a key role in the pathogenesis of bovine paratuberculosis; furthermore, there is evidence for a genetic basis for susceptibility to Mycobacterium avium subspecies of bovine paratuberculosis in cattle33. This gene was downregulated in the metacarpal adipose tissue of Yakutian cattle in comparison to that of Mirandesa cattle. ACKR1 has been shown to be related to mastitis, as it was found to be differentially expressed in comparison with healthy udder tissue in other studies34. This is interesting because local Yakutian cattle farmers ‘help’ their cows to stay resistant by exposing the udder to cold but also when it gets too cold, protecting the udders via thermal insulation. This also involves traditional culture-based observations as to which cows need such additional shelter and which are resilient enough to withstand this. We could hypothesize that this practice may have contributed to artificial selection to increase mastitis resistance. SLC2A1 is associated with the immunological response to bovine respiratory disease in Washington and Colorado beef cattle35. On the other hand, Yakutian cattle presented metabolism-related genes (e.g., TPRG1 and GLP1R).

Within-breed sex comparisons of metacarpal adipose tissues revealed that Yakutian cattle females exhibit a gene expression profile associated with higher metabolic efficiency in comparison with males. According to sex comparisons of all tissues of Yakutian cattle, TPRG1 was highly upregulated in females whereas GDF5 was the highest upregulated gene in males. TPRG1 is associated with feed efficiency traits, specifically high weight gain under low-feed conditions in cattle36 and GDF5 is associated with body traits in Bos taurus and Bos indicus37. In Mirandesa cattle females, upregulated genes were associated with latent tuberculosis (CD209)33, ovarian morphology and milk traits (ADCY5), mastitis and mastitis immunity (OSMR, PTX3)30, and temperament (BARHL2)38. In Mirandesa cattle males, the top upregulated genes were linked to immune function (SLC7A8)39 as well as susceptibility to BSE (PRND) and BVD (DHCR24)40,41. Similarly, the most upregulated genes in Northern Finncattle females and males were GNAO1 and CACNA1G, respectively. GNAO1 is known to affect lactation traits42, and CACNA1G is associated with feed efficiency43. Interestingly, in perirenal adipose tissue, the gene with the greatest upregulation in Yakutian cattle females was still TPRG1, highlighting the importance of this gene and the general trend toward feed efficiency in Yakutian cattle. The most highly upregulated gene in Yakutian cattle males was IL20RA, which is associated with susceptibility to bovine tuberculosis in cattle44. In Mirandesa cattle females, the highest known upregulated gene was MT1E, which is associated with postpartum oxidative stress45. Moreover, this gene was also present in the breed comparisons. This was followed by ACKR1, MT2A, and S100A12, which are associated with mastitis and mastitis response34,46,47. The highest known upregulated gene in Northern Finncattle females was JCHAIN, associated with mastitis response48. This was followed by UCHL1 and DYNC1I1, which are associated with splayed forelimbs at birth and limb development as well as body conformation traits30,49,50. In tailhead adipose tissue, the most upregulated gene in Mirandesa cattle females was MT1E. Similarly, the most upregulated DEG in Yakutian cattle females was GLP1R, which correlates with the trend in Yakutian cattle females. In addition to being associated with feed efficiency, this gene is also associated with obesity in humans51, and GLP-1 analogues such as semaglutide are used for the treatment of obesity and diabetes52. These cattle do not eat industrially produced food. During summer, they feed on fresh grass, but in spring, they also eat some green plants from previous years on the shores of rivers. Moreover, in the late summer/autumn, they also travel to the forest and feed on shrubs, leaves and twigs. According to herders, this makes the animals more resistant to different kinds of feed. In the most important ‘bottleneck’ season for feeding, in spring, the feed that herders give to their animals is most different by sex: lactating females obtain the freshest hay. Herders try to ensure that they always receive green hay, while males obtain different qualities of hay. For example, hay from later harvests or of lower quality that was already brown. This practice, according to herders, is a regular measure in spring when good quality fodder is scarce, so they prioritize females that feed calves and people with milk. During fieldwork when one sample animal for this study was slaughtered in the village of Kustur, the local veterinarian measured the length of the intestines of the slaughtered animal and commented that it was much longer than that of other breeds of cattle in other areas. According to him, this is a result of the long-term adaptation of the animals to a more diverse diet and the need for the animals to digest “tougher” fodder such as twigs, branches and shrubs in the forest.

Breedwise comparisons showed that genes associated with lactation53, milk production54, immunity and disease55, were upregulated in Northern Finncattle. Similarly, the expression of genes related to lactation and other mammary traits (SLC38A2, SLC35B1), fertility (SLC38A2, MT1E), susceptibility to paratuberculosis, mastitis and trypanosomiasis (IL6), and feed intake/efficiency (MT1E, SLC38A2) was upregulated in Mirandesa cattle45,53,54,55. This is somewhat related to the characteristics of this breed, which is well adapted to a rough forage and the local environment, but also has an incidence of paratuberculosis. In addition, this breed is known for its ease of calving and longevity. SLC35B1 is associated with glucose transportation and the lactose biosynthesis pathway. This gene was found to be expressed in lactose synthesis56. SLC38A2 has been found to be associated with postmortem carcass traits in Nellore cattle, specifically with the ribeye area and the amount of meat in the carcass57. This gene has also been associated with pregnancy maintenance in cattle58 and peak lactation in sows59. This gene specifically plays an important role in amino acid transport during lactation in sows and may influence dairy quality. Interestingly, MT1E was found to be upregulated in Mirandesa cattle in all comparisons with other native breeds in perirenal adipose tissue. This protein-coding gene is thought to enable zinc ion binding activity; zinc ions have a limited ability to bind to metallothionein, which is sensitive to oxidative stress. The oxidative stress results in elevated concentrations of free zinc and induces a pro-oxidative state. Therefore, this gene is associated with postpartum oxidative stress45 in cattle and was also found to be associated with high residual feed intake in heifers adjusted for backfat thickness60. The same gene was found to be highly upregulated in all adipose tissues of Mirandesa cattle females compared with males, possibly indicating that this feature is superior in females.

Based on the comparisons with Holstein cattle, it can be concluded that there are significant differences in gene expression between Holstein and Yakutian cattle in adipose tissues. The most differentially expressed genes were observed in tailhead adipose tissue from Holstein cattle compared to Yakutian cattle, followed by perirenal adipose tissue. Furthermore, GO term analysis of the DEGs revealed multiple upregulated immune-related terms in both tailhead and perirenal adipose tissues, indicating potential differences in disease resistance and immunity between the two breeds. Interestingly, the MT1E gene was again highly upregulated in Mirandesa cattle compared to Holstein cattle. The frequency of upregulation of this gene in Mirandesa cattle points to its importance in this breed’s traits. Highly upregulated genes in Northern Finncattle were associated with meat quality (PLTP, ACTN2, and ABCD2)61,62, immune response (CDHR1)63, and milk composition traits as well as thermotolerance in other species (DUSP1)64,65. These findings point to the cold adaptive qualities of this native breed, which has been reared for centuries in northern Finland and as far north as the Lapland.

Additionally, KEGG pathway analysis of tailhead adipose tissues from Holstein cattle and Yakutian cattle suggested that there may be differences in energy metabolism and immune system functions between the two breeds. In perirenal adipose tissue, relatively fewer genes were differentially expressed between Holstein and Yakutian cattle, but these differences were still significant. GO term analysis revealed upregulation of immune system processes and immune responses in Holstein cattle compared to Yakutian cattle, which is consistent with the results observed in tailhead adipose tissue.

Overall, the results suggest that there are significant differences in gene expression and biological pathways between Holstein and Yakutian cattle, particularly in terms of immune-related functions and energy metabolism. Further studies on the functional roles of these genes and pathways could provide insights into the underlying mechanisms of these differences and their potential implications for production and disease resistance in cattle.

Conclusions

The novelty of this study lies in its comprehensive investigation of differential gene expression in multiple cattle breed adipose gene transcriptomes. By examining different adipose tissue types (metacarpal, perirenal, and tailhead) across various cattle breeds (Yakutian, Mirandesa, Northern Finncattle, and Holstein), this study provides unique insights into the genetic adaptations that have evolved in response to diverse environmental conditions and selective pressures. There are likely differences in the functions and adaptations of these adipose tissues among the different breeds of cattle. For example, the genes upregulated in Yakutian cattle in metacarpal adipose tissue and perirenal adipose tissue, such as TPRG1 and IL20RA, suggest adaptations related to feed efficiency and susceptibility to tuberculosis, respectively. In contrast, the upregulated genes in perirenal adipose tissue of Mirandesa cattle, such as CD209 and MT1E, suggest adaptations related to tuberculosis susceptibility and postpartum oxidative stress, respectively. Furthermore, the study highlights the differences in adipose gene expression between males and females within each breed, which may provide insights into the underlying mechanisms of these differences and their potential implications for production, disease resistance, reproductive traits and metabolism. Overall, this research provides valuable information about the physiological and metabolic differences between cattle breeds, which could lead to the identification of genetic markers associated with specific traits. However, further research is needed to fully understand the functions and mechanisms of the genes and pathways identified in the different adipose tissues. By identifying breed-specific differences in gene expression related to immune functions, energy metabolism, and adaptation to local environments, this research highlights the potential for using genetic markers to improve livestock management, breeding strategies, and disease resistance.

Methods

Sample collection

A total of 81 adipose tissue samples were collected from 12 adult cows and 14 bulls (2–7 years old) of native breeds from different geographical locations and climates, namely, Northern Finncattle from Finland, Mirandesa cattle from Portugal, Yakutian cattle from Sakha in the Russian Federation, and Holstein cattle from Finland (Fig. 2; Table 6). The samples were randomly collected at slaughter from northern and central Finland, northern Portugal and northern and central Sakha (the Eveno-Bytantay and Magan regions, Sakha, Yakutia, the Russian Federation) during October 2015 and March 2016. Three different types of adipose tissue were collected: perirenal adipose tissue around the kidneys, tailhead adipose tissue between the tailhead and the tuber ischii, and metacarpal adipose tissue from the bone marrow in the diaphysis of the metacarpal bone (left front leg). In addition, three prescapular adipose tissue samples were collected from Yakutian cattle beneath the cervical muscle in front of the scapula. The samples were stored in RNAlater Solution (Ambion/QIAGEN, Valencia, CA, USA) after collection. Among the Yakutian samples, 3 were derived from castrated males and 3 from uncastrated males, whereas all the males from Finland (n = 5, Northern Finncattle) and Portugal (n = 3, Mirandesa) were uncastrated. All protocols and sample collections were performed in accordance with the legislations approved by the Russian authorization board (FS/UVN 03/163,733/07.04.2016) and the Animal Experiment Board in Finland (ESAVI/7034/04.10.07.2015).

Fig. 2
figure 2

The breeds analysed in this study. From top to bottom: Holstein cattle, Yakutian cattle, Mirandesa cattle and Northern Finncattle). Photos: Finnish Animal Breeding Association (Holstein cattle), Juha Kantanen (Yakutian cattle, Mirandesa cattle, Northern Finncattle).

Full size image
Table 6 Sample summary.
Full size table

Metacarpal adipose tissue = MAT, perirenal adipose tissue = PAT, tailhead adipose tissue = TAT and prescapular adipose tissue = PSAT. A total of 81 samples representing female (F) and male (M) individuals were sequenced. For instance, 3 + 5 refers to 3 females and 5 males. Out of the 6 Yakutian males, denoted by an asterisk “*”, 3 were castrated and 3 were noncastrated.

RNA extraction and sequencing

RNA extraction, library preparation, and sequencing were performed at The Finnish Functional Genomic Center (FFGC), Turku, Finland. Total RNA was extracted from adipose tissues (ca 30 mg/sample) using a Qiagen AllPrep DNA/RNA/miRNA kit. RNA extractions were performed according to the manufacturer’s protocol. The quality of the extracted RNA was confirmed with an Agilent Bioanalyzer 2100 (Agilent Technologies, Waldbronn, Germany), and the concentration of each sample was measured with a Nanodrop ND-2000 (Thermo Scientific, Wilmington, USA) and a Qbit Fluorometric Quantification Kit (Life Technologies). All the samples had an RNA integrity number (RIN) above 7.5.

Library preparation was performed according to the Illumina TruSeq Stranded mRNA Sample Preparation Guide (part #15031047). Unique Illumina TruSeq indexing adapters were ligated to each sample to pool several samples later in one flow cell lane. Library quality was inferred with an Advanced Analytical Fragment Analyser, and library concentration was inferred with a Qubit fluorometer; only good-quality libraries were sequenced.

The samples were normalized and pooled for automated cluster preparation, which was carried out with the Illumina cBot station. The libraries were analysed on an Illumina HiSeq 3000 platform. Paired-end sequencing with a 2 × 75 bp read length was performed with an 8 + 8 bp dual index run. Base calling and adapter trimming were performed using Illumina’s standard bcl2fastq2 software.

Sequence data analysis

The overall quality of the raw RNA-seq reads in fastq and aligned reads in BAM format were examined using FastQC v0.11.766. The FastQC reports were summarized together with the results from other analyses using MultiQC v.1.767. We used Spliced Transcripts Alignment to a Reference (STAR) (version 2.7.8a)68 for sequence alignment. Prior to alignment, a StarGenome was created using star v 2.7.8a with the parameters –runThreadN 4, –runMode genomeGenerate, and –sjdbOverhang 74 using the cattle reference genome (ARS-UCD1.2) and transcriptome (ARS-UCD1.2.103.gtf), which were downloaded from Ensembl. The same version of STAR was also used for alignments with the following parameters: –readFilesCommand zcat, –runThreadN 4, –outSAMtype BAM SortedByCoordinate, –quantMode GeneCounts, –twopassMode Basic. Quality control checks of the resulting BAM files resulting from STAR alignment were performed with RSeQC v4.0069. Gene body coverage showed uniform coverage, and junction saturation ranged from 107,569 to 323,347 at 100% of the reads. RSEM v 1.3.370 quantification was performed to establish gene expression levels for RNA-seq data with the parameter –paired-end.

We used DESeq2 1.32.071 to perform differential gene expression analysis for tissues, gender and breeds. We implemented differential gene expression analysis between breeds by incorporating gender as a second factor because we have clearly noted that gender does have an effect on gene expression. There were no male Holstein cattle; thus, separate gene expression analysis was performed between commercial Holstein females and females of the native breeds. Lowly expressed transcripts (< = 5 reads) were filtered out prior to running the DESeq command. We performed principal component analysis (PCA) using PlotPCA to obtain an overview of the sample distribution. Differentially expressed genes (DEGs) were identified by using several pairwise comparisons within and between populations. We used Log2Foldchange (log2FoldChange > 1.49) and adjusted P-value (padj < 0.05) to select the list of DEGs. The additional gene information, including the gene name and gene description, was retrieved for all DEGs using the biomaRt 2.48.3 bioconductor package72. Finally, we performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses73,74,75 on all tissue-specific uniquely expressed genes and later DEGs using GAGE v3.1676 for GO and KEGG analyses using biological pathway datasets.