To determine the changes in tibia-related parameters in broilers in the CON and TD groups, tibia weight, tibia length and tibia diameter were measured (Figures 1A–F). The tibia weight and length of broiler chickens in the TD group were markedly (p < 0.001 and p < 0.001, respectively) decreased compared to the CON group on days 11 and 15. Similarly, the tibia mean diameter of broiler chickens with TD was also significantly (p = 0.001, p < 0.001, and p = 0.005, respectively) decreased during the experiment compared to normal broiler chickens of CON group. Our previous study showed the dynamics of cartilaginous growth plates in broiler chickens play central roles in the proper development and growth of the tibiae (Huang et al., 2018). The morphological observation found that the width of the tibia growth plate (TGP) in the TD group was significantly widened on days 7 and 11 compared to the CON group (Figure 1G). On day 15, the TGP width of the TD group was significantly increased and the area of the TD lesion was increased (Figure 1G). Quantitative histomorphology analysis of the bone showed that the TGP width and its index were markedly higher in the TD group than the CON group on days 7, 11 and 15 (Figures 1E,F). These results collectively demonstrated that tibia-related parameters were obviously altered in TD broilers.

Sankey diagram is a specific type of flow diagram that usually displays the flows and quantities from one set of values to the other proportion, and is used to further determine which experimental stage of TD broiler chickens tibia-related parameters changed most significantly. The results showed that the width of tibia growth plate and its index accounted for a high proportion of tibia-related parameters, mainly derived from the TD group. In addition, the proportion of tibia-related parameters was much more significant on days 10 and 15 than that on day 7 (Figure 1H), indicating that the tibia growth plates of TD broilers showed the most observable changes on days 11 and 15 than on day 7.

The analysis of tibia-related parameters showed significant changes in TD symptoms in broilers on days 11 and 15. Faecal samples are a non-invasive way to screen biomarkers for poultry (Zhu et al., 2021). Therefore, we further screen the potential biomarkers of TD broilers from 11 days of age. 16S rRNA gene sequencing examined the gut microbiota in the feces of both the CON and TD groups. The α-diversity analyzed using different indexes, including the Chao1 index, ACE index, Simpson index, and Shannon index, showed that there was decrease (18.46, 19.72, 2.85, and 22.66%, respectively) in the indexes of gut microbiota listed in the TD group, but the difference was not significant as compared with the CON group (p > 0.05) as shown in Figure 2A. The Venn diagram (Figure 2B) shows that 365 unique OTUs in the TD group compared with the CON group, accounting for approximately 22.9%. To analyze the differences between the CON and TD groups, the bioinformatics analysis methods of unweighted unifrac PCA and OPLS-DA model were performed to evaluate the gut microbiota in feces (Figures 2C,D). The constructed PCA and OPLS-DA score plots showed a clear separation between the CON and TD groups, indicating a significant difference in fecal microbial structure.

Next, we compared the relative abundance features of each group of bacterial taxa to distinguish specific alterations in microbiota. As shown in Figures 3A,B, five phyla with the highest average relative abundance, including the Firmicutes, Proteobacteria, Tenericutes, Actinobacteria and Bacteroidetes, and the proportion of the phyla Proteobacteria (2.82%), Actinobacteria (0.20%) and Bacteroidetes (0.18%) in the CON group was lower than in the TD group (7.09, 0.70, and 0.28%, respectively). Still, the phyla Firmicutes (95.82%) and Tenericutes (0.90%) in the CON group were higher than those in the TD group (91.72 and 0.13%, respectively) (see Supplementary Material S1). In addition, the genera between the CON and TD groups were distinctly enriched in Lactobacillus, Ruminococcus [Ruminococcus], Blautia, Oscillospira, Dorea, Candidatus_Arthromitus, Kurthia, Coprococcus, Ochrobactrum, Coprobacillus, Agrobacterium, Sediminibacterium and Eggerthella (Figure 3B, Supplementary Material S2).

