Comparison of changes in fecal microbiota of calves with and without dam

Animals and experimental design

A total of 16 healthy calves with similar birthday (0 ± 2 d) and body weight (13.1 ± 1.13 kg) (Yak ♂ × Pian cattle ♀) were selected and randomly assigned to one of two groups according to their cultivation method. The control group was cultivated with heifers, whereas the other group was fed milk replacer. Milk replacer and warm boiled water (50–60 °C) in a ratio of 1:8. The milk replacer was washed and thoroughly stirred, then cold boiled water added and the water temperature adjusted to 37 °C. The daily feeding amount of milk replacer is 3% of body weight, and adjusted at any time with the increase of body weight. After 10 days of age, the control group followed the cows to pick small amounts of the TMR feed, and the treatment group was provided free access to the same diet. The diet was formulated according to China Feeding Standard of Beef Cattle (NY/T 815-2004) and the animal was fed twice daily at 07:00 and 17:00. After 3 days of adaptation period, the formal experiment lasted 90 days. On day 35, 65 and 95, the feces of calves were collected in 50 mL test tubes with lids. The samples were frozen immediately at −80 °C until use.

DNA extraction

A total genomic DNA was extracted from fecal samples, according to the manufacturer’s instructions of DNeasy PowerSoil Kit (Qiagen, Valencia, CA, USA). After extraction, DNA integrity was detected by 0.8% agglutinate gel electrophoresis, and DNA concentration and purity were detected by NanoDrop ND-1000 spectrophotometer (Nyxor Pharmacia, Paris, France). All extracted DNA samples were stored at –20 °C until be used as templates for real-time PCR and Illumina sequence analyses.

PCR amplification and 16s rRNA gene sequencing

A template, fecal total DNA, was amplified by the V4 variable of 16s rRNA genes. The universal primer set, 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) are used (Caporaso et al., 2011). Specific primers with Barcode were used for PCR according to the selection of sequencing region. Each 25 μL system consisted of a 1 × PCR buffer, 1.5 mM MgCl2, each deoxynucleoside triphosphate at 0.4 μM, each primer at 1.0 μM, 0.5 U of KOD-Plus-Neo enzyme (Toyobo, Tokyo, Japan) and 10 ng of template DNA. And the PCR amplification procedure consisted of starting of denaturation at 94 °C for 1 min, then 30 cycles (including denaturation at 94 °C for 20 s, annealing at 54 °C for 30 s and elongation at 72 °C for 30 s) and a final extension 72 °C for 5 min. The PCR reactions were repeated three times for each sample and combined together. The PCR procedure was mixed with 1/6 volume of 6X loading buffer and detected by 2% agarose gel electrophoresis. We chose the sample between 200~400 bp, then the PCR production was purified using OMEGA Gel Extraction Kit (Omega Bio-Tek, Norcross, GA, USA), the gel purified barcoded amplicons were pooled with equal molar amount and quantified on a Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA). The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific, Waltham, MA, USA) and Agilent Bioanalyzer 2100 system. Then, the pooled amplicons were paired-end sequenced (2 × 250 bp) on the Hiseq Illumina Sequencing Platform (Rhonin Biosciences Co., Ltd, Chengdu, China).

Bioinformatics and statistical analysis

Paired-end reads from the original DNA fragments were merged using FLASH (Magoc & Salzberg, 2011). Paired-end reads was assigned to each sample according the unique barcode. The sequence with quality (length > 200 bp, without ambiguous base ‘N’, and average base quality score >30) were screened for chimmeras checking using Uchime (Edgar et al., 2011). Sequence were clustered into OTUs at 97% identity threshold using UPARSE algorithm (Edgar, 2013), and picked representative sequences for each OTU. Taxonomy was assigned using Silva database (Edgar, 2013) and the representative sequences were aligned using PyNAST (Caporaso et al., 2010). The analysis of alpha, which included calcuation of Chao1, ACE, Shannon and Simpson indices were conducted with Vegan (Version 2.0-2.R CRAN packet) and PD index was conducted with Picante (Kembel et al., 2010). Rarefaction curves were generated based on these metrics. Principal Co-ordinates analysis (PCoA) was performed to assess significant different between samples. Linear discriminant analysis (LDA) effect size (LEdSe) method was performed to identify the bacterial taxa differently represent between groups at genus or higher taxonomy level, which would help discover biomarkers (Segata et al., 2011). All statistical analysis was run using SAS v9.4 (SAS Institute Inc, Cary, NC, USA). The level of statistical significance was set at P < 0.05.

