Skip to main content

Fecal microbiota and their association with heat stress in Bos taurus



Humans have been influencing climate changes by burning fossil fuels, farming livestock, and cutting down rainforests, which has led to global temperature rise. This problem of global warming affects animals by causing heat stress, which negatively affects their health, biological functions, and reproduction. On the molecular level, it has been proved that heat stress changes the expression level of genes and therefore causes changes in proteome and metabolome. The importance of a microbiome in many studies showed that it is considered as individuals’ “second genome”. Physiological changes caused by heat stress may impact the microbiome composition.


In this study, we identified fecal microbiota associated with heat stress that was quantified by three metrics – rectal temperature, drooling, and respiratory scores represented by their Estimated Breeding Values. We analyzed the microbiota from 136 fecal samples of Chinese Holstein cows through a 16S rRNA gene sequencing approach. Statistical modeling was performed using a negative binomial regression. The analysis revealed the total number of 24 genera and 12 phyla associated with heat stress metrics. Rhizobium and Pseudobutyrivibrio turned out to be the most significant genera, while Acidobacteria and Gemmatimonadetes were the most significant phyla. Phylogenetic analysis revealed that three heat stress indicators quantify different metabolic ways of animals’ reaction to heat stress. Other studies already identified that those genera had significantly increased abundance in mice exposed to stressor-induced changes.


This study provides insights into the analysis of microbiome composition in cattle using heat stress measured as a continuous variable. The bacteria highly associated with heat stress were highlighted and can be used as biomarkers in further microbiological studies.

Peer Review reports


Global warming and the resulting long-term increase in temperatures are the main cause of heat stress in mammals [1]. Moreover, selection towards high production yield in livestock associated with high metabolic load is an additional factor that makes livestock especially prone to overheating. Heat stress negatively affects health, reproduction, and other biological functions [2]. Specifically, in dairy cattle, heat stress impedes milk production, welfare, and growth [3]. Unfortunately, the phenomenon of heat stress is common in current ages, and we should understand how its long-term susceptibility affects organisms. On the genomic level, heat stress is manifested by transcriptional and post-transcriptional regulation of heat stress-associated genes [4]. It is known that Bos indicus has greater heat tolerance than Bos taurus [5], which indicates a genetic component of heat resistance. A few genes responsible for thermotolerance in dairy cattle – HSF1, MAPK8IP1, and CDKN1B have been recently identified [6]. However, the effect of heat stress on animal-associated microbiotas is not well known. In cattle genomics, bacteria are the main cause of mastitis – one of the most prevalent diseases of dairy cattle [7]. In general, the composition of gut microbiota depends on multiple factors – genetic [8], dietary [9] and environmental [10].

Heat stress belongs to environmental factors that may change the composition of the microbiota. Studies in cattle reported microbial species which abundance depends on heat stress conditions. Chen and colleagues [11] reported the effect of heat stress on physiological characteristics and circulation levels of immune activity and the microbiome. In an experimental study Zhao and colleagues [12] identified bacterial species in the rumen microbiome associated with heat stress. In their study it was found that heat stress has no effect on both alpha and beta diversity, however the effect on the richness of microbiota was identified, especially significant increase in the abundance of Streptococcus, Enterobacteriaceae, Ruminobacter, Treponema and Bacteroidaceae. Sales and colleagues [13] in his study also reported that heat stress influenced microbiota in beef cattle rumen. Particularly, they found genera Flavonifractor, Treponema, Ruminococcus, and Carnobacterium significantly associated with heat stress. However, assessing the composition of microbiota in farm animals’ environments is important to study its association with heat stress under breeding conditions. Moreover, the categorization between heat stress and normal conditions is a simplification. The level of an animal’s heat stress is a continuous variable and thus can be assessed using quantitative metrics. This however implies non-standard statistical modeling of the association between heat stress traits and microbiome as compared to the experimental-based, case-control setup. In our study, we focused on the identification of bacteria associated with heat stress measured by drooling score, rectal temperature, and respiratory score, expressed by the estimated breeding values.

Material and methods

