Subjects
Healthy children aged 6–9 years were recruited from a kindergarten and primary school in Guangzhou between 2015 and 2017 as reported previously [29, 30]. Two recruitment methods were used: recruitment leaflets distributed in five primary schools in different areas of Guangzhou, and recruitment through parental mutual acquaintances and WeChat public accounts. The following exclusion criteria were adopted: (1) use of antibiotics within 3 months before fecal sample collections; (2) twins and preterm births; (3) incomplete general data or incomplete collection of feces or bone density testing; or (4) history of serious disease. Eventually, 236 participants (145 boys and 91 girls) were included in this study.
This study was approved by the Ethics Committee of the School of Public Health at Sun Yat-sen University (No. 201549). The parents or legal guardians of all participants provided written informed consent prior to enrollment.
Dual X-ray absorptiometry (DXA) measurements for bone mineral content (BMC) and BMD
The weight of all participants was measured to the nearest 0.1 kg using a Tanita MC-780A scale (Tanita Corporation, Tokyo, Japan), with the participants wearing light clothes and no shoes. Height was measured to the nearest 0.1 cm using a portable fixed stadiometer, with the participants in an upright position with no shoes. Whole body DXA scans were performed using the Hologic Discovery W System (Discovery W; Hologic Inc., Waltham, MA, USA) in accordance with the manufacturer’s instructions. The BMC and bone area (BA) values for total body (TB) and total body less head (TBLH) were determined. The BMDs for the TB and TBLH were calculated by dividing the BMC by the BA. The coefficients of variation between two consecutive measurements of BMC (BMD) for TB and TBLH, with repositioning among 33 children on the same day, were 1.09% (1.58%) and 1.37% (2.04%), respectively.
Measurement of gut microbiota
About 90% of the children provided their stool specimens during the process of physical examinations. The specimens were frozen at −80°C within 10 minutes and stored until DNA extraction. The parents or legal guardians of the remaining 10% of the children were asked to collect stool samples at home and place them in a sterile plastic container. The samples were refrigerated at home and transported to the research facility within 12 hours in coolers with ice packs.
The gut microbiota was measured as described previously [29]. Briefly, fecal microbial DNA was extracted using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany). The bacterial genomic DNA was then used as the template for amplifying the V3–V4 hypervariable region of the 16S rRNA gene with the primer pair 341F (Illumina adapter sequence 1 + CCTAYGGGRBGCASCAG) and 806R (Illumina adapter sequence 2 + GGACTACHVGGGTWTCTAAT) by Phusion® High-Fidelity PCR Master Mix with GC Buffer (New England Biolabs). The reaction volume (30 μl) comprised Phusion High-Fidelity PCR Master Mix (New England Biolabs) (15 μl), forward primer (0.2 pmol /μl), 5 μM reverse primer (0.2 pmol /μl), and template DNA (10 ng). Cycling proceeded as follows: 1 min at 98 °C thirty cycles (10 sat 98 °C, 30 sat 50 °C, 60 sat 72 °C); 5 min at 72 °C. The PCR products were purified by electrophoresis on 2% agarose gel in 1× TAE buffer. The purity and concentration of the sample DNA were determined using a Nanodrop 2000 Spectrophotometer (Thermo Fisher Scientific, Wilmington, Delaware, USA).
An Illumina TruSeq DNA PCR-free kit was used to construct the library. After quantification and library testing using Qubit@ 2.0 Fluorometer (Thermo Scientific), the Illumina HiSeq2000 platform was used for high-throughput sequencing. Operational taxonomic units (OTUs) were assigned by clustering the sequences with a threshold of 97% pairwise identity, and chimeras were removed using UPARSE. OTUs were taxonomically assigned at a confidence threshold of 80% based on the Ribosomal Database Project database by the Mothur classifier [31]. All samples were sequenced together and at the same laboratory. Bioinformatics analysis was performed by Novogene Bioinformatics Technology Co., Ltd. (Tianjin, China).
Measurement of SCFAs
Fecal SCFA concentrations were measured as described previously [30] using high-performance liquid chromatography (HPLC) [32]. Approximately 0.3 g of feces were added to 5.0 mL of 70% ethanol and centrifuged at 20°C, 2500 rpm for 10 min. The supernatants were filtered through a 0.45 µm syringe filter. A mixture of 300 µL of the supernatant and 50 µL of 2-ethylbutyric acid (800 mM), an internal standard, was pre-labeled with 2-nitrophenylhydrazide using a Short- and Long-Chain Fatty Acid Analysis Kit (YMC Co., Ltd., Kyoto, Japan). The SCFA derivatives were extracted with n-hexane and diethyl ether and subsequently evaporated to dryness. The obtained fatty acid hydrazide was dissolved in 200 μL of methanol, and 30 μL of the resulting solution was injected into an HPLC system (Agilent 1200, CA, USA) with a YMC-Pack FA column (250 × 6.0 mm; YMC Co., Ltd.). HPLC was performed under the following conditions. The column oven temperature was 50°C; the mobile phase consisted of acetonitrile-methanol-water [30:16:54 v/v/v, pH 4.5 adjusted using 0.1% trifluoroacetic acid (Wako Pure Chemical industries, Japan)]; the flow rate was 1.1 mL/min; and the eluate absorbance was monitored online at a wavelength of 400 nm. A sample of pooled fecal supernatants was analyzed with certain batches of study samples to monitor the analytic precision, and the resulting coefficients of variation were calculated as 2.4%–9.8% and 3.9%–11.4% for intra-day and inter-day reproducibility, respectively, with recovery rates of 80–120% for the SCFAs. A total of 90 samples were tested in triplicate and the intra-class correlation coefficient of fecal SCFAs was from 0.797 to 0.967.
Covariates
Data on household income, parental educational level, and delivery method were collected by trained interviewers through face-to-face interviews. The parental education level was classified into three categories: primary or lower, secondary, and graduate or above. Household income per month was classified into two categories: ≤15,000 yuan and >15,000 yuan. The delivery mode was defined as a binary variable: cesarean or vaginal. Dietary intake for the past year was assessed by a quantitative food frequency questionnaire, and the dietary consumption of energy and fiber was calculated using the 2009 China Food Composition Table [33]. Information on physical activity was obtained using a three-day physical activity questionnaire, and physical activity was calculated by combining the metabolic equivalent score for each type of physical activity after multiplying it by its duration (hours) per day [34].
Bioinformatics analysis
Alpha- and beta-diversity estimates, which indicate within-sample richness and between-sample dissimilarity, respectively, were computed using the R software. The indicators of alpha-diversity included the abundance-based coverage estimator (ACE), Chao1 index, and Shannon index. Principal coordinate analysis (PCoA) was used to depict the beta-diversity at the OTU level and the phylogenetic tree. The Adonis statistical analysis method in permutational multivariate analysis of variance was used to test the significance of the difference in the gut microbial colony structure between the boys and girls.
Statistical analysis
Continuous variables are expressed as means ± standard deviations if normally distributed, or as medians with interquartile ranges if not normally distributed. The Box-Cox method was used to normalize the data. Categorical variables are presented as percentages. Differences in baseline characteristics and the relative abundances of taxa between boys and girls were assessed using Student’s t-test for continuous variables and using the χ2 test for categorical variables.
The BMC is affected by bone size, so we used the residual method to calculate the BMC after size adjustment [35]. Size-adjusted BMC (SA-BMC) for TB and TBLH were derived using the following models [28]:
BMC = β1 × Weight + β2 × Height + β3 × BA + Constant + Residual
The βi, constant, and residual were derived from the linear regression model:
SA-BMC = β1 × Weightmean + β2 × Heightmean + β3 × BAmean + Constant + Residual
Spearman correlation analysis was conducted to evaluate the associations of microbial abundance and SCFAs with DXA-derived bone mass measurements. Multiple linear regression models using the enter method were applied to examine whether the correlations of the gut microbiota abundance and SCFA concentrations with bone mass were independent after adjusting for other potential covariates. In the regression models, we adjusted for weight, age, sex, delivery mode, household income, parental education, physical activity, use of calcium and multivitamin supplements, dietary energy intake, and dietary fiber intake. The false discovery rate (FDR) was used for the P-value correction upon multiple comparisons using the Benjamini-Hochberg method [36]. All analyses were conducted using R version 4.0.0. The significance level was conventionally set at 0.05.