Analysis of Illumina MiSeq sequencing data in feces of calves with and without dam

A rarefaction test was performed at the OTU level and the results are presented (Fig. 1). All the curves asymptotically approached a plateau, which suggested that new phylotypes could not be detected even if additional sequencing was performed, and the majority of the bacterial phylotypes in feces were identified.

Illumina MiSeq sequence of the 16S rRNA yielded 2,624,000 valid sequences from the 48 samples (including eight calves with dam and eight calves without dam, and samples were collected on day 35, 65 and 95 respectively), with a mean of 54,667 sequences per sample. These sequences were clustered into 62,010 operational taxonomic units (OTUs) at the 97% similarity level, with an average of 1,292 OTUs per sample.

According to the Venn diagram at a 97% similarity level (Fig. 2), the shared OTUs between the calves with and without dam was 1,333 (85.41% of sequence, Fig. 2A), the OTUs unique to the calves without dam was 1,165 (11.12% of sequence, Fig. 2A) on day 35. On day 65, the shared OTUs between calves with and without dam was 1,522 (90.57% of sequence, Fig. 2B), the OTUs unique to the calves without dam was 774 (5.92% of sequence, Fig. 2B). On day 95, the shared OTUs between calves with and without dam was 1,537 (86.01% of sequence, Fig. 2C), the OTUs unique to the calves without dam was 825 (10.06% of sequence, Fig. 2C). A total of 1,848 and 1,794 OTUs (99.16% and 98.78% of the total sequences) were common to calves with and without dam in three different timepoints, respectively. In calves with dam, the OTUs unique to day 35, 65 and 95 were 55 (0.04% of sequence), 47 (0.04% of sequence) and 52 (0.04% of sequence), respectively (Fig. 2D). In the meantime, the OTUs unique to day 35, 65 and 95 were 417 (0.2% of sequence), 154 (0.07% of sequence) and 183 (0.1% of sequence) in calves without dam (Fig. 2E).

The OTUs, species richness and diversity were estimated using the Chao1, ACE, Shannon and Simpson indices are shown in Table 1. On day 35, The PD indices of calves without dam was higher than calves with dam (P < 0.05). No significant difference was observed in the Chao1, ACE, Simpson and Shannon indices between the two groups (P > 0.05). On day 65, all the α-diversity indices including the Chao1, ACE, Simpson, Shannon and PD indices of calves without dam were lower than calves with dam (P < 0.01). However, on day 95, only the Chao1 and ACE indices of calves without dam were lower than calves with dam (P < 0.01).

In the time dimension, the α-diversity indices including the Chao1, ACE Simpson and PD indices of day 65 were greater than day 35 (P < 0.01), but no significant difference with day 95 in the calves with dam (P > 0.05). On the contrast, no significant difference on the Chao1, ACE and PD indices among the three timepoints was observed in the calves without dam (P > 0.05); however, the Shannon indices of day 65 and 95 were higher than that of day 35. These results indicated that the bacterial species richness and diversity in feces of calves with dam was greater than calves without dam in different timepoints, and the bacteria in the gastrointestinal of calves with dam is established earlier than calves without dam.

Clustering of the samples based on unweighted analysis of OTUs at 97% identity is shown in Fig. 3. The closer the samples and the shorter the branches, indicating more similar the species composition of the two samples. Surprisingly, the samples were absolutely clustered into two different branches on day 35, 65 and 95 between calves with and without dam (Figs. 3A–3C), which means that the microbial composition was totally different in feces on the three different timepoints between the two groups.

The relationship of bacterial community among the three different timepoints showed that the samples were clustered into three branches (Fig. 3D). All eight samples of day 65 were clustered into the same clade, while two samples of day 35 (sample nos. 1-4 and 1-6) were clustered together with a sample of day 95 (sample nos. 3-8) into a seperate clade and then merged with samples of day 65 and 95 at the distance of less than 0.4 in calves with dam. This means that the microbial composition in calves with dam changed from day 35 on and kept stable thereafter. On the contrary, a totally different alteration trend was observed in the composition of calves without dam (Fig. 3E). All eight samples of day 95 and two samples of day 65 (sample nos. 2-1 and 2-2) were clustered into a clade at the distance of 0.4, while the samples of day 35 and 65 were clustered into one clade first and then mingled together at the distance over 0.4. Two samples of day 35 (sample nos. 1-1 and 1-2) were clustered into one clade and mingled with the left samples at the distance of over 0.5. This means that the microbial composition in calves without dam changed from day 65 on and kept stable thereafter. The evenness of microbial composition of calves with dam is superior than that of calves without dam.

