Estimation of silent phenotypes of calf antibiotic dysbiosis

Reducing antibiotic usage among livestock animals to prevent antimicrobial resistance has become an urgent issue worldwide. This study evaluated the effects of administering chlortetracycline (CTC), a versatile antibacterial agent, on the performance, blood components, fecal microbiota, and organic acid concentrations of calves. Japanese Black calves were fed with milk replacers containing CTC at 10 g/kg (CON group) or 0 g/kg (EXP group). Growth performance was not affected by CTC administration. However, CTC administration altered the correlation between fecal organic acids and bacterial genera. Machine learning (ML) methods such as association analysis, linear discriminant analysis, and energy landscape analysis revealed that CTC administration affected populations of various types of fecal bacteria. Interestingly, the abundance of several methane-producing bacteria at 60 days of age was high in the CON group, and the abundance of Lachnospiraceae, a butyrate-producing bacterium, was high in the EXP group. Furthermore, statistical causal inference based on ML data estimated that CTC treatment affected the entire intestinal environment, potentially suppressing butyrate production, which may be attributed to methanogens in feces. Thus, these observations highlight the multiple harmful impacts of antibiotics on the intestinal health of calves and the potential production of greenhouse gases by calves.


Introduction
The beneficial effects of using antibiotics administered at non-therapeutic concentrations in feed as growth promoters (termed antimicrobial growth promoters, AGPs) were first recognized in the 1940s when chickens fed streptomycin exhibited enhanced growth and feed efficiency 1 .Since then, the growth-promoting effects of several antimicrobial agents have been documented in cattle, poultry, and swine 2 .In recent years, however, antimicrobial growth promoters (AGPs) have been under scrutiny from public health, food safety, and regulatory perspectives due to concerns about the potential for developing antimicrobial resistance (AMR) 3,4 .In addition, resistant bacteria against chlortetracycline (CTC), an antibacterial agent widely used in the livestock industry, have been discovered worldwide 5 .However, in the modern livestock industry, which depends on AGPs, removing AGPs from all feeds on the growth performance, nutrient metabolism, and intestinal environment in calves is relatively unknown.
The gut microbiota is known to have a role in shaping key aspects of postnatal life, such as the development of the immune system 6,7 , and influencing the host's physiology, including energy balance.In addition, antibiotic usage has already been shown to affect fecal microbiota profiles in humans 8 and swine 9 .However, the effects of antibiotics on the gut microbiota, characterized by the use of a cultureindependent metagenomic approach in ruminant species, particularly neonatal calves, have not been examined yet.Therefore, the objective of this study was to evaluate the effects of the lack of AGPs on the growth performance, health condition, and microbial diversity preweaning in beef calves.
Here, the physiological properties and bacterial population of the calves with and without antibiotics treatment were evaluated in detail.These data were evaluated by a procedure jointed with machine learning and statistical causal inference.As a result, statistical causal analyses based on machine learning data estimated the symbiotic bacterial groups that could potentially control the production of butyrate and methane.These observations provide critical viewpoints for developing sustainable livestock technologies to protect animal health and the global environment.