The material comprises 136 fecal samples of 136 Chinese Holstein cows, which were collected in 2017, 2018, and 2019 directly in herds belonging to Beijing Shounong Livestock Development Co., Ltd. The cows from the same year had been fed with exactly the same total mixed ration diet for over 1 month and the cows from different years were fed with different total mixed ration diets with small change; however, all diets were based on corn silage and concentrate, and all the cows were fed ad libitum. The experimental design deliberately did not involve formal case and control groups but was carried out on the production population.

Fecal samples were collected directly from the cow’s rectum using a method a bit similar to rectal inspection. Around 7 AM, before the new feed is provided to the cow, is the time point we selected to take fecal samples. As cows are calmer after a good rest during the night, samples are easier to keep during the morning cooler period in summer. Moreover, we can be more sure that cows are in similar digestion stage without stimulation from feed for a relatively long time, and feces accumulated in the cow’s rectum. By wearing a disposable plastic long-armed glove, the sampler inserts his hand and arm into the cow’s rectum, first removed the outer part of feces accumulated in the rectum, and then grabbed a certain amount of feces from the inner part by hand, after a few feces mixing actions in the rectum. A disposable plastic long-armed glove can be used once for each cow. After the sampler’s hand holding feces moves out of the cow’s rectum, one can turn the glove outside in. Feces will naturally accumulate into the finger parts of the glove. By cutting a small hole at the tip of the finger parts of the glove, a fecal sample can be easily transferred into a properly labeled sterile 5 ml cryopreservation tube. Since big particles within feces may precipitate at the bottom of the finger parts of the glove, the very first part of the fecal sample can be discarded. Fecal samples are normally then placed on dry ice maximum for 3-4 h before they stored at −80C at the laboratory.

The thermal environment during the sampling process was measured by temperature, humidity, and Temperature Humidity Index (THI) presented in Table 1.

Table 1 Characteristic of the thermal environment

The procedures of DNA extraction, amplification, and sequencing were completed by Wekemo Tech Co., Ltd. (Shenzhen, China). Microbial DNA was extracted from fecal samples using the E.Z.N.A. soil DNA Kit (Omega Bio-tek, Norcross, GA, U.S.) according to manufacturer’s protocols. The final DNA concentration and purification were determined by NanoDrop 2000 UV-vis spectrophotometer (Thermo Scientific, Wilmington, USA), and DNA quality was checked by 1% agarose gel electrophoresis. The V3-V4 hypervariable regions of the bacterial 16S rRNA gene were amplified with primers 338F(5’-ACTCCTACGGGAGGCAGCAG-3’) and 806R(5’-GGACTACHVGGGTWTCTAAT-3’) (for samples picked in 2017) as well as 341F(5’-CCTAYGGGRBGCASCAG-3’) and 806R(5’-GGACTACNNGGGTATCTAAT-3’) (for samples picked in 2018 and 2019) by thermocycler PCR system (GeneAmp 9700, ABI, USA). The PCR reactions were conducted using the following program: 3 min of denaturation at 95 C, 27 cycles of 30 s at 95 C, 30s for annealing at 55 C, and 45s for elongation at 72 C, and a final extension at 72 C for 10 min. PCR reactions were performed in triplicate 20 μL mixture containing 4 μL of 5 μL FastPfu Buffer, 2 μL of 2.5 mM dNTPs, 0.8 μL of each primer (5 μM), 0.4 μL of FastPfu Polymerase and 10 ng of template DNA. The resulted PCR products were extracted from a 2% agarose gel and further purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA) and quantified using QuantiFluor™-ST (Promega, USA) according to the manufacturer’s protocol. Amplicons were sequenced using the HiSeq-PE250 (samples picked in 2017) and MiSeq-PE300 (samples picked in 2018 and 2019) Illumina platforms in paired-end modes. Part of the sequence data analyzed previously by Zhang and colleagues [14] was used in this study.

All the heat stress phenotypes were measured as it was described by Luo and colleagues [15]. In particular, each lactating cow was recorded twice a day for 2 consecutive days. In order to correct for the environmental effects that may affect phenotype values, cows’ response to heat was expressed by breeding values for rectal temperature (RT), drooling (DS), and respiratory scores (RS) estimated using the following model:

