Skip to main content
Download PDF

Causal estimation of the relationship between reproductive performance and the fecal bacteriome in cattle

Animal Microbiome volume 7, Article number: 33 (2025) Cite this article

Abstract

Background

The gut bacteriome influences host metabolic and physiological functions. However, its relationship with reproductive performance remains unclear. In this study, we evaluated the relationship between the gut bacteriome and reproductive performance in beef cattle, such as Japanese black heifers. Artificial insemination (AI) was performed after 300 days of age, and the number of AI required for pregnancy (AI number) was evaluated. The relationship of the fecal bacteriome at 150 and 300 days of age and reproductive performance was visualized using statistical structural equation modelling between traits based on four types of machine-learning algorithms (linear discriminant analysis, association analysis, random forest, and XGBoost).

Results

The heifers were classified into superior (1.04 ± 0.04 cycles, n = 26) and inferior groups (3.87 ± 0.27 cycles, n = 23) according to the median frequency of AI. The fecal bacteria of the two groups were examined and compared using differential analysis, which demonstrated that the genera Rikenellaceae RC9 gut group and Christensenellaceae R-7 group were increased in the superior group. Subsequently, correlation analysis evaluated the interrelationships between bacteriomes, which demonstrated that the patterns exhibited distinct characteristics. Therefore, four machine-learning algorithms were employed to identify the distinctive factors between the two groups. The directed acyclic graphs carried out by DirectLiNGAM based on these extracted factors inferred that the family Erysipelotrichaceae and the genera Clostridium sensu stricto 1 and Family XIII AD3011 group at 150 days of age were strongly associated with an increase in AI number. Furthermore, a pathway involved in creatinine degradation (PWY-4722) at 150 days of age was related to an increase in AI number. However, bacteriomes and/or pathways at 300 days of age were not necessarily related to AI number.

Conclusions

In this study, a causal inference methodology was applied to investigate AI-dependent gut bacterial communities in pregnant cattle. These findings suggest that AI numbers, which are crucial for beef cattle production management, could be inferred from the fecal bacterial patterns nearly six months before the AI, rather than immediately before. This study provides a novel perspective of the gut environment and its role in reproductive performance.

Graphical Abstract

Background

To align with the Sustainable Development Goals (SDGs) [1], environmentally sustainable practices must be implemented in contemporary livestock production. The reduction of feed and production costs is imperative for environmentally sustainable livestock management. In order to achieve this, it is necessary to improve reproductive performance and reduce the period during which feed is fed solely for the growth and maintenance of female cows [2]. Moreover, the early exclusion (shipping) of unsuitable individuals from the reproductive programme in heifers will eliminate the need to feed the heifers waste feed. Consequently, enhancing reproductive performance is critical for the long-term sustainability of the livestock industry.

In general, its well-known that reproductive traits have a low heritability due to environmental factors. Previous studies have reported that the heritability of fertility is low in pigs [3], dairy cattle and beef cattle [4, 5]. This may be attributed to the fact that fertility is largely influenced by environmental factors rather than genetic factors. For example, excessive obesity has been shown to reduce fertility in humans [6] and ruminants [7]. Recently, high plane of nutrition programs in which a large amount of milk replacer (MR) is fed to calves has become increasingly popular. This has been shown to improve calf growth [8, 9]. However, there is concern that excessive amounts of MR may cause over-fattening of calves and reduce reproductive performance. Indeed, Taguchi et al. reported that differences in the amount of MR fed during the suckling period affected the concentrations of nutritional metabolites in the plasma and reproductive performance at the first lactation of heifers [10]. In addition, it has been reported to have adverse effects on fertility in mice with inactivated immune cells [11] and in cattle with immune dysfunction [12]. As such, the reproductive performance of cattle is likely to be greatly affected by their nutritional status and immune system. In humans, it has been established that alterations in gut microbiota are associated with obesity [13] and immune system [14]. In addition, in a microbiological point of view, there is a belief that symbiotic environmental microorganisms affect the reproduction of animals [15, 16]. Therefore, the gut bacteriome plays a role in reproduction [16,17,18] via the alteration in nutrient metabolism and immune system [19]. The role of the gut bacteriome in reproductive performance is a topic of ongoing research [20]. In recent years, the influence of gut microbiota on reproductive dysfunction and the viewpoint of the gut microbiota-gonadal axis have attracted attention [17, 18]. In ruminants, it has been suggested that bovine vaginal and fecal microbiome could be used as biomarkers of bovine reproduction [21]. These findings suggest the importance of macroscopic assessment of the characteristics of the gut microbiota of cattle before conception. However, the relationship between the gut microbiota during the preinsemination period and reproductive performance is not yet fully understood.

This study aimed to evaluate a relationship between the gut bacteriome and reproductive performance for Japanese black cattle, as a model of industrial animals. We conducted a comprehensive analysis of bacteriomes in the feces of heifers with different genetic backgrounds at 150 and 300 days of age (Fig. 1). It was hypothesized that this would allow us to evaluate some of the effects of symbiotic microorganisms that do not depend on genetic factors using heifers with different genetic backgrounds (Father, Grandfather, and Great Grandfather, Table S1). Artificial insemination (AI) was performed after 300 days of age and continued until pregnancy. The relationships between the number of AI required for pregnancy (AI number) and bacteriomes were examined via the machine learning (ML) algorithms and statistical causal inference. Based on these results, we estimated the causal structure of the bacterial community related to the reproductive performance of the heifers. The findings of this study offer insights into the prediction of livestock reproductive performance and efficient breeding management.

Fig. 1
figure 1

Experimental design of heifers from birth to parturition. The heifers in the experiment were fed a TMR throughout the study period, concentrate up to 300 days of age, milk replacer up to 90 days of age (days to weaning), and hay from weaning to 300 days of age. Following the collection of fecal samples, correlation analysis and machine learning were conducted based on the bacterial population data. Subsequently, causal inference was performed. Artificial insemination (AI) during gestation period was performed only on heifers meeting the AI physique criteria: body weight ≥ 270 kg and body height ≥ 116 cm [10]

Full size image

Results

Physical and reproductive performance

Table 1 shows the growth and reproductive performance of each group at 300 days of age. There was no significant difference in body height (BH) or body length (BL) between the groups at p < 0.05, and body weight (BW) tended to be greater in the superior group (p = 0.07). However, there was no difference (p = 0.44, p = 0.55, p = 0.44, p = 0.85, p = 0.16, and p = 0.38; Tables S2–S7) in BW between the groups when the feeding patterns of the different MRs were compared. These results pointed out that the two groups with different reproductive performance were not affected by a special feeding pattern. This also means that it is difficult to evaluate the difference with a small number of animals. The AI number was significantly lower (p < 0.01) in the superior group (1.04 ± 0.04) than in the inferior group (3.87 ± 0.27). The breakdown revealed that 65.22% of inferior heifers became pregnant after four or more AI cycles, whereas 96.15% of the superior heifers became pregnant during their first AI cycle. The results revealed that the superior group had a successful AI at 339.26 ± 5.98 days of age, approximately 130 days earlier than the inferior group did. The two groups presented significantly different reproductive performance, independent of growth performance. Therefore, assuming that factors other than heifer growth influence reproductive performance, we evaluated fecal bacteriomes as follows.

Table 1 Growth and reproductive performance of the heifers with superior or inferior reproductive performance
Full size table

Relative abundance and diversity analysis of the fecal bacteriomes