Results and Discussion
To investigate how the lack of AGPs in the preweaning period affects the growth performance, health condition, and microbial diversity, calves were fed milk replacer containing CTC, an antibacterial agent widely used in the livestock industry, at 10 g/kg (CON) or 0 g/kg (EXP) (Fig. 1a).As a result, none of the calves developed major illnesses throughout the study.Administered feed to beef calves with or without AGPs did not adversely affect the phenotype of calves, as indicated by no difference in feed intake and body weight between treatments (Figs.1b and S1ab).On the other hand, the contents of serum components, cytokines, hormones, fecal short-chain fatty acids, and phosphate did not appear to differ significantly between the two groups.(Figs.S2-S4).Differences in bacterial diversity and the major bacterial population were also not always observed between the treatment groups (Figs.S5-S7).However, the correlations between major fecal bacteria, organic acids, and phosphate were clearly different between the two groups (Figs.1c and  S8).In addition, linear discriminant analysis (LDA) effect size (LEfSe) analysis, a supervised machine learning method, revealed that CTC treatment affected the population of various types of fecal bacteria (Figs.1d and S9).In addition, the correlation between the LDA-selected fecal bacterial population and organic acid and phosphate levels was altered by CTC treatment (Fig. S10), although these acid contents did not significantly differ.
Furthermore, association analysis, an unsupervised machine learning method, was performed for the items related to the presence or absence of CTC (Fig. S11).Although phenotypes were not affected by treatment, association analysis showed increased growth performance indicators, including BW, heart girth, waist circumference, and plasma IGF-1 concentration.These results did not contradict previous studies reporting the positive impact of antibiotics on growth in ruminants 10 .Statistical causal inference by LiNGAM (a Linear non-Gaussian Acylic Model) for each of the two groups estimated that the relationship of bacteria was different between the groups with and without CTC (Figs.S12 and  S13).In particular, it suggested that CTC positively affected serum T-Cho concentration and negatively affected fecal butyrate and serum IgA production (Fig. S12a).Moreover, energy landscape analysis showed the entire stability of the physiological components in the experimental duration (Figs.2a) but classified the unstable groups (four groups) after CTC exposure (Fig. S14).In addition, interaction network analyses estimated a negative effect of the genus Methanobrevibacter on the genus Lactobacillus abundance and fecal butyrate concentration (Fig. 2b).LiNGAM inferred a negative causality for butyrate production by antibiotic treatment and the genus Methanobrevibacter (Fig. 2c), together with a positive causality by genus Dorea belonging to family Lachnospiraceae.
As a result of association analysis, increased populations of the genera Methanosphaera and Methanobrevibacter, which belong to the family Methanobacteriaceae, were associated with the CTC treatment (Fig. S11), consistent with the results of LEfSe and energy landscape analysis.The increased population of methanogens is reported to be involved in developing obesity in humans 11 .In humans, methanogenic archaea were reported to be highly resistant to antibiotics 12 , which may explain the increased methanogens in calves administered CTC.Enteric methane emissions are a problem due to energy loss, adverse effects on animal productivity, and environmental issues 11,13 .In addition, methane emissions from beef and dairy cows cause a loss of enteric methane energy, accounting for 2% to 12% of gross energy intake 14 .Thus, lipid accumulation due to an increased methanogen population can explain the relationship between antibiotic administration and growth.
Notably, we previously reported that Caldibacillus hisashii, a member of the class Bacilli, markedly reduces the abundance of the genus Methanobrevibacter 11 .Here, LiNGAM estimated that the class Bacilli works to diminish the effects of antibiotics (Fig. S15) and was causally related to the increase in fecal butyrate level and the abundance of the family Lachnospiraceae, a butyrate-producing bacterial family 15 .In addition, a rise in the fecal family Lachnospiraceae had positive causality with butyrate concentration.Butyrate has various beneficial effects, including improvement of intestinal development and barrier function, mitigation of inflammatory responses, and T cell-independent IgA response 16,17 .Based on these observations, it was inferred that the negative effect of antibiotic administration on fecal butyrate concentration was partly due to an increase in the abundance of methanogens, although the precise mechanism is unknown.
Thus, it is speculated that AGP administration has a negative impact on the gut function and performance of calves via a possible decrease in fecal butyrate production.These observations imply an instability of homeostatic plasticity for growth and environmental load by antibiotics.Overall, the current study first highlighted the possibility that the administration of antibiotics to promote the growth rate decreases intestinal butyrate production through alterations in the population of butyrate-producing bacteria and methanogens, which may have significant implications for the improvement of livestock productivity and reduction in the environmental load by antibiotics and global warming gas.

Animal management
Twelve Japanese black calves (8 male and 4 female calves) were managed with the respective dams until 3 d of age.Thereafter, calves were separated from dams and moved to calf pens at 4 d of age.Milk replacer (Calf Top EX Black, Zen-Raku-Ren, Tokyo, Japan) was offered from 4 d of age using an automated calf feeder (Forster Technique, Co. Ltd., Germany).
Calves were randomly assigned to one of two treatment groups [CON group: n = 6 (4 male and 2 female calves), initial BW = 31.5 ± 5.1 kg; EXP group: n = 6 (4 male and 2 female calves), initial body weight = 30.8± 5.1 kg]: calves in the CON group were fed milk replacer containing CTC at 10 g/kg, whereas those in the EXP group were fed antibiotic-free milk replacer.All milk replacers were formulated for 28.0% crude protein (CP), 18.0% crude fat (CF), and 108.0%total digestible nutrients (TDNs).The amount of milk replacer offered was increased from 0.5 kg/d to 1.0 kg/d from 4 to 21 d of age and maintained at this level until 60 d of age.All calves were fed the same calf starter containing 18.0% CP, 2.0% CF, and 72.5% TDNs (Hello Starter, Zen-Raku-Ren, Tokyo, Japan) and hay containing 12.4% CP, 4.1% CF, and 62.0% TDNs ad libitum from 4 d of age.All diets, except for the CON milk replacer, did not contain antibiotics.All calves had free access to water and mineral blocks (Cowstone A, Nippon Zenyaku kogyo Co. Ltd., Fukushima, Japan) throughout the study.Individual feed intake was recorded daily, and BW was measured every month.The present study was conducted at a Kuju agricultural research center (Oita, Japan), where Japanese black cattle were raised.All experimental protocols were approved by the Kyushu University Laboratory Animal Care and Use Committee (approval no.A30-355-1).The procedures used in the present study were performed according to ARRIVE guidelines and the Guidelines for Animal Experiments by the Faculty of Agriculture at Kyushu University.