$$ \begin{aligned} y_{ijklqno}\! =\! \mu + fym_{i} + p_{j} + s_{k} + m_{l} + t_{q} + thi + a_{n} + pe_{n} + \epsilon_{ijklqno}, \end{aligned} $$

where yijklqno refers to phenotype (RS, DS or RT), μ is the population mean, fymi is the fixed effect for farm-year, pj is the fixed effect of parity, sk is the fixed effect of lactation stage, ml is the fixed effect of the indication if the animal is before or after milking, tq is the fixed effect of testing time (morning or afternoon), thi is the fixed effect of temperature-humidity index, an is the animal additive genetic effect, pen is the permanent environmental effect, and εijklqno is the random residual. The covariance matrix of random effects has the following structure:

$$ var\left(\left[\begin{array}{c} a \\ pe\\ \epsilon \end{array}\right]\right) = \left[\begin{array}{ccc} A \bigotimes \sigma^{2}_{a} & 0 & 0 \\ 0 & I \bigotimes \sigma^{2}_{pe} & 0 \\ 0 & 0 & I \bigotimes \sigma^{2}_{\epsilon} \end{array}\right]. $$

The total number of 155 cows were used to estimate the breeding values. The reliability of calculate EBVs was presented in Table 2.

Table 2 Descriptive statistics of the reliability of the estimated breeding values

Furthermore, the breeding values were additionally corrected by deregression [16] in order to remove the ancestral information from the EBVs.

Processing of sequencing data

The first step of the analysis included quality control of sequenced data. For this purpose, the FastQC [17] software was used. Then, poor quality reads and adapter sequences were removed using Trim Galore [18]. Followingly, cleaned reads were processed using the QIIME 2 [19] software. First of all, data were dereplicated – reads that are 100% the same were pooled together. Next, reads were denoised – reads that occur very rarely were considered to be PCR errors and removed, as well as chimeric sequences and singletons. Those steps were done using DADA2 algorithm [20]. implemented in QIIME 2. All the sequencing runs were processed separately. Afterward, the Amplicon Sequence Variants (ASVs) table that represents counts of occurrence of a given sequence in a sample was created. Diversity within samples (α-diversity) was calculated using Simpson’s evenness and Shannon’s diversity indices using the phyloseq [21] R package. The association of microbes composition with heat stress factors was tested using aGLMM-MiRKAT test implemented in GLMM-MiRKAT R package [22].

The SILVA database (SILVA SSU 138.1) [23] was used to classify ASVs taxonomically. For the classification, the naive Bayes algorithm implemented in scikit-learn Python package was used [24].

Since taxa originally assigned by the SILVA database represent different levels of taxonomy, they were aggregated to genera and phyla levels. Genera and phyla with a variance below one and that occurred in less than three samples were excluded from downstream analyses. Filtered tables were used for the further differential abundance analysis. Additionally, for organoleptic testing of batch effect occurrent, the Uniform Manifold Approximation and Projection (UMAP; [25]) dimensional reduction technique was used to find potential sources of unwanted variability. The phylogenetic tree was generated using the align-to-tree-mafft-fasttree pipeline [26] implemented in QIIME2 software.

Differential abundance analysis

The edgeR [27] R package was used for the normalization of the processed ASV table as well as for statistical modeling of the association between the abundance of microbiota and heat stress indicators. In particular, the Trimmed Mean of M-values (TMM) based normalization [28] was applied. It identified and excluded highly abundant and highly variable genera and phyla, whereupon weighted mean of an abundance of remaining groups was used for the actual normalization [29]. The association between genera/phyla abundance and heat stress indicators was modeled using the negative binomial distribution:

$$ \begin{aligned} K &= \beta_{0} + \beta_{1}X_{1} + \beta_{2}X_{2} + e \\ \end{aligned} $$

where K represents the counts of reads for a given genus/phylum, β0 is the intercept, β1 is the effect of the DRP (expressed by log fold change), X1 is the design matrix for DRP, β2 is the effect of the sampling year class, X2 is the incidence matrix for sampling year, e is the random error.

