Characteristics of Gut Microbiome and Its Metabolites, Short-Chain Fatty Acids, in Children With Idiopathic Short Stature

ISS patients who visited the pediatric outpatient department of Ruijin Hospital affiliated with Shanghai Jiao Tong University School of Medicine between October 2018 and February 2020 were enrolled with written consent. Patients were recruited based on the following criteria: (1) the height of pre-adolescent children aged 4-8 years was two standard deviations (SD) lower than the normal reference value of the same region, sex, age, and race; (2) normal birth weight and body length, annual growth rate of height < 5 cm, and delayed or normal bone age; (3) the peak values of growth hormone (GH) stimulation tests were all greater than 10 µg/L; and (4) no other malformations or organic diseases (chronic diseases of the bone, liver, kidney, heart, and lung), no serious psychological and emotional disorders. The control group consisted of healthy children from East China. Inclusion criteria: 4-8 years old, and height within the range of plus or minus one standard deviation (-1SD to +1SD) of the mean normal reference values of the same sex, age, region, and race. The exclusion criteria were as follows: (1) family history of familial short stature and constitutively delayed pubertal development or any other genetic disease; (2) suffering from congenital malformation, mental retardation, special facial features, obesity, or digestive tract diseases; (3) received antibiotic and probiotic intervention in the last three months. All the patients and healthy controls were from east China (with similar dietary habits) to minimize the impact of diet differences on gut microbiota.

This study was approved by the Ethics Review Committee of Ruijin Hospital affiliated to Shanghai Jiao Tong University School of Medicine and was conducted in accordance with the principles of the Helsinki Declaration. All participants and their guardians provided informed consent and signed informed consent forms (Project identification code: KY2020-121).

Stool samplers and sterile stool sampling tubes (Sarstedt, Germany) were used to collect fresh stool samples from the subjects. Each tube of stool sample weighed 5–10 g and was stored in a refrigerator at -80°C. Four to five milliliters of elbow venous blood was collected from the subjects, into a pro-coagulant vacuum vein collection, and placed in a 4°C refrigerator for 2 h. The blood samples were centrifuged at 1300g/min for 15 min. The serum was saved to sterile tubes and stored at -80°C until further use.

Serum IGF-1 concentration was detected using an IMMULITE2000 IGF-1 kit (L2KGF2, Siemens, UK), following the manufacturer’s instructions (19). IGF-I standard deviation scores were calculated using the laboratory reference data for a survey of normal children (20). The CNP-22 ELISA kit (Cat. NO: EKE-012-03 Phoenix Pharmaceuticals Inc., USA) was used to detect the serum NT-proBNP concentration, following the manufacturer’s instructions.

Fecal bacterial DNA was extracted using TIA Namp Stool DNA Kit (Invitrogen, California, USA) according to the manufacturer’s protocol, DNA concentration and integrity were determined using a Nanodrop ND 1,000 Spectrophotometer (Thermo Fisher Scientific, United States) and the quality of the extracted DNA was evaluated by 1.2% agarose gel electrophoresis, respectively. PCR amplification of the bacterial 16S rRNA genes V3–V4 region was performed using the forward primer 338F (5’- ACTCCTACGGGAGGCAGCA -3’) and the reverse primer 806R (5’- GGACTACHVGGGTWTCTAAT -3’). Thermal cycling consisted of initial denaturation at 98°C for 2 min, followed by 25 cycles consisting of denaturation at 98°C for 15 s, annealing at 53°C for 30 s, and extension at 72°C for 30 s, with a final extension of 5 min at 72°C. PCR products were purified with Vazyme VAHTS™ DNA Clean Beads (Vazyme, Shanghai, China) and quantified using a Quant-iT™ PicoGreen dsDNA Assay Kit (Invitrogen, United States). Sequencing libraries were prepared using Illumina’s TruSeq Nano DNA LT Library Prep Kit (Illumina, Delaware, USA), and sequenced on an Illumina Hiseq2500 platform. The sequencing service was provided by Personal Biotechnology Co., Ltd (Shanghai, China).