Sample collection
On the day before sampling, access to the automated calf feeder was blocked.Blood samples were collected at 3, 30, and 60 d of age immediately before the morning feeding at 0900 h; this was achieved by using Vacutainers to collect serum (Venoject®II, Terumo Co. Ltd., Tokyo, Japan).All samples were centrifuged at 1200 ×g for 30 min at 20 °C, and serum was stored at -80 °C until analysis.Single-use gloves disinfected with 70% ethanol and sampling tubes with spoons were used to collect rectal feces from calves immediately after blood sampling at 3, 30, and 60 d of age.Fecal samples were stored at −30 °C and thereafter stored at −60 °C -−80 °C until analysis.

HPLC analysis for fecal metabolic acids
The fecal samples (200-400 mg) were prepared according to a previous protocol 20 with some modifications.Briefly, the samples were mixed with a 9-fold volume of Milli-Q water for 10 min.After centrifugation at 15,000 rpm, all of the supernatants were filtered with 0.45 μm (Millex-HA Filter Unit SLHA025NB; Merck).The filtered solutions were subjected to highperformance liquid chromatography (HPLC) analysis.To determine the concentrations of lactic acid, acetic acid, propionic acid, butyric acid, valeric acid, isovaleric acid and phosphoric acid, frozen fresh fecal samples were analyzed by using an HPLC Prominence instrument (Organic Acid Analyzer; Shimadzu, Kyoto, Japan), which was equipped on an ion-exclusion column (Shim-pack SCR-102H; Shimadzu) and an electric conductivity detector (CDD-10AVP; Shimadzu).The analytical conditions were as follows: mobile phase, 5 mM p-toluenesulfonic acid; buffer, 5 mM p-toluenesulfonic acid, 20 mM Bis-Tris, and 0.2 mM EDTA-4H; temperature, 40 °C; and flow rate, 0.8 ml/min.