Differential bacterial taxa between the two groups were further identified using the LEfSe analysis (Figures 3C,D), an algorithm for high-dimensional biomarker discovery and interpretation. The LDA was employed to determine the data and effects on the different species. LEfSe analysis identified 50 discriminative features (genus level; LDA score >2, p < 0.05), and the relative abundance was significantly different between the CON and TD groups. Among them, Klebsiella, Arthrobacter, Lvsinibacillus, Citrobacter, Sphingobacterium, Myroides, Oceanobacillus, Wautersiella and Serratia were concentrated in TD broilers; while Ruminococcus, Oscillospira, Dorea, butyrate-producing genera Ochrobacturm, Coprococcus, Brevibacterium, Sediminibacterium, Eggerthella, _Eubacterium, Anaerostipes and Agrobacterium were concentrated in normal broilers (see Supplementary Material S3). These results indicated that the fecal microbial composition of TD broiler chickens has been significantly changed compared to that of normal broiler chickens.

The random forest analysis was performed to screen the top20 important species, showing the highest important species, including the Lactobacillus, Clostridiaceae, Ruminococcaceae, Clostridiales, Lachnospiraceae, RF39 [Ruminococcus], Dorea, Ruminococcus, Cupriavidus, Enterobacteriaceae, Oscillospira, Ochrobactrum [Mogibacteriaceae], Sediminibacterium, Blautia, Methyolbacterium, Clostridium, Acinetobacter and Klebsiella (Figure 4A). In addition, Venn diagram analysis of the LEfSe identified 50 characteristic species (LDA score >2; p < 0.05) and the top 20 species in importance scores by random forest analysis showed that 12 species were screened, and we noted the higher abundances of Clostridiales (p = 0.013), Lachnospiraceae (p = 0.023), Ruminococcaceae (p = 0.030) [Ruminococcus] (p = 0.020), Dorea (p = 0.028), Ruminococcus (p = 0.032), RF39 (p = 0.016), Oscillospira (p = 0.028), Ochrobactrum (p = 0.043), and Sediminibacterium (p < 0.001) in the TD group were clearly lower than those in the CON group, while the abundances of Enterobacteriaceae (p = 0.206) and Klebsiella (p = 0.178) in the TD group were increased compared with the CON group (Figure 4B). Next, the results of Pearson correlation analysis revealed that the tibia growth plate width and its index were correlated with Clostridiales (r = −0.728 and −0.736; p = 0.040 and 0.037, respectively) [Ruminococcus] (r = −0.710 and −0.717; p = 0.048 and 0.045, respectively), Sediminibacterium (r = −0.913 and −0.917; p = 0.002 and 0.001, respectively), Ochrobactrum (r = −0.813 and −0.811; p = 0.014 and 0.015, respectively) as shown in Figure 4C. In addition, tibia-related indicators were positively correlated with Sediminibacterium (r = 0.689, 0.742 and 0.605; p = 0.059, 0.035 and 0.112, respectively), and negatively correlated with Enterobacteriaceae (r = −0.799, −0.634 and −0.743; p = 0.017, 0.091 and 0.035, respectively) and Klebsiella (r = −0.806, −0.461 and −0.846; p = 0.016, 0.250 and 0.008, respectively) (Figure 4C). These results indicated that tibial growth and tibia growth plate quality are closely related to the differential faecal microbiota of TD broiler chickens.

The gut microbiota identified by 16S rRNA gene sequencing can only reach the genus level and has limitations as a potential biomarker, but metabolites derived from differential faecal microbiota become the optative choice (Kim et al., 2020). Next, to explore the composition of fecal metabolites, metabolomic analysis of the faecal samples was performed using UPLC-MS/MS. The results showed that a total of 531 metabolites were identified between the CON and TD groups (see Supplementary Material S4). PCA analysis revealed that the three biological replicates from the two groups clustered together in different regions, indicating significant differences in metabolism between the two groups (see Supplementary Figure S1A). Additionally, the supervised OPLS-DA method was performed to assess the faecal samples. The OPLS-DA score plot showed a clear separation with high reliability using the permutation test (Q2Y = 0.998) between the CON and TD groups (see Supplementary Figure S1B). The differential metabolites of the comparison groups were screened by combining the criteria of fold-change ≥ 2 or ≤0.5 and VIP value ≥1 of the metabolites. The results presented a total of 189 differentially produced metabolites, including 80 downregulated and 109 upregulated, as illustrated in the Volcano plot and Pie chart (Figure 5A, see Supplementary Material S5). Meanwhile, these differential metabolites between the two groups were mapped to the KEGG database for functional clustering analysis, and the results showed that the differentially expressed metabolites were mainly enriched in purine metabolism, vitamin digestion and absorption, glycerophospholipid metabolism, cGMP-PKG signaling pathway, biotin metabolism, etc. (Figure 5B).