$$ \begin{aligned} K &\sim NB(\mu_{k}, \phi) \\ \end{aligned} $$

where μk represents the mean of counts reads, and ϕ is the dispersion parameter such that \(Var(K) = \mu _{k} + \phi \mu _{k}^{2}\) calculated using Cox-Reid approximate conditional inference moderated towards the mean [30].

The significance of the effect of DRP on the relative abundance of the genus/phylum was tested using the Likelihood Ratio Test [31]. The test statistic is as follows:

$$ LR = -2ln\Bigg(\frac{L(m_{1})}{L(m_{2})}\Bigg) \sim \chi^{2}(1) $$

where m1 is the reduced model (i.e. the formula (5) without the effect of DRP), and m2 is the full model (5).

Since each genus/phylum is tested separately, multiple testing correction method was applied using a False Discovery Rate (FDR) [32].


Processing and classification of sequence variants

46,825 unique sequences of V3 and V4 regions with a total of 6,486,706 reads were identified and classified. Table 3 summarizes the classification of ASVs based on the different taxonomic levels. In general, reads were classified into two domains – archaea (0.01%) and bacteria (99,99%). Almost all reads could be taxonomically assigned up to order, but species could be assigned only to 2.35% of reads. Further analysis was carried out using genus-level and phylum-level resolution.

Table 3 Amplicon Sequence Variants classification results

Microbiota composition

The general composition of microbiota in all samples was presented on bar plot using genus-level and phylum-level resolution. Figure 1 presents the relative abundance of genera with average proportions of more than 0.5%. We can see that Clostridium is a genus with the highest relative abundance (15.14%). There were 209 genera with a relative abundance of less than 0.5%. Figure 2 presents the analogous visualization for phyla. Firmicutes is a phylum with the highest relative abundance (63.66%). Regarding the less abundant phyla, there were identified 22 phyla with less than 0.5% of the relative abundance.

Fig. 1
figure 1

The relative abundance of genera with average proportions of more than 0.5%

Fig. 2
figure 2

The relative abundance of phyla with average proportions of more than 0.5%


Genera table was then clustered using UMAP algorithm. The projection of the UMAP coordinates calculated from the ASVs counts matrix on the genus level demonstrates three distinct clusters (Fig. 3), which reflect the sampling year. In further analysis, the effect of the sampling year was corrected.

Fig. 3
figure 3

UMAP projection of the ASVs counts matrix on the genus level

Correlation analysis of diversity metrics

In order to check whether the general diversity of microbiota within samples is correlated with the DRPs, a correlation analysis was performed. Correlations were generally positive, but non-significant (Table 4) with the highest correlation estimated between DRPs for the respiratory score and the Simpson’s evenness index (0.27). Overall, non-significant correlations indicate that there is no linear dependence between DRPs and sample diversity calculated based on the abundance of genera.

Table 4 Pearson correlation coefficients between DRPs and alpha diversity measures expressed by Simpson’s evenness and Shannon diversity

Relationship of EBVs with microbiomes composition

aGLMM-MiRKAT test was performed to test the association between the microbial community composition and EBVs. None of the analyzed EBVs showed statistically significant association with the microbial composition. It means that individual genera and phyla should be considered in a statistical model.

Differential abundance analysis

Based on the results of the negative binomial model and considering FDR≤0.05 22 genera were significantly associated with rectal temperature with all but one (Helococcus) of them showing decreased abundance with the increase of rectal temperature. Rhizobium – that represents soil bacteria – was the most associated genus with the rectal temperature. The occurrence of this bacteria might be observed perhaps due to the specific metabolism or the specific plant diet.

Succinivibrio was the only genus associated with respiratory score and Pseudobutyrivibrio – with the drooling score. There was no overlap between genera significant for the three heat stress indicators (Table 5). Differential abundance analysis of phylum (Table 6) showed that 6 phyla were significiantly associated with rectal tempearture. All of them showed decreased abundance with the increase of rectal temperature. Fibrobacteres was the only phylum associated with respiratory score. Surprisingly, for drooling score, 5 differentially abundant phyla were identified. Five of them showed increase abundance with the increase of rectal temperature. Only Fibrobacteres showed the descrease abundance. Fibrobacteres was significantly associated with both drooling and respiratory scores, while Nitrospiarae, Gemmatimonadetes, Acidobacteria, and Planctomycetes were significantly associated with both rectal temperature and drooling score.