The alpha and beta diversity of fecal bacteriomes at 150 and 300 days of age for each group are shown in Fig. 2. Chao1 values and beta diversity were not significantly different between groups (Fig. 2a and c). In contrast, Simpson values in the superior group tended to be higher compared with the inferior group at 150 days of age (Kruskal p < 0.01, Dunn p = 0.02, Fig. 2b). Although there was no significant difference in fecal bacteria at the phylum level between the groups (Fig. 3a), the genus Rikenellaceae RC9 gut group at 150 days of age and the genus Christensenellaceae R-7 group at 300 days of age were had an increased abundance in the superior group than in the inferior group (p = 0.02 and 0.04, respectively; Fig. 3b). To identify differences in the relationships between the fecal bacteria of heifers within a group at different ages, a correlation analysis was conducted on the bacterial genera with relative abundances greater than 10% between 150 and 300 days of age (Fig. 4). A total of eighteen × eighteen relationships were observed between 150 and 300 days of age in the groups, with no differences in correlations for the same bacterial genera between the groups. However, among the different bacteria, the correlation coefficients for the superior and inferior groups were significantly different (p < 0.05) for nineteen relationships (Fig. 4 and Table S8). These results indicate that the relationships between days of the same bacteria are consistent across groups, whereas the relationships between separate bacteria may vary. Thus, the interrelationships of bacteria differed between the groups.

Fig. 2
figure 2

Fecal bacterial diversity in each stage. The values of (a) Chao1 and (b) Simpson as α-diversity of the bacteriome in the feces of heifers with superior or inferior reproductive performance were shown. Data are presented as the mean ± standard error of the mean. Asterisks (*) indicate significant differences between groups (p < 0.05). c The plot for non-metric multidimensional scaling (NMDS) as evaluation for β-diversity were visualized. Abbreviations are as follows: 150d_Sup (pink), superior group in 150d of the age; 150d_Inf (green), inferior group in 150d of the age; 300d_Sup (orange), superior group in 300d of the age; 300d_Inf (light green), inferior group in 300d of the age. The table shows the statistical values calculated using the R library ‘pairwiseAdonis’

Full size image
Fig. 3
figure 3

Relative abundances of (a) phyla and (b) genera in the feces of heifers with superior or inferior reproductive performance. Asterisks (*) indicate significant differences between groups (p < 0.05)

Full size image
Fig. 4
figure 4

Correlation heatmaps of fecal bacterial genera between 150 and 300 days of age in heifers with superior or inferior reproductive performance. Eighteen bacterial genera with relative abundances greater than 10% were selected for analysis between 150 and 300 days of age. Yellow squares indicate the correlation coefficients that were significantly different between the superior and inferior groups

Full size image

Bacterial classification by ML

Although the differential analysis captured only the changes in a single taxon, the results of the correlation analysis at 150 and 300 days of age indicated that the relationships between the bacteria differed. Accordingly, given the potential for differences in the interrelationships among various bacterial groups, ML was employed to facilitate the classification of these groups. We used the bacterial data on the presence ratio of 1% or more for these analyses.

First, linear discriminant analysis (LDA) of the physical score and fecal bacterial taxa at 150 and 300 days of age was performed, as shown in Fig. 5, to select factors that may influence reproductive performance. LDA revealed that the families Ruminococcaceae and Rikenellaceae, and the genus Rikenellaceae RC9 gut group increased, and the genus Holdemanella decreased significantly in the superior group compared with the inferior group at 150 days of age (Fig. 5a and b). The family Christensenellaceae, the genus Christensenellaceae R-7 group, the family Barnesiellaceae, and the genus Bacteroidales RF16 group increased in the superior group compared to the inferior group at 300 days of age (Fig. 5c and d).

Fig. 5
figure 5