The difference analysis was used for comparisons of differential metabolites and the log2 fold change (log2FC ≥ 2.5 or ≤ −2.5) represented the comparison against the reference group based on OPLS-DA (Figure 6A). Compared with the CON group, the differential metabolites with the highest abundance in the TD groups were the thiamine (log2FC = 18.46), N-alpha-acetyl-L-asparagine (log2FC = 17.31), oxaloacetic acid (log2FC = 16.12), indole-3-acetic acid (log2FC = 15.31), oxindole (log2FC = 15.28), guanosine 3′,5′-cyclic monophosphate (log2FC = 14.50), hydrocinnamic acid (log2FC = 14.27), 5-methyluridine (log2FC = 13.84), acetaminophen glucuronide (log2FC = 13.17) and 2′-deoxycytidine-5′-monophosphate (log2FC = 12.27), while the content of urocanic acid, hydroquinone, 2,6-dihydroxypurine and dodecanedioic acid were lower (Figure 6A). These results indicated that the fecal metabolites of broilers with TD were significantly fluctuated. Then, the top 20 important metabolites were screened using random forest analysis, including the dethiobiotin, 4-hydroxybenzaldehyde, anthranilic acid, 1-aminopropan-2-ol, 2, 4-dihydroxybenzoic acid, lysope (16: 0), trimethylamine N-Oxide, lysope (18: 0), glycylphenylalanine, 3-(4-hydroxyphenyl)-propionic acid, 1-methylhistidine, Nα-acetyl-L-glutamine, dulcitol, N-lactoyl-phenylalanine, M-coumaric acid, lysopc (18: 0), D-(+)-malic acid, 2, 6-dihydroxypurine, daltol, DL-o-tyrosine, N-acetyl-L-tyrosine and biopterin (Figure 6B). Moreover, Venn diagram analysis of OPLS-DA (identified 61 characteristic metabolites (log2FC ≥ 2.5 or ≤ -2.5)) and random forest analysis (the top 20 important metabolites) showed that three potential biomarkers were screened, including 4-Hydroxybenzaldehyde, Dethiobiotin, and 2,6-Dihydroxypurine (Figure 6C).

Pearson analysis was performed to unveil the correlations between the abundance of 12 differential fecal microbiota, three differential fecal metabolites, and tibia-related parameters in the host (Figure 7A). The results revealed that dethiobiotin and 2,6-dihydroxypurine were positively correlated with the relative abundances of several important bacterial genera, including the Clostridiales, Lachnospiraceae, Ruminococcaceae [Ruminococcus], Dorea, Ruminococcus, RF39, Oscillospira, and Sediminibacterium. However, 4-hydroxybenzaldehyde was negatively correlated with these important bacterial genus (Figure 7A). More interestingly, the damage indexes (TGP and TGPI) of the tibial growth plate of TD broilers were negatively correlated with the differential metabolites dethiobiotin and 2,6-dihydroxypurine, whereas the damage indexes of the tibial growth plate were positively correlated with 4-hydroxybenzaldehyde (Figure 7A). These results indicated that damage to the tibia growth plate in TD broilers is associated with changes in differential metabolites, which is also closely related to the gut-microbiota in feces.

Receiver operating characteristic (ROC) curve analysis was performed to evaluate the predictive ability and diagnostic performance of differential fecal metabolites in distinguishing TD from normal broiler chickens. Three metabolites with higher importance scores among the differential fecal metabolites, namely 4-hydroxybenzaldehyde, dethiobiotin, and 2,6-dihydroxypurine, were screened and their ROC analysis showed high values of sensitivity, specificity and the area under the curve (AUC) (AUC = 1, p = 0.049; Figure 7B).