PCoA was conducted based on unweighted Unifrac distance to distinguish fecal microbial structure between the calves with and without dam on day 35, 65 and 95 respectively and the results are shown Fig. 4. The closer the distance between points the more similar the microbial community structure of the sample. The PCoA plot demonstrated that the microbial community between the calves with and without dam was different, and could be distinguished clearly (Figs. 4A–4C), and PCo1 accounted for 51%, 56.7% and 55% of variance on day 35, 65 and 95, respectively.

The characteristics of the PCoA of the three timepoints was that all 16 samples on day 65 and 95 were closely clustered together, while six samples on day 35 were vertically isolated from samples of day 65 and 95 in calves with dam (Fig. 4D), and the PC1 account for 16% of the variance. In terms of the samples of calves without dam, the samples of day 35 and 65 were horizontally distinguished from samples of day 95, but two exceptional plots of day 35 and 65 respectively existed (Fig. 4E), and the PC1 account for 12.9% of the variance. This again suggested that the establishment process of microbiota in calves with dam was earlier that calves without dam.

Bacterial composition in feces between calves with and without dam

Species annotation and statistics were performed on different taxonomic levels using the effective sequences obtained Fig. 5. At the phylum level, the most common phyla in feces were Firmicutes, Bacteroidetes, Proteobacteria and Actinobacteria in calves with and without dam. Compared with calves with dam, only the abundance of Actinobacteria was higher than in calves without dam in the three different timepoints (P < 0.05). As time passed by, the abundance of Firmicutes in calves with dam increased, while Proteobacteria and Actinobacteria decreased (P < 0.05). The abundance of Bacteroidetes and Proteobacteria in calves without dam increased from day 35 to 65, and then decreased from day 65 to 95 (P < 0.05).

At a lower taxonomical level, the top 10 were shown in Fig. 6. On day 35, the dominating genera (relative abundance ≥4%) in calves with dam were Bacteroides (32.91%), Escherichia-Shigella (9.79%), Lactobacillus (8.45%) and Parabacteroides (4.75%) while the genera Peptostreptococcus (21.05%), Lactobacillus (7.61%), Myroides (6.05%), Bacteroides (4.38%) and Acinetobacter (4.23%) were dominating in calves without dam. Compared with calves with dam, the abundance of Bacteroides, Escherichia-Shigella and Parabacteroides were lower and the abundance of Peptostreptococcus, Myroides and Paraeggerthella were higher in calves without dam (P < 0.05).

On day 65, the dominating genera (relative abundance ≥4%) in calves with dam were Bacteroides (7.96%), Rikenellaceae RC9 (4.03%), Treponema 2 (4.55%), Ruminococcaceae UCG-010 (4.64%) and Ruminococcaceae UCG-005 (11.33%), while the genera Bacteroides (15.59%), Ruminococcaceae UCG-005 (4.62%), [Eubacterium] coprostanoligenes group (5.47%) and Ruminococcaceae UCG-005 (5.47%) was dominating in calves without dam. Compared with calves with dam, the abundance of Bacteroides, Lactobacillus, and Ruminococcaceae UCG-002 were higher (P < 0.05), while the genera Ruminococcaceae UGG-005, Rikenellaceae RC9 gut group, Treponema 2 and Ruminococcaceae UCG-010 were lower in calves without dam (P < 0.05).

On day 95, the dominating genera (relative abundance ≥4%) in calves with dam were Rikenellaceae RC9 (6.92%), Prevotellaceae (4.03%) and Ruminococcaceae UCG-005 (17.8%), while the genera Bacteroides (4.33%), Ruminococcaceae UCG-005 (7.12%), [Eubacterium] coprostanoligenes group (5.54%), Lachnospiraceae NK4A136 group (5.54%) and Treponema 2 (5.97%) were dominating in calves without dam. Compared with calves with dam, the genera of [Eubacterium] coprostanoligenes group and Treponema 2 were higher (P < 0.05) in calves without dam, while the genera Ruminococcaceae UCG-005 and Rikenellaceae RC9 gut group were lower (P < 0.05).