Table 5 Significant differentially abundant genera
Table 6 Significant differentially abundant phyla

In order to check the genetic relationship between those associated genera, the phylogenetic tree was created based on the 16S rRNA sequences. The genetic relationship of the associated genera was shown in Fig. 4. Colors indicate the association between the genera and the phenotype. We can see, that Succinivibrio that was associated with the respiratory score phenotypes created the single clade with all the genera associated with the rectal temperature. Only Pseudobutyrivibrio that was associated with the drooling score creates a single, separate clade.

Fig. 4
figure 4

Phylogenetic tree of genera identified in fecal samples. Color indicates significantly associated genera with a given phenotype


This study aimed to identify genera that are associated with the rectal temperature, drooling score, and respiratory score, and in the consequences, associated with heat stress. The quantitative pseudophenotypes were used in order to model animals’ microbiomes under conventional production conditions, without setting up a case (heat stress conditions) – control (standard conditions) experiment. Such an approach allows for the estimation of genera effect on heat stress under real conditions underlying dairy herd management.

The general composition of microbiota was not altered by heat stress. Therefore we focussed on single genera as potentially involved in heat stress response. Most of the genera were significantly associated with rectal temperature which might be caused by the fact that samples and measurement came from the same environment (rectum). Since most of the significantly associated genera showed decreased abundance with the increase of heat stress, we can assume, that heat stress favors the inhibition of growth of some microbial populations.

Based on the current literature, Bailey [33] observed a reduced abundance of bacteria in genus Pseudobutyrivibrio in mice exposed to stressor-induced changes. Such reduced abundance was also observed by us the association with a drooling score. Baek [34] in his study observed that Succinivibrio shows increased abundance in cows under heat stress. In our study, this genus was also associated with the respiratory score metric. Interestingly, Helcococcus, the only genus that abundance increased with increasing rectal temperature, has not been reported in studies focused on heat stress and any stress-induced conditions, but it was reported as associated with postpartum endometritis by Miranda CasoLuengo [35]. Moreover, [36] showed that Streptomyces was reported as a genus with enriched relative abundance in Jersey cows in the normal condition compared to the heat stress condition. It is worthwhile to mention that many genera reported with the association to the rectal temperatures show the high fold change, suggesting that increased rectal temperature has a high impact on microbiota composition. Proteobacteria phylum that represents most of the associated genera in our study seems to be the most important phylum in heat stress conditions. Yu [37] already reported that Proteobacteria and Firmicutes are the most common phyla associated with heat stress conditions. Interestingly, analysis based on the phylum resolution showed that there were overlapping phyla. Fibrobacteres turned out to be the significantly associated phylum with respiratory and drooling scores. This phylum was already reported as significant in heat stress analysis of pigs reported by He [38]. Chloroflexi and Planctomycetes significantly associated with rectal temperature were also reported as a significant phyla in the analysis of short-term acute heat stress on the rumen microbiome of Hanwoo steers [34].

Differences found in microbial compositions and in genera/phyla abundance suggest that those changes might occur due to adapting to climate change. In this study, the abundance of Fibrobacteres was decreased due to heat stress. The role of this bacteria is the degradation of plant-based cellulose in ruminants and acetate production. Ransom-Jones and colleagues [39] reported that glycosyl hydrolases of Fibrobacteres may produce carbohydrate activators, including cellulose enzymes and in consequence, cows may produce more energy with acetate in the rumen that can be associated with heat production. Some bacteria (e.g. Pseudobutyrivibrio) were described as a part of the microbiome, but their impact on host physiology is not yet known.