Serum ALP and TRACP are important indicators for the assessment of normal bone function. Serum ALP secreted by osteoblast has been extensively studied as a biomarker of bone disease, which is essential for osteoid formation and mineralization. TRACP secreted by osteoclasts is a sensitive index of bone resorption, and its activity correlates with resorptive activity in bone diseases (Li et al., 2010). ROC curves of bone funtion biomarkers are shown in Figure 7C, the AUCs of the serum ALP and TRACP were 0.83 and 0.89, respectively; these serum biomarkers can distinguish the TD from normal broilers, which is lower than the predictive ability of the differentially fecal metabolites. To further validate the predictive ability of these three biomarkers, we sequenced the feces from normal and TD broilers on days 7 and 15 based on metabolomics (Figure 8A). Verification and evaluation of the predictive ability in the three potential biomarkers were assessed based on ROC curve analysis, the results showed that the abundance of 4-hydroxybenzaldehyde, dethiobiotin, and 2,6-dihydroxypurine were significantly different (p < 0.05; Figures 8B–D) and their AUCs were 0.90 (95% CI = 0.76–1.00; p = 0.004), 0.52 (95% CI = 0.23–0.81; p = 0.895), and 0.67 (95% CI = 0.39–0.98; p = 0.233) during the entire trial, respectively (Figures 8E–G). These results indicated that differential fecal metabolites, especially 4-hydroxybenzaldehyde, which is closely related to differential gut-microbiota, could be used as a valuable biomarker for the assessment and diagnosis of broiler chickens with TD.

A total of 180 one-day-old healthy Arbor Acres (AA) broilers were purchased from an intensive commercial broiler farm located in Kaifeng City (Henan, China). They were housed in cages in the livestock laboratory of the Henan Agricultural University (Henan, China). The chickens were reared and performed according to the standard AA broiler management manual as previously described (Huang et al., 2018).

All animal experiments and procedures were performed in strict accordance with the guidance of the Animal Care and Use Committee of Henan Agricultural University and the Care and Use of Laboratory Animals of the College of Veterinary Medicine (approved no: 17–0126), none of the chickens exhibited signs of distress before their death during the experiment. Besides, all the possible measures were carried out to ensure the welfare of the broiler chickens.

All broilers were allowed to adapt to the environment before the further experiment. Broilers were randomly divided into two separate groups: 15 chicks per cage, total four replicates per group. The grouping of the chickens was as: the control group (CON) and the thiram-induced TD chicken group (TD; the addition of 100 mg/kg of thiram from day 4 to day 7, then followed by feeding with a normal diet until the end of the experiment). All experiments lasted for 15 days, with ad libitum provision of feed and water.

Ten broilers were randomly selected from each group for weighing on days 7, 11 and 15 of the experiment, respectively. Blood samples were obtained from the wing vein using heparinized syringes and the serum sample was centrifugated at 3,000 rpm for 10 min at 4°C and stored at −20°C.

Chicks were sacrificed by cervical dislocation in the CON and TD groups. Then, tibia specimens were retrieved for measurement of tibia-related parameters, and the growth plate was stripped from the articular cartilage of the proximal tibia using a surgical knife for morphological analysis of the tibial growth plate as in our previous studies (Huang et al., 2017; Huang et al., 2018). Furthermore, faecal samples were collected from the chickens of both groups on days 7, 11, and 15 to analyze the microbiome and metabolome. The sterile cotton swabs were used to take the upper half of freshly removed chicken feces (within 5 min) in a 2 ml sterile eppendrof (EP) tube, mixed and stored at −80°C. To reduce the difference in sampling time, the first group was sampled in turn until the last group was sampled. We mixed the faeces of 4–6 different broilers into one faecal sample and used four mixed faecal samples per group for subsequent 16S rRNA sequencing and metabolomic analyses to reduce sampling error.

