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.


Introduction
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, prevalent diseases of dairy cattle [7]. In general, the composition of gut microbiota depends on multiple factorsgenetic [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 -80°C at the laboratory.
The thermal environment during the sampling process was measured by temperature, humidity, and Temperature Humidity Index (THI) presented in Table 1.
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  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: where y ijklqno refers to phenotype (RS, DS or RT), μ is the population mean, fym i is the fixed effect for farm-year, p j is the fixed effect of parity, s k is the fixed effect of lactation stage, m l is the fixed effect of the indication if the animal is before or after milking, t q is the fixed effect of testing time (morning or afternoon), thi is the fixed effect of temperature-humidity index, a n is the animal additive genetic effect, pe n is the permanent environmental effect, and ijklqno is the random residual. The covariance matrix of random effects has the following structure: The total number of 155 cows were used to estimate the breeding values. The reliability of calculate EBVs was presented in Table 2.
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: 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), X 1 is the design matrix for DRP, β 2 is the effect of the sampling year class, X 2 is the incidence matrix for sampling year, e is the random error.
where μ k represents the mean of counts reads, and φ is the dispersion parameter such that Var(K) = μ k + φμ 2 k 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: where m 1 is the reduced model (i.e. the formula (5) without the effect of DRP), and m 2 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.

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.

Clustering
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.

Correlation analysis of diversity metrics
In order to check whether the general diversity of microbiota within samples is correlated with the DRPs, a Fig. 1 The relative abundance of genera with average proportions of more than 0.5% Fig. 2 The relative abundance of phyla with average proportions of more than 0.5% 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.

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. 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.

Discussion
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.