Microbial diversity analysis Both alpha and beta diversity indices were used to estimate the microbiome diversity within—and between microbial communities. Calculations were done with the “summary.single” command in the software package mothur (http://www.mothur.org/; Kozich et al., 2013). Alpha diversity indices analysis included Chao1 index (richness estimate), abundance-based coverage estimator (ACE, richness estimate), Shannon's diversity, and Simpson's diversity index. For the Bray-Curtis similarity, the dataset was rarefied to the minimum number of sequences per sample. Rarefaction curve was constructed based on the observed number of OTUs and nearly reached asymptotes for all samples (data not shown). Principal component analysis was done with JMP® (Version 10.0.0, SAS Institute Inc., Cary, NC). Shared OTUs among gut sites at different age were plotted as Venn diagrams using the R environment (package “VennDiagram,” version 1.6.17.) (Chen, 2016). Heatmaps were created using JColorGrid (Joachimiak et al., 2006).