Energy landscape analysis
The energy landscape analysis (ELA) was performed as previously described 23 .Energy landscape analysis is a data-driven method for constructing landscapes that explain the stability of community compositions across environmental gradients.Here, ELA is based on an extended pairwise maximum entropy model that explains the probability for the occurrence of the ecological state of sample  ,  (") given the environmental condition  (") ; as the ecological state, we combined the presence/absence status of selected taxa and levels of physiological factors, ") …  & (") ( where  ) is the number of bacterial taxa and  * is the number of chemicals; as the environmental condition, two environmental factors representing with (1)  or without (0) antibiotic treatment ( +,- (") ) and the growth stages converted to the 0-1 range ( .,- (") ; 3, 30 and 60 days are converted as 0, 0.53 and 1, respectively) were combined as  (") = ( +,- (") ,  .,-") ) .The model can be written as: ) ).
Here, ( (") |) is the probability for the occurrence of an ecological state  (") .Eq. (I) shows that the probability is high when energy % (") /( is low and vice versa.In eq.(II), % (") /( is defined as the sum of the effect of interaction among components, antibiotic treatment, growth stages, and respectively.We used a permutation test (N=2,000) 24 to obtain p-values for each parameter.
All the components in  (") were converted to the 0-1 range as follows.
We first interpret the relative abundance of each bacterial genus in samples as presence (1) or absence (0) states by setting a threshold value of 0.001.
Then, we selected 36 genera that appeared in more than 2 samples but fewer than 35 samples.We combined them into the density of eight selected chemicals, which were converted to the 0-1 range.Accordingly, for each of 36 samples, we obtained the set of explanatory variables  (") with 44 components, which accompanies the environmental condition  (") representing the status for antibiotic treatment and growth stage.

DirectLiNGAM
To estimate a structural model beyond the distribution of limited experimental data 21 , the linear non-Gaussian acyclic model (LiNGAM) approach 25 involves independent component analysis and a non-Gaussian method for estimating causal structures.DirectLiNGAM was established with Python code on the website (https://github.com/cdt15/lingam)(Python version 3.6).The data calculated by the LiNGAM were visualized as LiNGAM networks in Gephi 0.9.2.

Statistical analyses
In addition to LDA, association analysis, energy landscape analysis, and LiNGAM, individual data were analyzed using the MIXED procedure of JMP14 (SAS Institute Inc. Cary, NC, USA) according to the following model: where Yijk is the dependent variable, μ is the overall mean, Gi is the fixed effect of treatment, Tj is the fixed effect of time (age) used as a repeated measure, GTij is the fixed effect of treatment by time after birth interaction, Ck is the random effect of the calf, and eijk is the error term.A simple main effect test was performed to detect the differences between treatment groups at the same point.Significance was declared at P < 0.05, and a tendency was assumed at 0.05 ≤ P < 0. 20.The relative values of dominant and/or characteristic bacteria were visualized through construction of a correlation graph and heatmap after the Pearson correlation coefficient was calculated for the selected bacteria (> 1% of the detected community and > 0.1% of the LDA selected detected community) by using R software (version 4.0.5).The data are presented as the means ± SDs.            ) (X axis) and responsiveness to antibiotics treatment ( !# ) (Y axis) were plotted.Four categories were shown as Group I-IV.The bacteria categorized within Group I have low level of population at 30-60d but increased by antibiotics treatment.The ones within Group II have high level of population at 30-60d and increased by antibiotics treatment.The ones within Group III have low level of population at 30-60d but decreased by antibiotics treatment.The ones within Group IV have high level of population at 30-60d, and/or independent population upon the stage, but the populations decreased by antibiotics treatment.See Fig. 2

Figure and Table legends
0.5).The number of observed OTUs and Chao1, Shannon, and Simpson index values were assessed as measures of α-diversity.The β-diversities were estimated by principal coordinate analysis (PCoA) using weighted or unweighted UniFrac distances based on the OTU distribution across samples.All 16S rRNA gene datasets were deposited in the GenBank Sequencing Read Archive database (accession number: DRA010973; BioSample numbers: SAMD00252058-SAMD00252093 / SSUB016265).LDA Linear discriminant analysis (LDA) is an elementary method of supervised machine learning.Here, LDA effect size (LEfSe) was used to identify genomic taxa characterizing the differences between the experimental conditions.LDA score plots and the cladogram based on LEfSe were visualized by Galaxy (https://huttenhower.sph.harvard.edu/galaxy/),which is described in a previous overview 22 .The populations of predominant bacterial community members (more than 0.1% of the total bacterial population) were analyzed.The threshold of the logarithmic LDA score for discriminative features was set at 3.0.A cladogram based on LEfSe and LDA score plots were made by the value for the factorial Kruskal-Wallis test among classes and the value for the pairwise Wilcoxon test between the net effect of unobserved environmental factors.Parameters in eq.(II), namely,  -1 ,  - + ,  - ., and ℎ -, indicate the effect of the relationship among components (  -1 > 0 favors and  -1 < 0 disfavors the cooccurrence of components  and  ), the effect of the antibiotics on component  (the antibiotic treatment positively ( - + > 0) or negatively ( - + < 0) affects the occurrence of component ), the effect of the growth stages on component  (component  favors the later ( - .> 0) or early ( - .< 0) growth stage) and how likely component  occurs when the other factors are equal,

Fig. 1
Figure and Table legendsFig. 1 (a) The experimental design.(b) Changes in BW during the period.(c) shows the correlation between fecal organic acids and bacterial genera during the period based on the sampling data at 30 d and 60 d.(d) Significant changes in the bacterial population calculated by LDA (p<0.05;> 3-fold change).CON: a group treated with antibiotics (n = 6); EXP: a group treated without antibiotics (n = 6).

Fig. 2
Fig.2 (a) The entire energy landscape associated with antibiotic treatment are shown.The axis formed the energy landscape with compositional energy, community state, and treated time (days).(b) The interaction network shows the significant relationships in the extended pairwise maximum entropy model fitted to the observational data.The blue and red lines show positive and negative effects between the components, respectively.The components are selected by LDA and association analyses, respectively.The bacteria selected by both analyses were underlined.The abbreviations shows as follows: E_: family Erysipelotrichaceae; L_: family Lachnospiraceae; GLU: serum glucose; NEFA: serum nonesterified free fatty acid; Butyrate: fecal butyric acid;

Fig. 1 (
Fig. 1 (a) The experimental design.(b) Changes in BW during the period.(c) shows the correlation between fecal organic acids and bacterial genera during the period based on the sampling data at 30 d and 60 d.(d) Significant changes in the bacterial population calculated by linear discriminant analysis (LDA) (p<0.05;> 3-fold change).CON: a group treated with antibiotics (n = 6); EXP: a group treated without antibiotics (n = 6).
Fig. 2 (a) The entire energy landscape associated with antibiotic treatment are shown.The axis formed the energy landscape with compositional energy, community state, and treated time (days).(b) The interaction network shows the significant relationships in the extended pairwise maximum entropy model fitted to the observational data.The blue and red lines show positive and negative effects between the components, respectively.The components are selected by LDA and association analyses, respectively.The bacteria selected by both analyses were underlined.The abbreviations shows as follows: E_: family Erysipelotrichaceae; L_: family Lachnospiraceae; GLU: serum glucose; NEFA: serum nonesterified free fatty acid; Butyrate: fecal butyric acid; Propionate: fecal propionic acid; T-Cho: serum total cholesterol.(c) The calculated causal relationship of the components strongly linked with butyrate by DirectLiNGAM is visualized.The amounts of changes (Days 3, 30, and 60) with respect to values on Day 3 of the components lined with butyrate (green letters) were used for the calculation.The arrow shows a trend of the causal relationship.The number shows the value of the causal contribution calculated by DirectLiNGAM.The minus and plus value shows a negative and positive causal contribution, respectively.

Figure
Figure S6.(a,b) Relative abundance of fecal microbiota at phylum (a) and genus level (b) for calves fed milk replacer containing Chlortetracycline (CTC) at 10g/kg (CON) or 0g/kg (EXP).(c-e) Relative abundance of genus Bacteroides (c), Prevotella (d), and Prevotellarelated bacteria (e) for the CON and EXP.
Figure.S12 DirectLiNGAM-caluculated results for antibiotics-positive components selected by association analysis are shown.All data (family and genus) in antibiotics group (CON) (a) and in non-antibiotics group (EXP) (b) of Days 3, 30, and 60 (> lift value 1.3 in the association analyses) were used.The paths were visualized by the Gephi.The arrow shows a trend of the causal relationship.The number shows the value of the causal contribution (cost) calculated by DirectLiNGAM.The plus and minus value shows positive and negative causal contribution, respectively.The abbreviation show as follows: g_, genus; f_: family; L_: family Lachnospiraceae; bast, bast_Physique; waist, waist_Physique ; T-Cho, total choresterol; IgA: immunoglobulin A.

Figure S14
Figure.S13 DirectLiNGAM-caluculated results for antibiotics-negative components selected by association analysis are shown.All data (family and genus) in antibiotics group (CON) (a) and in non-antibiotics group (EXP) (b) of Days 3, 30, and 60 (> lift value 1.3 in the association analyses) were used.The paths based on all the data were visualized by the Gephi, respectively.The arrow shows a trend of the causal relationship.The number shows the value of the causal contribution (cost) calculated by DirectLiNGAM.The plus and minus cost shows positive and negative causal contribution, respectively.The abbreviation show as follows: g_, genus; f_: family; E_: family Erysipelotrichaceae.
Figure.S15 DirectLiNGAM-caluculated results of the components associated with class Bacilli (> lift value 1.3 in the association analyses) in group without antibiotics are shown.All data (family and genus) in antibiotics group (CON) (a) and in non-antibiotics group (EXP) (b) of Days 3, 30, and 60 (> lift value 1.3 in the association analyses) were used.The arrows of the right and left side shows the trends of the targeted components in the case with and without antibioitcs, respectively.The upper (red) and lower (blue) arrows show the increased/decreased trend as shown in the association analysis (> lift value 1.3), respectively.The abbreviation show as follows: f_, family, c_ class.