Sequencing data analyses were mainly performed using QIIME2 (https://docs.qiime2.org/2019.4/citation/↗) (21) and R packages (v3.6.0). The raw sequence data were demultiplexed and quality filtered using the q2‐demux plugin followed by denoising with DADA2 (22). All amplicon sequence variants (ASVs) were aligned with mafft (23) (via q2‐alignment) and used to construct a phylogeny with fasttree2 (24) (via q2‐phylogeny). ASV-level alpha diversity indices, such as Chao1 richness estimator, were calculated using the ASV table in QIIME2, by Kruskal-Wallis test, and visualized as box plots. Wilcoxon rank-sum test was utilized to detect differentially enriched microbiota at genus between ISS and healthy controls groups (P < 0.05). Beta diversity analysis was performed to investigate the structural variation of microbial communities among different samples using Bray-Curtis metrics and visualized via principal coordinate analysis (PCoA) (25). Non-parametric permutational multivariate analysis of variance (PERMANOVA) was conducted for analyzing significant difference of microbiota structure between the two groups. LEfSe (Linear discriminant analysis effect size) was performed to detect differentially abundant taxa in both groups (26). The raw data were deposited into the NCBI Sequence Read Archive (SRA) database (accession number PRJNA819198, https://www.ncbi.nlm.nih.gov/Traces/study/?acc=PRJNA819198↗).

With appropriate dose of acetic acid, propionic acid, butyric acid, isobutyric acid, valeric acid, isovaleric acid(Sigma, Germany) and caproic acid(Aladdin, The United States), diethyl ether (Sinopagic Chemical Reagent Company, China) was added to make 10 mixed-standard concentration gradients(0.02 μg/mL, 0.1 μg/mL,0.5 μg/mL, 2 μg/mL,10 μg/mL,25 μg/mL,50 μg/mL,100 μg/mL,250 μg/mL,500 μg/mL). Feces were thawed on ice, and approximately 50 mg of feces were homogenized in 50 μl 15% phosphoric acid. The suspensions were homogenized with a vortex and centrifuged at 15285g for 10min at 4°C. The supernatants were processed for gas chromatography-mass spectrometry on a Thermo TRACE 1310-ISQ gas mass spectrometer (Thermo, Massachusetts, USA). The injection volume was 1 μl, and the split ratio was 10:1. Samples were separated with an Agilent HP-INNOWAX capillary GC column (30 m × 0.25 mm ID × 0.25 µm). The initial temperature was 90°C and was increased to 120°C at 10°C/min, after which the temperature was increased to 150°C at 5°C/min and then to 250°C at 25°C/min, where it remained for 2 min. The carrier gas was helium (1.0 ml/min). The temperatures of the injection port and ion source were 250°C and 230°C respectively under SIM model. Agilent MSD ChemStation software was used for quantitative analysis of chromatographic peak area and retention time (27).

Data documentation and statistical analyses were performed using SPSS software version 23.0. Parametric data are expressed as the mean and standard deviation (SD). SD scores were calculated for height (Ht SDS) and IGF-1 (IGF-1 SDS). The differences between the groups were examined using the independent sample t-test. Continuous nonparametric data are represented as median values (interquartile spacing [IQR]) and were analyzed using the Wilcoxon rank-sum test. Receiver operating characteristic (ROC) analysis was performed to assess the diagnostic value of the identified indicators by using pROC package and ggplot2 package of R language. Spearman rank coefficient was used to test the correlation between the two study variables. The correlation coefficient (r) > 0.3 and < 0.5 was considered to be moderate correlation. The correlation coefficient (r) > 0.5 was considered to be strongly correlated, and P < 0.05 was considered to be statistically significant.

Thirty-two subjects aged 4-8 years were enrolled in the study. Sixteen patients with ISS (median age 6.4 years, 44% males) were enrolled as the ISS group and 16 healthy children of normal height (median age 6.6 years, 50% males) were enrolled as the control group. The height standard deviation score (HtSDS) of the ISS group was -2.33 ± 0.29, while that of the control group was 0.52 ± 0.27, with a statistical difference between the two groups (P < 0.001) (Figure 1A). The body mass index standard deviation score (BMI SDS) of the ISS group was 0.04 ± 0.74, and the BMI SDS of the control group was 0.69 ± 0.56. There was a statistically significant difference between the two groups (P = 0.009) (Figure 1B). Age, sex, and target height (TH) did not significantly differ between the ISS and control groups (P > 0.05) (Figures 1C, D). The GH peak of the growth hormone stimulation test in the ISS group was 12.77 ± 1.9 ng/mL.

As mentioned, IGF-1 and amino-terminal propeptide of C-type natriuretic peptide (NTproCNP) deficiencies are known factors contributing to the short stature. To evaluate the growth of children with ISS and healthy controls, we measured the levels of IGF-1 and NTproCNP, which are important regulators of serum longitudinal bone growth in childhood. The results showed that the serum NTproCNP of ISS group was 0.46 ± 0.08 ng/mL, and that of the control group was 0.60 ± 0.04 ng/mL. There was a statistically significant difference between the two groups (P < 0.001) (Figure 1E). The integral standard deviation of IGF-1 (IGF-1 SDS) of ISS group was 0.06 ± 0.77, and that of control group was 1.03 ± 0.82, indicating a statistically significant difference between the two groups (P = 0.002) (Figure 1F).

We next evaluated the gut microbiota in ISS, as evidence from animal studies suggests that gut microbiota can regulate host growth (15, 16, 18). We collected stool samples from our ISS patients and control cohort for 16S rRNA profiling. The alpha diversity index was used to analyze the richness and diversity of gut microbiota. The analysis results showed that Chao1 index and Observed species index in ISS group were significantly lower than those in control group (P = 0.018, P = 0.038), suggesting that the species richness of ISS group was lower than that of control group (Figures 2A, B). Simpson and Shannon indices, reflecting the diversity between the two groups, showed no significant differences (Figure 2).

Beta diversity analysis showed that the gut microbiota of the ISS and control groups were significantly separated, and the dispersion degree between samples was large; however, gut microbiota among samples from the same group appeared to be relatively concentrated. Adonis analysis was used to verify the difference between the two groups. The results confirmed that there were significant differences in gut microbiota between the two groups (Adonis test, Bray-Curtis distance, R2 = 0.092, P < 0.001; unweighted UniFrac distance, R2 = 0.069, P < 0.001), indicating significant changes in the composition of gut microbiota in ISS patients (Figures 2E, F). The main strains at the phylum and genus levels in the ISS and control groups are shown in Figures 2G, H. The main strains at the phylum level were Firmicutes, Bacteroidetes, Actinobacteria, Proteobacteria, and Verrucomicrobia. The main strains at the genus level were Bacteroides, Faecalibacterium, Bifidobacterium, Clostridium, and Blautia.

Linear discriminant analysis effect size (LEfSe) was used to analyze the difference in abundance of species between the two groups and to determine which species contributed to the significant difference between groups. Five different species were screened at the family level, as shown in Figure 3A. Compared to the control group, the relative abundance of Ruminococcaceae and Erysipelotrichaceae in the ISS group was significantly decreased, while the relative abundance of Enterococcaceae, Turicibacteraceae, and Eubacteriaceae was significantly increased in ISS group. Thirteen different species were found significantly different at the genus level. Compared to the control group, the relative abundance of Faecalibacterium, Eubacterium, Gemmiger, Collinsella, Leuconostoc, Slackia, and Macroglobules was significantly decreased, while the relative abundance of Parabacteroides, Clostridium, Anaerofustis, Turicibacter, Eggerthella, and Anaerotruncus was significantly increased in the ISS group. Figure 3B shows the results of LDA Effect Size of species in the ISS group and control group at the family and genus levels. Figure 3C shows the increased genus of bacteria in ISS and Figure 3D the decreased bacteria in ISS. Of these, The genus Faecalibacterium is one of the most abundant and important symbiotic bacteria in the human gut microbiota, accounting for approximately 5% of the total bacteria in feces and belongs to the family Ruminococcaceae (28, 29). The genus Eubacterium is also one of the core genera of the human gut microbiota, which is widely found in the human intestinal tract and has high specificity and adaptability to the human intestinal tract (30, 31).

We measured if ISS is associated with altered SCFA levels in our patient cohort. Quantitative analysis of SCFAs concentrations in fecal samples by GC-MS showed significant differences between the ISS and control group. The total SCFAs concentration of ISS group was 3450.33 ± 963.94 µg/g; butyrate: 698.35 ± 261.96 µg/g; valerate: 91.98 ± 105.79 µg/g. The total SCFA concentration in the control group was 5200.43 ± 1440.66 µg/g; butyrate, 1534.99 ± 563.71 µg/g; valerate: 275.80 ± 186.08 µg/g. The concentrations of total SCFAs and individual components butyrate and valerate in the ISS group were significantly lower than those in the control group (P < 0.05) (Figure 4). Spearman correlation analysis was carried out between SCFAs and serum bone longitudinal growth regulatory factors, and the results showed that total SCFA and butyrate levels were significantly positively correlated with IGF-1 SDS levels (R > 0.4, P < 0.05) (Figure 5B).

Spearman correlation analysis was conducted between different bacterial species and clinical indicators to determine the correlation between gut microbiota composition and host health status (Figure 5A). The family Ruminococcaceae and genera Faecalibacterium and Eubacterium were significantly positively correlated with HtSDS (r = 0.502, P = 0.005; r = 0.466, P = 0.047; r = 0.405, P = 0.023, respectively), suggesting that their deficiency exerts an influence on the development of ISS. In contrast, family Ruminococcaceae was positively correlated with serum levels of NTproBNP (r = 0.414, P = 0.028). The genera Parabacteroides and Clostridium were negatively correlated with HtSDS (r = -0.468, P = 0.045; r = -0.500, P = 0.047, respectively).

To evaluate whether the altered abundancy of gut Ruminococcaceae and genera Faecalibacterium and Eubacterium and others regulates the SCFAs in the host, Spearman correlation analysis was conducted between different bacterial species and SCFA levels (Figure 5C), and the results showed that the family Ruminococcaceae and genera Faecalibacterium and Eubacterium were significantly positively correlated with the concentrations of total SCFAs and butyrate (0.6> r > 0.4, P < 0.05). These results suggest that the significant decrease in SCFA concentration in the ISS intestine is closely related to the changes in intestinal microflora function.

ROC curve analysis was used to further evaluate the differential bacteria and metabolites to explore their diagnostic value for ISS. The area under the curve (AUC) for genus Faecalibacterium, which showed the most significant reduction in relative abundance in ISS group, was 0.83 (95% CI: 0.68–0.99, P = 0.001), and the cut-off value was 0.60 (sensitivity: 0.94, specificity: 0.78) (Figure 6A). The AUC of genus Parabacteroides, which showed the most significant increment in relative abundance in ISS group, was 0.79 (95% CI: 0.63–0.95, P =0.005). The cut-off value was 0.38 (sensitivity: 0.89, specificity: 0.68) (Figure 6B). We further included both genera Faecalibacterium and Parabacteroides as indices for binomial logistic regression analysis to construct a joint diagnostic model. As shown in Figure 6D, the AUC value, sensitivity, and specificity of the diagnostic model were higher: AUC value was 0.88 (95% CI: 0.76–1, P < 0.001), and cut value was 0.55 (sensitivity: 0.94, specificity: 0.81), indicating an improved diagnostic model.

The diagnostic value of the metabolites in ISS was evaluated, as shown in Figure 6C. The AUC value of acetate was 0.70 (95% CI: 0.52–0.92, P = 0.46) and the cut-off value was 2006.67 (sensitivity: 0.88, specificity: 0.63). The AUC value of butyrate was 0.94 (95%CI: 0.86–0.99, P < 0.001), and the cut-off value was 994.42 (sensitivity: 0.94, specificity: 0.88). The AUC value of valerate was 0.83 (95% CI: 0.69–0.97, P = 0.001) and the cut-off value was 42.51 (sensitivity: 0.56, specificity: 0.82). The AUC value of caproate was 0.76 (95% CI: 0.62–0.95, P = 0.006), and the cut-off value was 9.24 (sensitivity: 0.88, specificity: 0.69). These results indicate that butyrate has a high diagnostic value and could be used as a reference in clinical practice.