All animal work was conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Agriculture of China. And all methods are reported in accordance with ARRIVE guidelines (https://arriveguidelines.org) for the reporting of animal experiments. The project was also approved by Animal Care and Use Committee (ACUC) in Institute of Animal Sciences, Chinese Academy of Agricultural Sciences (ethical permit: IAS2019-24). Our study included 118 healthy horses. There were including Debao (DeBa) ponies (n = 31), Ningqiang (NiQi) ponies (n = 47) and Guanzhong (GuZh) horses (n = 40) (Supplementary Table S1). All horses were house fed forage and concentrate supplement. The adult animals were selected based on the following criteria: no drugs affecting gastrointestinal microbes were used within 6 months, no reported illness within the past 6 months of the study and, no gut-related disorders recorded until the beginning of the study. All the samples are in the same growing stage. Furthermore, the animals were clinically healthy based on their parasite profiles. One hundred and eighteen fecal samples were collected from the rectum of horses using long arm gloves. Samples were aliquoted into 2 ml cryovials and immediately snap frozen in liquid nitrogen. Then take it to the laboratory and store it in −80°C refrigerator until DNA extraction.

Microbiome DNA was extracted with the Omega Stool DNA kit (Omega, Norcross, GA, USA) follow the appropriate instructions. The DNA concentration and purity were quantified with a Nanodrop 2000^® (ThermoFisher Scientific, USA) and Qubit3.0 (Life Invitrogen, USA), respectively. 1% agarose gel electrophoresis was used to examined DNA quality. According to the literature describe (20, 21), We amplified the v3-v4 (338F-806R) region of the 16S rRNA gene. Brief, the PCR components contained 5 μl of buffer (5×), 0.25 μl of Fast pfu DNA, Polymerase (5 U/μl), 2 μl (2.5 mM) of dNTPs, 1 μl (10 uM) of each Forward and Reverse primer, 1 μl of DNA Template, and 14.75 μl of ddH2O. The PCR amplification program is carried out according to the Hai-Bin Yang program (20), as follows: 98°C for 5 min, followed by 25 cycles consisting of denaturation at 98°C for 30 s, annealing at 53°C for 30 s, and extension at 72°C for 45 s, with a final extension of 5 min at 72°C (20, 21). PCR amplicons were purified with Vazyme VAHTSTM DNA Clean Beads (Vazyme, Nanjing, China) and quantified using the Quant-iT Pico Green dsDNA Assay Kit (Invitrogen, Carlsbad, CA, USA) (20, 21). After the individual quantification step, amplicons were pooled in equal amounts, and pair-end 2×250 bp sequencing was performed using the Illlumina MiSeq platform with MiSeq Reagent Kit v3 at Shanghai Personal Biotechnology Co., Ltd (Shanghai, China).

QIIME2 (Quantitative Insights into Microbial Ecology), version 2019.04 (22) were used to process sequences. Analyze according to the guidance of the official website (https://docs.qiime2.org/2019.4/tutorials/). Removal of primers, quality control, denoising, splicing of double-terminal sequences, removal of chimera and identification of amplicon sequence variants (ASVs) were performed using DADA (23). For taxonomic classification, we selected the Greengenes database (version 13.8) using the Naïve Bayes classifier in QIIME2 (24, 25). ASVs that were identified in only a single sample or classified as non-bacterial were discarded. The sequence of each horse was randomly selected to achieve a uniform sequencing depth for fair comparison (26).

Rarefaction curves were plotted for each sample to determine the abundance of communities and sequencing data of each sample (27, 28). The QIIME2 software has been used to calculate alpha diversity index and evaluated by Observed_species index and Shannon index (29). The Observed_species index is used to evaluate species of abundance, and the higher the value of the index, the more species are included in the sample. The Shannon index measure species diversity, affect the species richness and evenness in the sample microbial community. Differences in alpha diversity indices between groups were analyzed by Kruskal-Wallis Rank sum test, and P < 0.05 was considered significant (30). Beta diversity measurements, including gut microbiota trees, were calculated as described to compare species diversity between different samples (31). Bacterial taxonomic distributions of sample communities were visualized using the R package software. Bacterial biomarkers with markedly different abundance between horses and ponies were analyzed using the LEfSe (linear discriminant analysis effect size) method. LDA (Linear discriminant analysis) was performed on the identified divergent species to estimate the effect size of each divergent species abundance on the difference between groups, with LDA scores >3.5 (32). The relationship between genera and body height was statistically analyzed using Pearson correlation. P < 0.05 was considered statistically significant (32).

Microbial genomic DNA was obtained from 118 samples, and the V3–V4 region of the 16S rRNA was sequenced. The reads of 118 samples were denoised, and a total of 6,657,888 clean reads were obtained. In total, 211, 683 ASVs were obtained by DADA2 software and were annotated 24 phylum and 368 genera. The Shannon-Wiener and grade abundance curves generated by R software are shown in Supplementary Figures S1A–F. This result indicates that the number of ASVs per sample is relatively homogeneous.

The Good's coverage of DeBa ponies, NiQi ponies, and GuZh horses was 0.952 7, 0.942 6, and 0.958 6, respectively, indicating that the proportion of undetected species in the sample is relatively small (Supplementary Figure S2). The Observed species index of DeBa ponies and NiQi ponies was significantly higher than those in GuZh horses (Figure 1A). The Shannon index of DeBa ponies was 9.15, which was significantly lower than that in NiQi ponies (10.16, P < 0.01) and GuZh horses (10.07, P < 0.01) (Figure 1B). The statistical analysis showed that NiQi ponies had a richer and more various gut microbiota. Interestingly, DeBa ponies had more observed species than GuZh ponies, the shannon index value of DeBa ponies was smaller than that of GuZh. The reason for this phenomenon may have something to do with the distribution of horses. DeBa ponies are distributed in the south and GuZh horses in the north (8).

We investigated the relationship between the 118 feces samples from three different regions using Bray–Curtis distances. The UPGMA cluster tree of intestinal microbial structure of three horse breeds was drawn. Each subfield on the tree represented one regions of gut microbiota. Even more interesting is that the gut microbiota of the DeBa pony and NiQi pony clustered together, but those of the GuZh horse located on different subfields (Figure 2A).

We used PCoA (principal coordinate analysis) to examine the gut microbiotas community structures of the three breeds equines. On the PCoA plot (Figure 2B), the bacterial communities from the DeBa pony and NiQi pony clustered tightly and were separated from those from the GuZh horse along principal coordinate axis 1 (PC1), and cluster analysis was similar, which explained the largest amount of variation (16.5%). This result indicating that the composition of intestinal microorganisms in animals of the same body size tends to be similar, which is consistent with the findings of previous studies (33).

We annotate a total of 24 phyla and 368 genera (Supplementary Figure S3). Analysis of the intestinal microbial composition of three species found that the abundance of Firmicutes was the highest, accounting for 28.87–84.98%. The second most abundant phylum was Bacteroidetes with 7.49–66.71% (Figure 3A). At the genus level, Treponema, Oscillospira, BF311 and Ruminococcaceae_Ruminococcus were rich in all samples (Figure 3B). However, bacterial taxa were different between the pony and horse fecal samples. At the phylum level, Firmicutes and Bacteroidetes showed a noteworthy difference in the three groups (Figures 4A–F). In addition, the ratio of F and B (Firmicutes and Bacteroidetes) in the GuZh horses (6.66 ± 1.28) was significantly higher than that in the DeBa (1.04 ± 0.55) and NiQi (1.36 ± 0.41) ponies (P < 0.01) (Table 1). The abundance of Streptococcus was significant higher in the DeBa (12.21%) and NiQi (4.56%) ponies than in the GuZh horses (0.06%). In addition, the abundance of Coprococcus in the GuZh horses (2.16%) than in the DeBa (0.63%) and NiQi (0.53%) ponies (Figures 5A–J). LEfSe analysis found that Phascolarctobacterium, Paludibacter, and Fibrobacter were markedly enriched in NiQi ponies. The relative abundances of Streptococcus and Prevotella were dramatically higher in DeBa ponies than in NiQi ponies and GuZh horses. The relative abundances of Treponema, Treponema Mogibacterium, Adlercreutzia and Blautia were dramatically higher in GuZh horses than in the DeBa ponies and NiQi ponies (Supplementary Figure S4).

The intestinal microbial function of the three horses breeds was predicted using PICRUSt2. PCoA based on the KEGG module revealed differences in microbial function among the DeBa ponies, NiQi ponies, and GuZh horses (Figure 6A). Analogous to the results of PCoA used for assessing beta diversity, the DeBa and NiQi ponies had a similar microbial composition and parallel functions, which were quite different from those of the GuZh horse.

We selected the bacteria with relative abundance more than 1% to analyses their correlation with body height and found that six genera were significantly correlated with body height, including Streptococcus (r = −0.48, P < 0.01), Treponema (r = 0.35, P < 0.01), Coprococcus (r = 0.71, P < 0.01), Phascolarctobacterium (r = 0.31, P < 0.01), Prevotella (r = −0.45, P < 0.01) and Mogibacterium (r = 0.53, P < 0.01). Among the six genera, Coprococcus, Streptococcus, Phascolarctobacterium and Mogibacterium were classified as Firmicutes; Prevotella was classified as Bacteroidetes and Treponema was classified as Spirochaetes. We speculated that Streptococcus, Treponema, Coprococcus, Phascolarctobacterium, Prevotella, and Mogibacterium were the potential microbiota that may affected the body height (Figure 6B).