All tibial specimens from chinkens were collected to analyze the weight (g), length (mm), mean diameter (mm) and growth plate width (mm) of the tibia. These bone histomorphometry parameters were determined by an electronic balance sensitive to 0.001 g (#FA1204N; Jinghai Instrument Co., Ltd. Shanghai, China) and Digital Calipers (#SATA91511; TATA Company, Shanghai, China), respectively. In addition, the tibia growth plate index (TGPI, mm/mm) was determined as the width of the tibia growth plate in the broiler. Tibia mean diameter was calculated as the average of the TUD (upper 1/2 diameter of the tibia), TMD (tibia mid-diameter), and TLD (lower 1/2 diameter of the tibia) (Figure 1A).

To determine the levels of calcium (Ca) and phosphorus (P) in serum samples from the CON and TD groups, the electrolyte indicators were quantified by a veterinary reagent tray (Chengdu Puli Taisheng Technology Co., Ltd, Chengdu, China). Quantitative assays were performed with a portable, fully automated veterinary biochemistry analyzer (#SMT-120 V; Seamaty Technology Co., Ltd, Chengdu, China).

Serum concentrations of alkaline phosphatase (ALP) and tartrate-resistant acid phosphatase (TRACP) were measured by an enzyme-linked immunosorbent assay (ELISA) kit (Nanjing Jiancheng Biotechnology Co. Ltd, Nanjing, China) according to the manufacturer’s instructions. The optical density values of each well were determined within 15 min at a wavelength of 450 nm using microplate reader (Model: 30086376 Spark, Tecan Group Ltd, Mannedorf, Switzerland).

According to the manufacturer’s instructions, faecal microbial genome DNA was extracted using a QIAamp DNA tool mini kit (Qiagen, Hilden, Germany). The V3 and V4 hypervariable regions of the bacterial 16S rRNA were PCR-amplified using bacterial primers 338 (forward: 5′-ACT​CCT​ACG​GGA​GGC​AGC​A-3′) and 806 (reverse: 5′-GGACTACHVGGGTWTCTAAT-3′). The amplicons were purified and quantified using Agencourt AMPure Beads and the PicoGreen dsDNA Assay Kit (Invitrogen, Carlsbad, CA, United States), respectively. TheDNA libraries were constructed and then sequencing was performed on an Illumina Miseq 250 platform (Illumina Inc. San Diego, CA, United States) according to the standard protocols by Shanghai Personalbio Technology Co., Ltd, China. Raw Illumina read data for all samples were deposited in the NCBI Sequence Read Archive with the accession code PRJNA755075.

As previously described, the bioinformatics analysis procedure of raw reads was performed (Kong et al., 2020). The bacterial community was compared with their beta diversity using the distance matrices generated from the principal component analysis (PCA) and partial least squares discriminant analysis (PLS-DA). The key bacterial taxa responsible for discrimination between the two groups were identified using linear discriminant analysis (LEfSe) with linear discriminant analysis (LDA) score >2.0 based on the Galaxy online analysis platform (http://huttenhower.sph.harvard.edu/galaxy/). The random forest classifier was performed by 10-fold cross-validation to predict the discrimination between the CON and TD groups in broilers.

Metabolomic analysis was conducted by a commercial service company (Wuhan Metware Biotechnology Co., Ltd, Wuhan, China) following previously reported methods (Huang et al., 2021). Briefly, 50 mg faecal samples were processed with 1,000 μL methanol/water (ice-cold, 70%, v/v) and 5 μl 2-chlorophenylalanine (1 μg/ml). Two pre-cooled sterile steel balls were added to the sample mixture at −20°C for 2 min and homogenized at 30 Hz for 3 min. The sample mixture was centrifuged at 12,000 g for 10 min at 4°C. The supernatant was collected and filtered through a 0.22 μm filter membrane, and quality control (QC) analysis of samples was performed before UPLC-MS/MS analysis (see Supplementary Material S2).

All data analyses were based on the self-built MWDB database (Metware Biotechnology Co., Ltd. Wuhan, China). PCA analysis and orthogonal partial least squares discriminant analysis (OPLS-DA) were performed for the identified metabolites. UPLC-MS/MS analyses were performed, and differentially accumulated metabolites with fold change (FC) ≥ 2 or ≤0.5 and variable importance in project (VIP) ≥ 1 were used as criteria for screening potential biomarkers. In addition, KEGG annotation and enrichment analysis were performed according to our previously described methods [64].

GraphPad Prism version 8 (GraphPad Software, La Jolla, CA, United States) was used to perform statistical analysis and drawing. All data are presented as the means ± SD. Comparisons between groups were conducted using two tailed unpaired student’s t-test. Correlations between the altered gut microbes, significantly different faecal metabolites, and tibia related-indexes (e.g., tibia weight, tibia length, tibia diameter and tibia growth plate width) in the CON and TD groups were analyzed by Pearson correlation test and meanwhile visualized by using GraphPad Prism software. Receiver operating characteristic (ROC) curve analysis obtained (GraphPad Prism V.8) for the display of the constructed models, then the area under the ROC curve (AUC) was used to designate the ROC effect. p < 0.05 was considered to significant.