Heat stress modeled as a binary variable (i.e. normal vs. stress conditions) provides valuable insights into the understanding of the microbiome association to heat stress, however, it should be beard in mind that the real, production environment of a dairy cow markedly deviates from the experimental conditions. The most obvious differences comprise duration, intensity, and variation in ambient temperatures, which are typically not modeled in experiments. Therefore, our study, despite being more challenging from the analytical perspective, provided an attempt to analyze the microbiome dynamics directly in a production herd. In such a situation, an important aspect of the analysis is the heat stress “phenotype”. In order to pre-correct for a whole series of genetic (i.e. familial relationship) and environmental effects (such as parity or lactation stage) possibly affecting the heat stress indicator measurements, prior to the actual heat stress modeling, we decided to use breeding values as pseudophenotypes, which were then deregressed in order to remove ancestral and familiar contributions. Such an approach provided a novel approach for the investigation of bacteria in dairy cattle under heat stress condition.

Availability of data and materials

The sequence data are available in the NCBI Sequence Read Archive with accession number SRP202074 ( Other datasets generated and/or analyzed during the current study are not publicly available due to institutional constraints but are available from Yachun Wang on reasonable request.


  1. Brivio F, Zurmühl M, Grignolio S, von Hardenberg J, Apollonio M, Ciuti S. Forecasting the response to global warming in a heat-sensitive species. Sci Rep. 2019;9(1).

  2. Polsky L, von Keyserlingk MAG. Invited review: Effects of heat stress on dairy cattle welfare. J Dairy Sci. 2017; 100(11):8645–57.

    CAS  Article  Google Scholar 

  3. Nardone A, Ronchi B, Lacetera N, Ranieri MS, Bernabucci U. Effects of climate changes on animal production and sustainability of livestock systems. Livest Sci. 2010; 130(1-3):57–69.

    Article  Google Scholar 

  4. Li Q, Yang C, Du J, Zhang B, He Y, Hu Q, Li M, Zhang Y, Wang C, Zhong J. Characterization of miRNA profiles in the mammary tissue of dairy cattle in response to heat stress. BMC Genomics. 2018;19(1).

  5. Lees AM, Sejian V, Wallage AL, Steel CC, Mader TL, Lees JC, Gaughan JB. The impact of heat load on cattle. Animals. 2019; 9(6):322.

    Article  Google Scholar 

  6. Sigdel A, Abdollahi-Arpanahi R, Aguilar I, Peñagaricano F. Whole genome mapping reveals novel genes and pathways involved in milk production under heat stress in US holstein cows. Front Genet. 2019;10.

  7. Hoque MN, Istiaq A, Clement RA, Sultana M, Crandall KA, Siddiki AZ, Hossain MA. Metagenomic deep sequencing reveals association of microbiome signature with functional biases in bovine mastitis. Sci Rep. 2019;9(1).

  8. Li F, Li C, Chen Y, Liu J, Zhang C, Irving B, Fitzsimmons C, Plastow G, Guan LL. Host genetics influence the rumen microbiota and heritable rumen microbial features associate with feed efficiency in cattle. Microbiome. 2019;7(1).

  9. Scott KP, Gratz SW, Sheridan PO, Flint HJ, Duncan SH. The influence of diet on the gut microbiota. Pharmacol Res. 2013; 69(1):52–60.

    CAS  Article  Google Scholar 

  10. Karl JP, Hatch AM, Arcidiacono SM, Pearce SC, Pantoja-Feliciano IG, Doherty LA, Soares JW. Effects of psychological, environmental and physical stressors on the gut microbiota. Front Microbiol. 2018;9.

  11. Chen S, Wang J, Peng D, Li G, Chen J, Gu X. Exposure to heat-stress environment affects the physiology, circulation levels of cytokines, and microbiome in dairy cows. Sci Rep. 2018;8(1).

  12. Zhao, Min, Zheng, Wang. Effect of heat stress on bacterial composition and metabolism in the rumen of lactating dairy cows. Animals. 2019; 9(11):925.

    Article  Google Scholar 

  13. Sales GFC, Carvalho BF, Schwan RF, de Figueiredo Vilela L, Meneses JAM, Gionbelli MP, da Silva Ávila CL. Heat stress influence the microbiota and organic acids concentration in beef cattle rumen. J Therm Biol. 2021; 97:102897.

    Article  Google Scholar 

  14. Zhang G, Wang Y, Luo H, Qiu W, Zhang H, Hu L, Wang Y, Dong G, Guo G. The association between inflammaging and age-related changes in the ruminal and fecal microbiota among lactating holstein cows. Front Microbiol. 2019;10.

  15. Luo H, Brito LF, Li X, Su G, Dou J, Xu W, Yan X, Zhang H, Guo G, Liu L, Wang Y. Genetic parameters for rectal temperature, respiration rate, and drooling score in holstein cattle and their relationships with various fertility, production, body conformation, and health traits. J Dairy Sci. 2021; 104(4):4390–403.

    CAS  Article  Google Scholar 

  16. Jairath L, Dekkers JCM, Schaeffer LR, Liu Z, Burnside EB, Kolstad B. Genetic evaluation for herd life in canada. J Dairy Sci. 1998; 81(2):550–62.

    CAS  Article  Google Scholar 

  17. Andrews S. FASTQC. A quality control tool for high throughput sequence data. 2010. Accessed 21 Nov 2021.

  18. Krueger F. Trim Galore!. Accessed 21 Nov 2021.

  19. Bolyen E, Rideout JR, Dillon MR, Bokulich NA, Abnet CC, Al-Ghalith GA, et al.Reproducible, interactive, scalable and extensible microbiome data science using QIIME 2. Nat Biotechnol. 2019; 37(8):852–7.

    CAS  Article  Google Scholar 

  20. Callahan BJ, McMurdie PJ, Rosen MJ, Han AW, Johnson AJA, Holmes SP. DADA2: High-resolution sample inference from illumina amplicon data. Nat Methods. 2016; 13(7):581–3.

    CAS  Article  Google Scholar 

  21. McMurdie PJ, Holmes S. phyloseq: An r package for reproducible interactive analysis and graphics of microbiome census data. PLoS ONE. 2013; 8(4):61217.

    Article  Google Scholar 

  22. Koh H, Li Y, Zhan X, Chen J, Zhao N. A distance-based kernel association test based on the generalized linear mixed model for correlated microbiome studies. Front Genet. 2019;10.

  23. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, Peplies J, Glöckner FO. The SILVA ribosomal RNA gene database project: improved data processing and web-based tools. Nucleic Acids Res. 2012; 41(D1):590–6.

    Article  Google Scholar 

  24. Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, Blondel M, Prettenhofer P, Weiss R, Dubourg V, Vanderplas J, Passos A, Cournapeau D, Brucher M, Perrot M, Duchesnay E. Scikit-learn: Machine learning in Python. J Mach Learn Res. 2011; 12:2825–30.

    Google Scholar 

  25. McInnes L, Healy J, Melville J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. cite arxiv:1802.03426. 2018. Accessed 21 Nov 2021.

  26. Yamada KD, Tomii K, Katoh K. Application of the MAFFT sequence alignment program to large data—reexamination of the usefulness of chained guide trees. Bioinformatics. 2016; 32(21):3246–51.

    CAS  Article  Google Scholar 

  27. Robinson MD, McCarthy DJ, Smyth GK. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2009; 26(1):139–40.

    Article  Google Scholar 

  28. Robinson MD, Oshlack A. A scaling normalization method for differential expression analysis of RNA-seq data. Genome Biol. 2010; 11(3):25.

    Article  Google Scholar 

  29. Smid M, van den Braak RRJC, van de Werken HJG, van Riet J, van Galen A, de Weerd V, van der Vlugt-Daane M, Bril SI, Lalmahomed ZS, Kloosterman WP, Wilting SM, Foekens JA, IJzermans JNM, Martens JWM, Sieuwerts AM. Gene length corrected trimmed mean of m-values (GeTMM) processing of RNA-seq data performs similarly in intersample analyses while improving intrasample comparisons. BMC Bioinformatics. 2018;19(1).

  30. Cox DR, Reid N. Parameter orthogonality and approximate conditional inference. J R Stat Soc Ser B (Methodol). 1987; 49(1):1–18.

    Google Scholar 

  31. Solomon DL. A note on the non-equivalence of the neyman-pearson and generalized likelihood ratio tests for testing a simple null versus a simple alternative hypothesis. Am Stat. 1975; 29(2):101–2.

    Google Scholar 

  32. Benjamini Y, Hochberg Y. Controlling the false discovery rate: A practical and powerful approach to multiple testing. J R Stat Soc Ser B (Methodol). 1995; 57(1):289–300.

    Google Scholar 

  33. Bailey MT, Dowd SE, Galley JD, Hufnagle AR, Allen RG, Lyte M. Exposure to a social stressor alters the structure of the intestinal microbiota: Implications for stressor-induced immunomodulation. Brain Behav Immun. 2011; 25(3):397–407.

    CAS  Article  Google Scholar 

  34. Baek YC, Choi H, Jeong J-Y, Lee SD, Kim MJ, Lee S, Ji S-Y, Kim M. The impact of short-term acute heat stress on the rumen microbiome of hanwoo steers. J Anim Sci Technol. 2020; 62(2):208–17.

    CAS  Article  Google Scholar 

  35. Luengo RM-C, Lu J, Williams EJ, Miranda-CasoLuengo AA, Carrington SD, Evans ACO, Meijer WG. Delayed differentiation of vaginal and uterine microbiomes in dairy cows developing postpartum endometritis. PLoS ONE. 2019; 14(1):0200974.

    Google Scholar 

  36. Kim D-H, Kim M-H, Kim S-B, Son J-K, Lee J-H, Joo S-S, Gu B-H, Park T, Park B-Y, Kim E-T. Differential dynamics of the ruminal microbiome of jersey cows in a heat stress environment. Animals. 2020; 10(7):1127.

    Article  Google Scholar 

  37. Yu M-F, Zhao X-M, Cai H, Yi J-M, Hua G-H. Dihydropyridine enhances the antioxidant capacities of lactating dairy cows under heat stress condition. Animals. 2020; 10(10):1812.

    Article  Google Scholar 

  38. He J, Guo H, Zheng W, Xue Y, Zhao R, Yao W. Heat stress affects fecal microbial and metabolic alterations of primiparous sows during late gestation. J Anim Sci Biotechnol. 2019;10(1).

  39. Ransom-Jones E, Jones DL, McCarthy AJ, McDonald JE. The fibrobacteres: an important phylum of cellulose-degrading bacteria. Microb Ecol. 2012; 63(2):267–81.

    CAS  Article  Google Scholar 

Download references


Calculations have been carried out using resources provided by Wroclaw Centre for Networking and Supercomputing (, grant No. 509.

Computations were carried out at the Poznan Supercomputing and Networking Centre.


This work was supported by the Wroclaw University of Environmental and Life Sciences (Poland) as the Ph.D. research program "Bon doktoranta SD UPWr”.

The publication is financed under the Leading Research Groups support project from the subsidy increased for the period 2020–2025 in the amount of 2% of the subsidy referred to Art. 387 (3) of the Law of 20 July 2018 on Higher Education and Science, obtained in 2019.

Supported by China Agriculture Research System of MOF and MARA; The Program for Changjiang Scholar and Innovation Research Team in University (IRT_15R62); National Agricultural Genetic Improvement Program (2130135).

Special Item of Regional Collaborative Innovation in Tibet Autonomous Region (Grant No. QYXTZX-LS2021-01).

Author information

Authors and Affiliations



B.C., Y.W. and J.S conceived and conducted the experiment, B.C. analyzed the results and wrote the manuscript in consultation with J.S., K.W., H.L. and Y.W. All authors reviewed the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Bartosz Czech.

Ethics declarations

Ethics approval and consent to participate

The data collection process was carried out in strict accordance with the protocol approved by the Animal Welfare Committee of the China Agricultural University. All experimental protocols were approved by the Animal Welfare Committee of the China Agricultural University. All methods are reported in accordance with ARRIVE guidelines ( for the reporting of animal experiments.

Consent for publication

Not provided.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Czech, B., Szyda, J., Wang, K. et al. Fecal microbiota and their association with heat stress in Bos taurus. BMC Microbiol 22, 171 (2022).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI:


  • 16S rRNA gene
  • Heat stress
  • Fecal microbiome
  • Sequencing
  • V3-V4 regions
  • Differential abundance