Screening for potential relationships among factors via linear discriminant analysis (LDA). The LDA effect size (LEfSe) cladogram was visualized on the basis of the significant changes in the bacterial population calculated by LDA (p < 0.05; > twofold change) in heifers at (a) 150 and (c) 300 days of age. Significant changes in the bacterial population were calculated using LDA (p < 0.05; > twofold change) in heifers at (b) 150 and (d) 300 days of age. The sharp (#) indicates bacterial genera whose abundance increased in the superior group in Fig. 3. The cross (♰) indicates bacterial genera with significantly different correlation values in Fig. 4

Full size image

Next, bacterial feature factors in 150 and 300 days of age were selected via three types of ML algorithms based on the workflow shown in Fig. S1a. The results of the association analysis (AA) revealed an increase in six bacterial taxa and a decrease in eight taxa in the superior group (Fig. S1b), and an increase in ten taxa and a decrease in six taxa in the inferior group at 150 days of age (Fig. S1c). At 300 days of age, AA revealed an increase in fifteen taxa and a decrease in five taxa in the superior group (Fig. S1d), and an increase in the five taxa and a decrease in the fourteen taxa in the inferior group (Fig. S1e). As a result, the factors selected by LDA were included in AA, and feature factors at 150 and 300 days of age were not necessarily the same. To evaluate the feature values in detail, Random Forest (RF) and XGBoost (XGB) were used to calculate each feature of bacteria selected by LDA and AA in the group (Fig. 6). The RF and XGB can also be applied to classifying and extracting groups of important factors to classify and discriminate between multigroup [22, 23]. First, the 150d dataset was used as training data to analyze the prediction accuracy of 300d data as test data. As the result, as described in Fig. 6a, accurate prediction could not necessarily be made, and the error rates were 13.04% (300d_Inf) and 80.77% (300d_Sup) in RF, respectively (Prediction table I in Fig. 6a). In the XGB prediction, the classification of superior group (Sup) and inferior group (Inf) was successful, suggesting that those of 150d and 300d could not discriminate in the conditions. Although it is possible that the XGB algorithm identifies two groups (superior and inferior) from 150 days of data as a training dataset, the 150d and 300d data were discriminated by other LDAs and association analysis (Figs. 5 and S1). Therefore, in order to accurately obtain the feature factors at each age, we tried to utilize the entire dataset of each age. In other words, machine learning algorithms were here applied for the purpose to accurately select features for categorization under this experimental condition from small dataset, rather than predicting the similarity of test datasets from a huge amount of training dataset. The feature factors were selected under the calculation conditions (training and test data set is same) that achieved accuracy in distinguishing between two groups (superior and inferior) of each age stage. Therefore, the extraction of these factor groups is not intended to have a highly universal prediction probability, but the features under the experimental conditions are exactly selected. At 150 days of age, seventeen taxa were selected using RF (Fig. 6b), but XGB focused fourteen taxa of them as follows: the orders Bacteroidales and Corynebacteriales, the families Corynebacteriaceae, Erysipelotrichaceae, Lachnospiraceae, Rikenellaceae, and Ruminococcaceae, and the genera Alistipes, Clostridium sensu stricto 1, Coprococcus 3, [Eubacterium] coprostanoligenes group, Family XIII AD3011 group, Rikenellaceae RC9 gut group, and [Ruminococcus] torques group (Fig. 6d). At 300 days of age, nineteen taxa were selected using RF (Fig. 6c), but XGB focused twelve taxa of them as follows: the class Bacilli, the order Lactobacillales, the families Barnesiellaceae, Erysipelotrichaceae, and Family XIII, and the genera Bacteroidales RF16 group, Catenisphaera, Christensenellaceae R-7 group, Clostridium sensu stricto 1, Holdemanella, Lachnospiraceae AC2044 group, and Treponema 2 (Fig. 6e). These results allowed us to select characteristic bacteria (including the minority of bacteria, not just the predominant data on the presence ratio of 10% or more) in the superior and inferior groups that could not be detected by differential or correlation analysis.

Fig. 6
figure 6

The feature importance ranking of bacterial taxa in heifers with superior or inferior reproductive performance by the random forest (RF) and XGBoost (XGB) algorithms. a Treatment of training dataset and test dataset. Here, the real number datasets of AA-selected factors for 150d and 300d (Fig. S1) are prepared for RF and XGB. Since the AA-selected factors are slightly different between the 150d and 300d datasets, factors detected in at least one of these datasets are also included. The confusion matrix was shown using the 150d and 300d data as the training and test data sets, respectively. The tables (Prediction Tables I and II) show a confusion matrix via RF and XGB, and class. error in the table shows the error rate in each class of the test data set. Since the accuracy is not high, it can be shown below that an attempt is made to detect features using the whole data. The abbreviations are as follows: training dataset, a dataset for training via the MLs for a model; test dataset, a dataset for testing via the MLs to be predicted; 150d, 150 days of age after birth; 300d, 300 days of age after birth; Inf, an inferior group; Sup, a superior group; 150d_Inf, the data of inferior group of 150d; 150d_Sup, the data of superior group of 150d; 300d_Inf, the data of inferior group of 300d; 300d_Sup, the data of superior group of 300d; class. error, an error rate in each class in the table. The relative levels of the RF feature is shown in (b) 150d and (c) 300d, and the levels of the XGB feature can be found in (d) 150d and (e) 300d. The X-axis indicates the feature values (MeanDecreaseGini, a RF feature unit in (b) and (c); and Importance, a XGB feature unit in (d) and (e)). The Y-axis indicates the bacterial taxa screened by LDA and AA at 150 and 300 days of age. The bacterial taxa classified using LDA (Fig. 5) and AA (Fig. S1) are as follows: red bars, bacterial taxa classified into superior groups; green bars, bacterial taxa classified into inferior groups. The factors selected by the LDA have been underlined. The factors with dotted corrals in (d) and (e) show feature in case that the 150d and 300d datasets are used as the training and test dataset, respectively

Full size image

Causal inference for the ML-selected bacteriome

To estimate the causal structure between reproductive performance and fecal bacteriomes, the direct method for linear non-Gaussian acyclic model (DirectLiNGAM), which estimates the causal structure through independent component analysis and can handle different relationship between traits, was performed on bacterial taxa selected by ML. DirectLiNGAM is used to predict the causal and/or essential structure of a set of factors that do not necessarily have a normal distribution (Gaussian distribution). It is known that the method provides important information for understanding the potential highly association between factors and for structural inference as a group of important factors, not an independent factor, or not necessarily causal relationship [24]. The analysis revealed no relationship between AI number and the taxa selected at either the combined data of 150 and 300 days of age, nor the data only 300 days of age (Figs. S2 and S3). In contrast, DirectLiNGAM from the data of only 150 days of age revealed that the family Erysipelotrichaceae and the genera Clostridium sensu stricto 1 and Family XIII AD3011 group had a direct causal structure with AI number of 6.08, 9.45, and 157.88, respectively (Fig. 7a). The calculation internal standard at this time was 1.00 for the values between the family Rikenellaceae and the genera Rikenellaceae RC9 gut group and Alistipes, i.e. all of these bacteria were not inconsistent with the fact that they all belong to the family Rikenellaceae. The statical value of DirectLiNGAM for Figs. 7a, S2, and S3 was shown in Tables S10–12. In general, DirectLiNGAM also indicates that the relationship is weak if the causal value is less than 1 and strong if it is greater than 1 [23,24,25]. Based on the hypothesis that it is possible to construct a structural model for three essential bacteria (family Erysipelotrichaceae, genus Clostridium sensu stricto 1, and Family XIII AD3011 group), estimated by DirectLiNGAM that have a strong contribution to the AI number, the group comprising these four components was also validated via structural equations (Fig. S4). This calculation was performed using a robust estimation method that can accommodate non-Gaussian distribution data, and a fitting model was obtained to show the importance as a group. This model showed suitable values for the indices (chi-square statistic (Chisq) = 1.161, p = 0.28; robust comparative fit index (CFI.robust) = 0.99; root mean square error of approximation (RMSEA) = 0.00; standardized root mean square residual (SRMR) = 0.04; goodness-of-fit index (GFI) = 0.993). Furthermore, differential analyses comparing these important factors demonstrated no statistical significance between the two groups at corresponding age of day, yet statistical significance was observed between different age of day (Fig. 7b). The results of differential analyses of the other bacterial taxa were shown in Fig. S5. Based on these results, a group of that are computationally feature factors under the experimental conditions in this study have been estimated at 150 days of age.

Fig. 7
figure 7

The calculated causal relationship of the bacteria strongly associated with reproductive performance was visualized via linear non-Gaussian acyclic model (DirectLiNGAM) analysis. a The arrow indicates the trend of the causal contribution. The number indicates the value of the causal contribution calculated via DirectLiNGAM analysis. The negative and positive values indicate the negative and positive causal contributions, respectively. “AI number” is displayed in underlined black bold type. Bacteria written in red indicate bacteria for which a direct relationship with “AI number” was observed. The sharp (#) indicates bacterial genera whose abundance increased in the superior group in Fig. 3. The cross (♰) indicates bacterial genera with significantly different correlation values in Fig. 4. b Relative abundance of bacterial taxa directly associated with AI number. The data are presented as the mean ± standard error of the mean. Relative abundances of the other bacterial taxa are shown in Fig. S5

Full size image

Causal inference for the metabolic pathways predicted using Picrust2

Figure 8a and b show the pathways predicted using the 16S rRNA sequence data from the fecal bacteria with significant changes in the superior and inferior groups, respectively, via volcano plots. The pathways selected for a false discovery rate (FDR) < 0.1 at 150 days of age were PWY-6906, PWY-622, VALDEG-PWY, and PWY-4722, and those at 300 days of age were PWY-5178, PWY-5430, PWY-5266, PWY-5273, PWY-6397, PWY-6383, and PWY-1G-0 (Fig. 8a and b, Table S9). We then estimated the causal structure of the pathways selected from combined data of 150 and 300 days of age and AI number (Fig. 8c). The results of the Shapiro–Wilk test revealed that all these factor groups were statistically non-normally distributed (p < 0.05). Therefore, assuming a non-normal distribution, DirectLiNGAM was applied to AI number. The internal standard at this time was − 1.00 for the values of “150d” and “300d,” which exhibited precisely the opposite characteristics [24]. AI number had a strong positive relationship with PWY-6383 and PWY-4722 (coefficient = 67,441.82, 223,231.81) and a strong negative relationship with PWY-5178 and PWY-6397 (coefficient = − 462,237, − 100,669). Furthermore, differential analyses comparing these pathways showed that no statistically significant variations were identified between the two groups at corresponding age of day; however, statistically significant variations were identified between different age of day (Fig. 8d). The results of the analysis by each day of age revealed that PWY-4722 had a direct causal structure with AI number at 150 days of age (coefficient = 288,852.57, Fig. S6), whereas no causal structure with AI number was found at 300 days of age. The statical value of DirectLiNGAM for Figs. 8c and S6 was shown in Tables S13-14.

Fig. 8
figure 8

Pathways selected in heifers with superior or inferior reproductive performance. Volcano plot for pathways selected after preprocessing with p-adjusted cut-off of 0.05, as shown in (a) 150 and (b) 300 days of age. Pathways whose expression was significantly increased or decreased (p < 0.05, FDR < 0.1) are shown in red or green, respectively. c The calculated causal relationships of the pathways strongly linked with AI number were visualized via linear non-Gaussian acyclic model (DirectLiNGAM) analysis. The arrow indicates a causal relationship trend. The number shows the value of the causal contribution calculated via the DirectLiNGAM analysis. The negative and positive values represent the negative and positive causal contributions, respectively. d Relative abundance of pathways screened by DirectLiNGAM analysis

Full size image

Discussion

Although the gut bacteriome influences the metabolic and physiological functions of the host and is thought to affect reproductive performance in monogastric animals [20], its relationship with reproductive performance in cattle is unclear. In this study, we comprehensively evaluated the relationship between the fecal bacteriome at 150 and 300 days of age and reproductive performance after 300 days of age. Using ML and causal inference, the characteristics of the fecal bacteria that influence reproductive performance were estimated as follows.

In general, the first AI is performed once the heifers have achieved the physical criteria set by each farm to avoid calving accidents and dystocia. In addition, previous studies have reported a positive curvilinear relationship between BW and reproductive performance in dairy heifers [26]. However, in this study, although there was no significant difference in BW, BH, or BL (p > 0.05; Table 1), AI number was significantly lower in the superior group, and the days of successful AI were significantly earlier than those in the inferior group were (Table 1). Therefore, this study investigated the relationship between the gut bacteriome and reproductive performance among heifers with the same pre-breeding physical size but different pregnancy success rates per AI.

First, to evaluate whether the fecal bacterial population differed in terms of reproductive performance, we analysed bacterial diversity and the relative abundance in the fecal bacteriome at 150 days of age (the rearing period) and 300 days of age (the preinsemination period) (Fig. 2). The results revealed that the fecal bacterial diversity was greater in the superior group than in the inferior group. In addition, the relative abundance of the genus Rikenellaceae RC9 gut group at 150 days of age and the genus Christensenellaceae R-7 group at 300 days of age was greater in the superior group than that in the inferior group (Fig. 3). Fecal diversity has been shown to be lower in obese patients than in healthy individuals [27]. Interestingly, elderly Italians with a greater abundance of the families Rikenellaceae and Christensenellaceae in their fecal bacteriomes presented a reduction in BW, body mass index, visceral adipose tissue, and subcutaneous adipose tissue [28]. Obesity negatively affects the reproductive performance of monogastric [6] and ruminant animals [7]. Although the results of this study did not reveal significant differences in physical size at 300 days of age between the groups (Table 1), the increase in these fecal bacteria suggests an obesity status that cannot be determined from an external appearance and may have affected subsequent reproductive performance, i.e., AI number. Furthermore, correlation analysis among bacteria whose relative abundance was greater than 10% at the genus level revealed that some fecal bacterial relationships differed between the two groups (Fig. 4 and Table S8).

Therefore, to further characterize the relationship between reproductive performance and fecal bacteriome, we performed ML and causal inference. As a result, we were able to select a characteristic group of bacterial taxa that could not be identified via differential analysis as candidates for bacteria that may have a positive or negative effect on reproductive performance. We performed RF and XGB analyses for the bacterial taxa selected using LDA (Fig. 5) and AA (Fig. S1). The RF results revealed that seventeen and nineteen taxa were characteristic of the superior and inferior groups, respectively, at 150 days of age and 300 days of age (Fig. 6b and c). Further refinement by XGB indicated that fourteen and twelve taxa were selected at 150 and 300 days of age, respectively (Fig. 6d and e). DirectLiNGAM was performed to infer the statistical causal structure as a structural equation modelling (Fig. 7a). The essential structure estimated by DirectLiNGAM was consistent with the results of ML (Figs. 5, 6, and S1) and a structural equation model with the robust estimation (Fig. S4). Specifically, at 150 days of age, the family Erysipelotrichaceae and the genera Clostridium sensu stricto 1 and Family XIII AD3011 group presented strong causal structures with AI number (coefficients = 6.08, 9.45, and 157.88, respectively), computationally indicating that these bacteria increase AI number. The genus Clostridium sensu stricto 1 has been reported to be positively correlated with BW and serum lipid indices in Chinese adults [29]. The family Erysipelotrichaceae has been linked to weight gain and an increased risk of developing obesity, cancer, and autoimmune disorders of the nervous system [30, 31]. In addition, the genus Family XIII AD3011 group, which was selected as a negative factor at 150 days of age (Fig. 6b and d), is known to be involved in inflammation and disease owing to abnormalities in the gut bacteriomes [32]. Therefore, an increase in these bacterial taxa may negatively affect reproductive performance due to metabolic diseases, inflammation, and gut bacteriome-derived diseases. The genus Rikenellaceae RC9 gut group exhibited a higher relative abundance in the superior group compared to the inferior group at 150 days of age (Fig. 3) and included the indirect causal relationship with AI number (Fig. 7a), indicating that Rikenellaceae RC9 gut group may be associated with a reduction in AI number. Furthermore, the genus Clostridium sensu stricto 1 at 150 days of age, which demonstrates a strong causal structure with AI number (Fig. 7a), exhibits a significantly disparate correlation coefficient with the genus Lachnospiraceae NK3A20 group at 300 days of age between the groups (Fig. 4 and Table S8). Considering the absence of a causal structure between the genus Lachnospiraceae NK3A20 group at 300 days of age and AI number (Fig. S3), as well as the lack of a significant difference in the abundance of the genus Lachnospiraceae NK3A20 group at 300 days of age between the groups, it can be postulated that the genus Clostridium sensu stricto 1 at 150 days of age may be related to an increase in AI number. Interestingly, the family Ruminococcaceae, an essential factor selected in this study, is associated with reduced visceral fat in humans [28]. Notably, a direct causal relationship was not observed between the administration of the probiotic thermophile Caldibacillus hisashii during the pre-weaning period (at least from 4 to 90 days of age) and AI number, the family Erysipelotrichaceae and the genera Clostridium sensu stricto 1 and Family XIII AD3011 group (AI number-regulated core bacterial group). As previously reported [33], the administration of C. hisashii period modulates the fecal bacterial populations of post-weaning calves. In the present study, the calculated results for the conditions under which C. hisashii was introduced for a short period of time in the pre-weaning period excluded C. hisashii (Fig. S7 and Table S15) from the causal structure around the AI number in Fig. 7. In other words, it is possible that the short-term administration of C. hisashii only in the pre-weaning period has no impact on the relationship for the AI number-regulated core bacterial group; and then, the presence of a potential pre-weaning probiotic-independent central relationship between reproductive performance and the AI number-regulated core bacterial group was suggested. Based on these results, it is possible that the bacterial taxa selected by ML and causal inference influence reproductive performance through the metabolic state and gut environment of the host. Furthermore, it may be possible to estimate subsequent reproductive performance by identifying the proportion of these bacterial taxa present at 150 days of age. On the other hand, no relationship with AI number was observed for the bacterial taxa at 300 days of age. Thus, it was predicted that the bacteriome at 150 days of age would be more likely to have a strong association with AI number than that at 300 days of age.

Picrust2 was implemented to predict pathway data based on the sequence data of the bacterial 16S RNA genes. We compared pathway data between the two groups to evaluate the metabolic function of the fecal bacteriome. Volcano plot analysis revealed the importance of four pathways with significant differences at 150 days of age, and nine pathways differed between the two groups at 300 days of age (Fig. 8a and b, Table S9). The statistical causal inferences were performed by DirectLiNGAM for these pathways and AI frequency for the 150- and 300-day data, and the pathways that showed strong positive causality with AI number were PWY-6383, which synthesizes components of the mycobacterial cell wall (mono-trans and poly-cis decaprenyl phosphate), and PWY-4722, which degrades creatinine. In addition, the pathways that showed a negative correlation with AI number were PWY-5178, which degrades toluene, and PWY-6397, which biosynthesizes the mycolyl-arabinogalactan-peptidoglycan complexes necessary for the survival of Mycobacterium tuberculosis (Fig. 8c). Causal inference was then performed by day of age, and it was estimated that PWY-4722 was strongly associated with an increase in AI number at 150 days of age (Fig. S6). Since metabolic pathways in the gut bacteriome generally function when the starting factor for the pathway is present or when the bacteria responsible for metabolism are present, it is possible that fecal creatinine is reduced in the superior group. In humans, concentrations of creatinine in the blood decrease after pregnancy [34]. In addition, the Suffolk breed of sheep, which has a low reproductive capacity, has a low ability to excrete creatinine [35]. Therefore, creatinine is thought to have a negative effect on the maintenance of pregnancy. In other words, the reduced fecal creatinine in the superior group may have resulted in the reduced expression of PWY-4722, allowing fewer AI number. Based on the results of the predicted metabolic pathway, it is possible to estimate the subsequent reproductive performance of heifers by examining the expression levels of creatinine and PWY-4722 at 150 days of age. In contrast, no metabolic pathways associated with AI number were observed in the causal inference at 300 days of age or in the bacterial analysis. In addition, the bacteriomes and metabolic pathways predicted by fecal bacteria at 150 days of age, compared with those at 300 days of age prior to AI, are more strongly associated with an increase in AI number. Thus, these observations suggest that differences in the gut bacteriome of heifers and their predicted metabolic pathways may also affect the functional aspects of the host, with potential consequences for subsequent reproductive performance.

In this study, multiple ML algorisms followed by structural equation modelling were used to narrow down the important factor groups, and grouping with a high accuracy rate was made possible at the bacterial level and at the metabolic pathway level based on bacterial data. These results do not always indicate accurate causal relationships themselves, but essential groups that be beyond the increase or decrease of individual factors. Needless to say, since these are only computational features under the experimental conditions of this study, further experimental verification is necessary. If this experimental condition is considered as an aid to universal evaluation, this study suggested that it is possible to select breeding female heifers that are more likely to be fertile from the gut bacteriomes at an early stage. This approach could facilitate more efficient management of livestock. The results of this study will contribute to improvements in animal welfare, industrial animal management, and animal science in the future.

Conclusions

The present study successfully demonstrated that the fecal bacteriomes of heifers prior to AI were associated with AI number, using Japanese black cattle as a model of reproductive animals. To the best of our knowledge, this is the first study to use ML and causal inference to analyse the relationship between the fecal bacteria of heifers and subsequent reproductive performance and to extract bacterial taxa that may influence fertility. The results of this study suggest that gut bacteriomes during the postweaning growing phase are more strongly associated with reproductive performance than those just prior to AI. The results indicate that it may be possible to predict subsequent fertility at least five months before the first AI. On the other hand, these results suggest that rather than a single bacterium affecting reproductive performance, multiple groups of bacteria act together to alter reproductive performance. Consequently, it is essential to develop a methodology that holistically assesses fluctuations in the abundance of multiple bacterial species rather than focusing on a single bacterial strain to establish a reliable approach for predicting reproductive performance. Therefore, further research should be conducted with different target animals under different environmental conditions. This research approach provides a foundation for the development of new technologies to increase the productivity of livestock animals by elucidating the interrelationship between the gut bacteriome and reproductive performance and by predicting subsequent fertility.

Methods

Animals and sample collection

All heifers used in this study were reared at a commercial farm in Miyazaki, Japan, which specializes in the production of Japanese black cattle. The experiments were performed as shown in Fig. 1. Ninety-three heifers were separated from dams at 4 days of age and administered MR (DM877, Dotoh Shiryo Co. Ltd., Hokkaido, Japan) from 4–90 days of age. All heifers were randomly assigned to six MR feeding plans (P1 and P2, maximum 1,600 g/d × 40 d; P3, maximum 1,800 g/d × 60 d; P4, maximum 1,800 g/d × 41 d; P5, maximum 1,400 g/d × 40 d; and P6, maximum 1,360 g/d × 30 d). Heifers were allowed to freely feed on the following total mixed ration (TMR): TMR-A [20.0% (w/w) crude protein (CP), 2.0% (w/w) crude fat, 75.0% (w/w) total digestible nutrients (TDN)], TMR-B [12.5% (w/w) CP; 2.0% (w/w) crude fat; 60.0% (w/w) TDN], TMR-C [17.0% (w/w) CP, 2.0% (w/w) crude fat, 71.0% (w/w) TDN], and TMR-D [10.9% (w/w) CP; 1.6% (w/w) crude fat; 60.4% (w/w) TDN]. The TMRs were fed at 31–150 days, 151–210 days, 211–540 days of age, and 541 days of age and older. Concentrates [17% (w/w) CP, 2% (w/w) crude fat, and 71% (w/w) TDN] were fed ad libitum at 181–300 days of age, while oats and timothy were supplemented at 0.1–0.5 kg/head/d from 91–300 days of age. In addition, the heifers from P2, P3, P4, and P5 were supplemented with Caldibacillus hisashii from 4–90 days of age during preweaning period, as described previously [33]. Body weight, BH, and BL were measured at 305.3 ± 0.8 days of age. Artificial insemination was only performed on heifers that met the physique criteria of BW > 270 kg and BH > 116 cm [10]. Artificial insemination was repeated until pregnancy was confirmed. Heifers with different genetic backgrounds were selected for this study (n = 49, Table S1). The experimental groups were assigned according to AI number as follows (Tables S2–S7): superior (n = 26), i.e. P1 (n = 7), P2 (n = 4), P3 (n = 4), P4 (n = 4), P5 (n = 4), and P6 (n = 3) as the number of heifers that required less than the median AI number (one AI when the median was one) for each MR feeding pattern; inferior (n = 23), i.e. P1 (n = 4), P2 (n = 4), P3 (n = 4), P4 (n = 4), P5 (n = 4), and P6 (n = 3) as the number of heifers that required a median or greater AI number. To reduce the effects of different MR feeding regimens, at least six heifers (three heifers per group) were selected from each feeding regimen. When there were four or more heifers with the same AI score, heifers with the same pedigree as the other group were selected preferentially to prevent variation in the genetic background between the superior and inferior groups. Fecal samples from heifers were collected at 10:00 AM (at least 2.5 h after the morning meal) at 150 and 300 days of age via a fecal collection tube (Sarstedt AG & Co. KG, Germany). All fecal samples were stored at − 30 °C until further analysis. All data collection and sampling were performed as a routine part of daily farm management to monitor the health status of the cattle.

Meta-sequence analysis of bacterial 16S rRNA gene sequences

DNA was extracted from the fecal samples of heifers in the superior and inferior groups via the QIAamp PowerFecal DNA Kit (QIAGEN N.V., Inc., Hilden, Germany) according to the manufacturer’s protocol. DNA concentrations were evaluated via a Quant-iT™ PicoGreen dsDNA Assay Kit (Thermo Fisher Scientific, Tokyo, Japan). The nucleotide sequence of the V3-V4 region (314F-806R) of the bacterial 16S rRNA gene was determined using QIIME2 (https://docs.qiime2.org/) [36]. In brief, feature-classifier was performed after the control or denoising of the sequence data with DADA2. The classifier was performed using SILVA_132 as the reference genome database (SILVA_132) under the default condition. Filtering of 10 total counts or less (p-min-frequency 10) and elimination of mitochondria and chloroplasts (p-exclude mitochondria, chloroplast) were also performed. Subsequent analyses were proceeded using the obtained data files. Based on the distribution of OTUs among samples, β-diversity was estimated by non-metric multidimensional scaling (NMDS). The statistical analysis of these NMDS values was performed using the library “vegan” (distance = ”bray”) and “pairwiseAdonis” in the R software package. The plot was visualized by the library “ggplot2”.

Correlation analysis

Relative values of major fecal bacteria (more than 10% of abundances) were visualized by constructing a correlation heatmap via the package “gplots” in the R software (version 4.4.0) after Pearson correlation coefficients were calculated as previously reported [25, 33]. Significant differences in correlation coefficients between the superior and inferior groups were calculated via Fisher's Z transform [37], and p-values were subsequently calculated via the Microsoft Excel function: = 2 × (1 − NORMSDIST(Z)) as previously described [38, 39].

Machine learning methods

This study conducted four types of ML as follows.

Linear discriminant analysis (LDA) is an elementary method of supervised ML. Here, LDA score plots and the cladogram based on LEfSe were visualized by Galaxy (https://huttenhower.sph.harvard.edu/galaxy/) as previously described [25]. In brief, the value provided an LDA cladogram and LDA score plots for the factorial Kruskal‒Wallis test among classes and the value for the pairwise Wilcoxon test between subclasses (set at 0.05) as nonparametric analyses. The pairwise comparisons among subclasses to be performed only among subclasses with the same name were set to “Yes”. The strategy for multiclass analysis was set to “All-against-all (more strict)”.

Association analysis is a basic unsupervised learning technique that is used in market research and environmental analysis. This analysis was applied to classify binarized data beyond the logic of real numbers [40, 41]. Binarized data allow different categories of factor groups (e.g., physical performance and bacterial population) to be treated for the same calculation. Here, data on physical performance and fecal bacterial population were dichotomized on the basis of the median value (M) of the data and sorted as 0 (≤ M) or 1 (> M) as the binarized data for subsequent analysis. We classified causal influences as x and y. The probability (P) was defined as follows: support (x ⇒ y) = P(x ∩ y), confidence (x ⇒ y) = P(x ∩ y)/P(x), and lift (x ⇒ y) = P(x ∩ y)/P(x)P(y). Association rules were determined using reference values for support, confidence, and lift values (support ≥ 0.1, confidence ≥ 0.6, and lift value ≥ 1.2). The R software (https://cran.r-project.org) packages “arules” and “aruleViz” were used. Associative networks were drawn with Force Atlas 2 via Noverlap in Gephi 0.10.1 (https://gephi.org).

Furthermore, feature factors of the bacterial data were extracted via two types of supervised ML: RF, ML with bagging (bootstrap aggregation), and XGB, ML with extreme gradient boosting. Feature factors were selected via the R software library packages “randomForest” and “xgboost” [24, 42]. In this study, the RF parameter was set to “doBest = TRUE” to check the model. The parameter “proximity = TRUE” was used to calculate the values of MeanDecreaseGini, a RF index. The XGB parameters were set to “objective = multi:softmax” and “eval_metric = mlogloss” for multi-class classification to compute Importance, a XGB index.

Causal inference

The direct method for the linear non-Gaussian acyclic model (DirectLiNGAM) [24, 40, 42], which estimates the causal structure through independent component analysis, was used to calculate causal inference. As previously described, the families and genera selected by ML were not separated and used for the calculation even if they had the same attributes, since the calculated results could be evaluated as an internal standard [24, 25]. As a result, the calculation results were confirmed to be close to 1 for data of the same attribute. DirectLiNGAM (https://github.com/cdt15/lingam) was developed via Python (version 3.10.8) on Mac OS Sequoia (version 15.3) as follows: the Python libraries numpy (version 1.26.4), pandas (version 2.2.3), matplotlib (version 3.9.2), scikit-learn (version 1.6.1), lingam (version 1.8.0), and graphviz (version 0.19.1). The calculated data were visualized as directed acyclic graphs (DAGs) via Gephi (version 0.10.1) (https://gephi.org). In addition, the main structure of the complex causal model derived from the DirectLiNGAM approach was statistically evaluated via structural equations. Structural equation modelling (SEM) was performed on a group of components with a large number of directed acyclic graphs (DAGs). For the selected components, SEM was conducted via the library package “lavaan” of the R software as previously reported [24, 40,41,42]. The models used as hypotheses were statistically estimated via the robust maximum likelihood estimator ‘mlr’ of the “lavaan”. The optimal model was assessed by the p-value of Chisq (p > 0.05, nonsignificant), CFI.robust (> 0.95), RMSEA (< 0.05), SRMR (< 0.05), and GFI (> 0.95) as indices of good model fit. The path diagram of the optimal model was visualized via the layout function ‘tree’ of the package “semPlot” in the R software.

Pathway analysis

Metabolic pathway predictions were performed via the PICRUSt2 algorithm (https://github.com/picrust/picrust2) [43]. The relative values of predicted pathways were evaluated via the library packages “limma” and “edgeR” in the R software. Volcano plots were generated via the library package “ggplot2”. The estimated structures of the significant screened pathways (FDR < 0.1) were visualized via DirectLiNGAM as described above.

Statistical analysis

The alpha diversity, relative abundances of dominant and/or characteristic bacteria (> 1% of the detected community for relative abundance analysis), and relative values of predicted pathways were analysed via the R software (version 4.4.0). The Gaussian distribution of the abundance was evaluated via the Shapiro‒Wilk test to select between parametric and nonparametric analyses. Student’s t-test or Welch’s t-test was used as a parametric test following the evaluation of equal variances. The Wilcoxon rank-sum test was used for nonparametric tests. Tukey–Kramer test following ANOVA (Shapiro–Wilk test, p > 0.05) and Dunn’ test following Kruskal–Wallis test (Shapiro–Wilk test, p < 0.05) were used for parametric and non-parametric tests. Significance and tendencies were declared at p < 0.05 and 0.05 ≤ p < 0.10, respectively. All data are presented as the mean ± standard error of the mean.

Availability of data and materials

Raw files of the bacterial V3–V4 16S rRNA sequences were deposited in the DNA Data Bank of Japan (DDBJ) under NCBI BioProject accession number PRJDB18668 (BioProject submission ID: PSUB023775) (BioSample accession number: SAMD00816437-SAMD00816534). The datasets generated and analysed during this study are stored in a supplemental data file (file name: Revised Raw data for Animal Microbiome Revision#2.xlsx) and are available from the corresponding author upon reasonable request.

References

  1. Chasek PS, Wagner LM, Leone F, Lebada A-M, Risse N. Getting to 2030: negotiating the post-2015 sustainable development agenda. Rev European Compar Int Environ Law. 2016;25:5–14.

    Google Scholar 

  2. Nilforooshan MA, Edriss MA. Effect of age at first calving on some productive and longevity traits in Iranian holsteins of the Isfahan Province. J Dairy Sci. 2004;87:2130–5.

    CAS  PubMed  Google Scholar 

  3. Rydhmer L. Genetics of sow reproduction, including puberty, oestrus, pregnancy, farrowing and lactation. Livest Prod Sci. 2000;66:1–12.

    Google Scholar 

  4. Fathoni A, Boonkum W, Chankitisakul V, Duangjinda M. An appropriate genetic approach for improving reproductive traits in crossbred Thai-Holstein cattle under heat stress conditions. Vet Sci. 2022;9:163.

    PubMed  PubMed Central  Google Scholar 

  5. Berry DP, Pryce JE, Wall E. Genetics and genomics of reproductive performance in dairy and beef cattle. Animal. 2014;8:105–21.

    PubMed  Google Scholar 

  6. Silvestris E, de Pergola G, Rosania R, Loverro G. Obesity as disruptor of the female fertility. Reprod Biol Endocrinol. 2018;16:22.

    PubMed  PubMed Central  Google Scholar 

  7. Wathes DC, Fenwick M, Cheng Z, Bourne N, Llewellyn S, Morris DG, et al. Influence of negative energy balance on cyclicity and fertility in the high producing dairy cow. Theriogenology. 2007;68:S232–41.

    CAS  PubMed  Google Scholar 

  8. Ebara F, Inada S, Morikawa M, Asaoka SH, Isozaki Y, Saito A, et al. Effect of nutrient intake on intramuscular glucose metabolism during the early growth stage in cross-bred steers (Japanese Black male × Holstein female). J Anim Physiol Anim Nutr. 2013;97:684–93.

    CAS  Google Scholar 

  9. Matsubara A, Takahashi H, Saito A, Nomura A, Sithyphone K, McMahon CD, et al. Effects of a high milk intake during the pre-weaning period on nutrient metabolism and growth rate in Japanese Black cattle. Anim Sci J. 2016;87:1130–6.

    CAS  PubMed  Google Scholar 

  10. Taguchi Y, Inabu Y, Hayasaki K, Maeda N, Kanmera Y, Yamasaki S, et al. Effects of feeding high volumes of milk replacer on reproductive performance and on concentrations of metabolites and hormones in blood of Japanese black heifer calves. Anim Sci J. 2021;92: e13505.

    CAS  PubMed  Google Scholar 

  11. Robertson SA, Care AS, Moldenhauer LM. Regulatory T cells in embryo implantation and the immune response to pregnancy. J Clin Investig. 2018;128:4224–35.

    PubMed  PubMed Central  Google Scholar 

  12. Esposito G, Irons PC, Webb EC, Chapwanya A. Interactions between negative energy balance, metabolic diseases, uterine health and immune response in transition dairy cows. Anim Reprod Sci. 2014;144:60–71.

    CAS  PubMed  Google Scholar 

  13. Vrieze A, Van Nood E, Holleman F, Salojärvi J, Kootte RS, Bartelsman JFWM, et al. Transfer of intestinal microbiota from lean donors increases insulin sensitivity in individuals with metabolic syndrome. Gastroenterology. 2012;143:913-6.e7.

    CAS  PubMed  Google Scholar 

  14. Zheng D, Liwinski T, Elinav E. Interaction between microbiota and immunity in health and disease. Cell Res. 2020;30:492–506.

    PubMed  PubMed Central  Google Scholar 

  15. Comizzoli P, Power ML, Bornbusch SL, Muletz-Wolz CR. Interactions between reproductive biology and microbiomes in wild animal species. Animal Microbiome. 2021;3:87.

    CAS  PubMed  PubMed Central  Google Scholar 

  16. Silva MSB, Giacobini P. Don’t trust your gut: when gut microbiota disrupt fertility. Cell Metab. 2019;30:616–8.

    CAS  PubMed  Google Scholar 

  17. Chadchan SB, Singh V, Kommagani R. Female reproductive dysfunctions and the gut microbiota. J Mol Endocrinol. 2022;69:R81–94.

    CAS  PubMed  PubMed Central  Google Scholar 

  18. Ashonibare VJ, Akorede BA, Ashonibare PJ, Akhigbe TM, Akhigbe RE. Gut microbiota-gonadal axis: the impact of gut microbiota on reproductive functions. Front Immunol. 2024. https://doiorg.publicaciones.saludcastillayleon.es/10.3389/fimmu.2024.1346035.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Rooks MG, Garrett WS. Gut microbiota, metabolites and host immunity. Nat Rev Immunol. 2016;16:341–52.

    CAS  PubMed  PubMed Central  Google Scholar 

  20. Qi X, Yun C, Pang Y, Qiao J. The impact of the gut microbiota on the reproductive and metabolic endocrine system. Gut Microbes. 2021;13:1894070.

    PubMed  PubMed Central  Google Scholar 

  21. Deng F, McClure M, Rorie R, Wang X, Chai J, Wei X, et al. The vaginal and fecal microbiomes are related to pregnancy status in beef heifers. J Animal Sci Biotechnol. 2019;10:92.

    CAS  Google Scholar 

  22. Yao C, Spurlock DM, Armentano LE, Page CD Jr, VandeHaar MJ, Bickhart DM, Weigel KA. Random forests approach for identifying additive and epistatic single nucleotide polymorphisms associated with residual feed intake in dairy cattle. J Dairy Sci. 2013;96:6716–29.

    CAS  PubMed  Google Scholar 

  23. Taguchi Y, Kurotani A, Yamano H, Miyamoto H, Kato T, Tsuji N, et al. Causal estimation of maternal-offspring gut commensal bacterial associations under livestock grazing management conditions. Comput Struc Biotechnol Rep. 2024;1: 100012.

    Google Scholar 

  24. Kurotani A, Miyamoto H, Kikuchi J. Validation of causal inference data using DirectLiNGAM in an environmental small-scale model and calculation settings. MethodsX. 2024;12: 102528.

    PubMed  Google Scholar 

  25. Okada S, Inabu Y, Miyamoto H, Suzuki K, Kato T, Kurotani A, et al. Estimation of silent phenotypes of calf antibiotic dysbiosis. Sci Rep. 2023;13:6359.

    CAS  PubMed  PubMed Central  Google Scholar 

  26. Handcock RC, Lopez-Villalobos N, McNaughton LR, Back PJ, Edwards GR, Hickson RE. Body weight of dairy heifers is positively associated with reproduction and stayability. J Dairy Sci. 2020;103:4466–74.

    CAS  PubMed  Google Scholar 

  27. Pinart M, Dötsch A, Schlicht K, Laudes M, Bouwman J, Forslund SK, et al. Gut microbiome composition in obese and non-obese persons: a systematic review and meta-analysis. Nutrients. 2022. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/nu14010012.

    Article  PubMed  PubMed Central  Google Scholar 

  28. Tavella T, Rampelli S, Guidarelli G, Bazzocchi A, Gasperini C, Pujos-Guillot E, et al. Elevated gut microbiome abundance of Christensenellaceae, Porphyromonadaceae and Rikenellaceae is associated with reduced visceral adipose tissue and healthier metabolic profile in Italian elderly. Gut microbes. 2021;13:1880221.

    PubMed  PubMed Central  Google Scholar 

  29. Zeng Q, Li D, He Y, Li Y, Yang Z, Zhao X, et al. Discrepant gut microbiota markers for the classification of obesity-related metabolic abnormalities. Sci Rep. 2019;9:13424.

    PubMed  PubMed Central  Google Scholar 

  30. Miyauchi E, Kim S-W, Suda W, Kawasumi M, Onawa S, Taguchi-Atarashi N, et al. Gut microorganisms act together to exacerbate inflammation in spinal cords. Nature. 2020;585:102–6.

    CAS  PubMed  Google Scholar 

  31. Ye J, Zhao Y, Chen X, Zhou H, Yang Y, Zhang X, et al. Pu-erh tea ameliorates obesity and modulates gut microbiota in high fat diet fed mice. Food Res Int. 2021;144: 110360.

    CAS  PubMed  Google Scholar 

  32. Misiak B, Pawlak E, Rembacz K, Kotas M, Żebrowska-Różańska P, Kujawa D, et al. Associations of gut microbiota alterations with clinical, metabolic, and immune-inflammatory characteristics of chronic schizophrenia. J Psychiatr Res. 2024;171:152–60.

    PubMed  Google Scholar 

  33. Inabu Y, Taguchi Y, Miyamoto H, Etoh T, Shiotsuka Y, Fujino R, et al. Development of a novel feeding method for Japanese black calves with thermophile probiotics at postweaning. J Appl Microbiol. 2022;132:3870–82.

    CAS  PubMed  Google Scholar 

  34. Abril-Parreño L, Druart X, Fair S, Krogenaes A. Metabolic signature of cervical mucus in ewe breeds with divergent cervical sperm transport: a focus on metabolites involved in amino acid metabolism. Metabolomics. 2023;19:59.

    PubMed  PubMed Central  Google Scholar 

  35. Chapman AB, Abraham WT, Zamudio S, Coffin C, Merouani A, Young D, et al. Temporal relationships between hormonal and hemodynamic changes in early human pregnancy. Kidney Int. 1998;54:2056–63.

    CAS  PubMed  Google Scholar 

  36. 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:852–7.

    CAS  PubMed  PubMed Central  Google Scholar 

  37. Dai Z, Coker OO, Nakatsu G, Wu WKK, Zhao L, Chen Z, et al. Multi-cohort analysis of colorectal cancer metagenome identified altered bacteria across populations and universal bacterial markers. Microbiome. 2018;6:70.

    PubMed  PubMed Central  Google Scholar 

  38. Eko Wahyunanto P, Heri R, Fitria L, Wisnu Budi W. The quality of Indonesian teachers in the digital era: A meta-analysis. Jurnal Penelitian dan Evaluasi Pendidikan. 2022:186–200.

  39. Rajoria B, Zhang X, Yee D. IGF-1 stimulates glycolytic ATP production in MCF-7L Cells. Int J Mol Sci. 2023. https://doiorg.publicaciones.saludcastillayleon.es/10.3390/ijms241210209.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Miyamoto H, Shigeta K, Suda W, Ichihashi Y, Nihei N, Matsuura M, et al. An agroecological structure model of compost—soil—plant interactions for sustainable organic farming. ISME Commun. 2023;3:28.

    PubMed  PubMed Central  Google Scholar 

  41. Miyamoto H, Kawachi N, Kurotani A, Moriya S, Suda W, Suzuki K, et al. Computational estimation of sediment symbiotic bacterial structures of seagrasses overgrowing downstream of onshore aquaculture. Environ Res. 2023;219: 115130.

    CAS  PubMed  Google Scholar 

  42. Ito K, Miyamoto H, Matsuura M, Ishii C, Nakanishi Y, Suda W, et al. A thermoprotective probiotic function by thermostable lactic acid bacteria and its causal structure. J Functional Foods. 2024;113: 106001.

    CAS  Google Scholar 

  43. Douglas GM, Maffei VJ, Zaneveld JR, Yurgel SN, Brown JR, Taylor CM, et al. PICRUSt2 for prediction of metagenome functions. Nat Biotechnol. 2020;38:685–8.

    CAS  PubMed  PubMed Central  Google Scholar 

Download references

Acknowledgements

We thank Ms. Naoko Tsuji (Sermas Co., Ltd.) for providing technical support for fecal experiments.

Funding

This work was supported by the Japan Society for the Promotion of Science [Grant-in-Aid for Scientific Research no. 16K18814], Mirai Global Farm Co., Ltd., ITOHAM FOODS Inc., NOSAN Corporation, and Sermas Co., Ltd.

Author information

Author notes

  1. Yutaka Taguchi and Haruki Yamano are contributed as Co-first author.

Authors and Affiliations

  1. Kuju Agricultural Research Center, Graduate School of Agriculture, Kyushu University, Taketa, Oita, 878-0201, Japan

    Yutaka Taguchi, Haruki Yamano, Yudai Inabu, Tetsuji Etoh, Yuji Shiotsuka, Ryoichi Fujino & Hideyuki Takahashi

  2. Graduate School of Horticulture, Chiba University, Matsudo, Chiba, 271‑8501, Japan

    Hirokuni Miyamoto, Teruno Nakaguma & Makiko Matsuura

  3. Graduate School of Medical Life Science, Yokohama City University, Tsurumi, Yokohama, Kanagawa, 230-0045, Japan

    Hirokuni Miyamoto

  4. RIKEN Center for Integrative Medical Sciences, Yokohama, Kanagawa, 230-0045, Japan

    Hirokuni Miyamoto, Chitose Ishii & Hiroshi Ohno

  5. Japan Eco-Science (Nikkan Kagaku) Co., Ltd., Chiba, 260-0034, Japan

    Hirokuni Miyamoto & Teruno Nakaguma

  6. Sermas Co., Ltd., Chiba, 271-8501, Japan

    Hirokuni Miyamoto, Teruno Nakaguma, Chitose Ishii & Makiko Matsuura

  7. Mirai Global Farm Co., Ltd, Miyakonojo, Miyazaki, 885-0225, Japan

    Koki Hayasaki, Noriyuki Maeda & Yoshiro Kanmera

  8. ITOHAM FOODS Inc, Nishinomiya, Hyogo, 663-8586, Japan

    Seiji Yamasaki

  9. NOSAN Corporation, Yokohama, Kanagawa, 220-8146, Japan

    Noboru Ota & Kenji Mukawa

  10. Research Center for Agricultural Information Technology, National Agriculture and Food Research Organization, Tsukuba, Ibaraki, 305-0856, Japan

    Atsushi Kurotani

  11. Center for Advanced Photonics, RIKEN, Wako, Saitama, 351-0198, Japan

    Shigeharu Moriya & Satoshi Wada

  12. Keiyogas Energy Solution Co., Ltd., Ishikawa, Chiba, 272-0015, Japan

    Motoaki Udagawa

  13. RIKEN Center for Sustainable Resource Science, Yokohama, Kanagawa, 230-0045, Japan

    Jun Kikuchi

Authors
  1. Yutaka Taguchi
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Haruki Yamano
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Yudai Inabu
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Hirokuni Miyamoto
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Koki Hayasaki
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Noriyuki Maeda
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Yoshiro Kanmera
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Seiji Yamasaki
    View author publications

    You can also search for this author inPubMed Google Scholar

  9. Noboru Ota
    View author publications

    You can also search for this author inPubMed Google Scholar

  10. Kenji Mukawa
    View author publications

    You can also search for this author inPubMed Google Scholar

  11. Atsushi Kurotani
    View author publications

    You can also search for this author inPubMed Google Scholar

  12. Shigeharu Moriya
    View author publications

    You can also search for this author inPubMed Google Scholar

  13. Teruno Nakaguma
    View author publications

    You can also search for this author inPubMed Google Scholar

  14. Chitose Ishii
    View author publications

    You can also search for this author inPubMed Google Scholar

  15. Makiko Matsuura
    View author publications

    You can also search for this author inPubMed Google Scholar

  16. Tetsuji Etoh
    View author publications

    You can also search for this author inPubMed Google Scholar

  17. Yuji Shiotsuka
    View author publications

    You can also search for this author inPubMed Google Scholar

  18. Ryoichi Fujino
    View author publications

    You can also search for this author inPubMed Google Scholar

  19. Motoaki Udagawa
    View author publications

    You can also search for this author inPubMed Google Scholar

  20. Satoshi Wada
    View author publications

    You can also search for this author inPubMed Google Scholar

  21. Jun Kikuchi
    View author publications

    You can also search for this author inPubMed Google Scholar

  22. Hiroshi Ohno
    View author publications

    You can also search for this author inPubMed Google Scholar

  23. Hideyuki Takahashi
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

Hideyuki Takahashi conceived and designed the experiments; Koki Hayasaki, Noriyuki Maeda, Yoshiro Kanmera, Seiji Yamasaki, Noboru Ota, Kenji Mukawa, Tetsuji Etoh, Yuji Shiotsuka, and Ryoichi Fujino conducted the biological experiments; Yutaka Taguchi, Haruki Yamano, Yudai Inabu, Hirokuni Miyamoto, Atsushi Kurotani, Shigeharu Moriya, Teruno Nakaguma, Chitose Ishii, Makiko Matsuura, and Satoshi Wada were involved in data arrangement; Hirokuni Miyamoto, Atsushi Kurotani, Motoaki Udagawa, Jun Kikuchi, Hiroshi Ohno, and Hideyuki Takahashi: prepared to supplementary samples, materials and equipment; Yutaka Taguchi, Haruki Yamano, Yudai Inabu, Hirokuni Miyamoto, Atsushi Kurotani, Teruno Nakaguma, and Hideyuki Takahashi analyzed raw data; and Yutaka Taguchi, Haruki Yamano, Yudai Inabu, Hirokuni Miyamoto, Atsushi Kurotani, and Hideyuki Takahashi were involved in the data analysis and the writing and review of the manuscript and reviewed. All authors have read and approved the final manuscript.

Corresponding authors

Correspondence to Hirokuni Miyamoto or Hideyuki Takahashi.

Ethics declarations

Ethics approval and consent to participate

The procedures used in the present study were performed in accordance with the ARRIVE guidelines and the Guidelines for Animal Experiments by the Faculty of Agriculture at Kyushu University.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher's Note

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

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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 http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Taguchi, Y., Yamano, H., Inabu, Y. et al. Causal estimation of the relationship between reproductive performance and the fecal bacteriome in cattle. anim microbiome 7, 33 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s42523-025-00396-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s42523-025-00396-x

Keywords

Animal Microbiome

ISSN: 2524-4671