Alteration of the gut microbiota associated with childhood obesity by 16S rRNA gene sequencing

Ethics statement

The study was conducted in accordance to the Declaration of Helsinki revised in 2013 and approved by the ethics committee of Zhongda Hospital, Southeast University (approval number: 2017ZDSYLL109-P01). Participation in this study was voluntary, and all parents gave written informed consent.

Research object

A total of 51 volunteers of both sexes (27 males and 24 females) were recruited from Nanjing, China. Subjects were strictly categorized into the normal weight group (n = 23) or the obese group (n = 28) according to the inclusion and exclusion criteria.

Inclusion criteria: (1) The age range for participation was 6-11 years old. (2) In accord with the determination standard of normal weight and obesity of children, the standard was “Body mass index cut-offs for overweight and obesity in Chinese children and adolescents aged 2-18 years” (Li et al., 2010) formulated by the Department of Growth and Development, Capital Institution of Pediatrics (Table 1). (3) Willing to participate in the study and obtain the consent of the guardian, voluntarily taking the subject and signing the informed consent form.

Exclusion criteria: (1) Antibiotics have been used in the past 4 weeks. (2) Gastrointestinal dysfunction, previous gastrointestinal disease history or diarrhea, abdominal distension, abdominal pain or constipation within the past 4 weeks. (3) Trauma, serious infection, and infectious diseases. (4) Hereditary obesity. (5) Drug-induced obesity. (6) Endocrine disorders and metabolic diseases.

Sample collection

Fecal samples were collected and well-sealed in sterile boxes, immediately frozen in a refrigerator and transported to school the next morning and stored at −20 °C before being transferred in insulating polystyrene foam containers to the Key Laboratory of Environmental Medical Engineering and Education Ministry, where samples were stored at −80 °C until further analysis.

DNA extraction and PCR amplification

Genomic DNA was extracted according to the specifications of TIANamp Stool DNA kit (TIANGEN, China, Cat. DP328), applied to all feces samples. The integrity and purity of DNA were detected by 1% agarose gel electrophoresis, while the concentration and purity of DNA were detected by NanoDrop One. The V4 region of the bacterial 16S rDNA was amplified by PCR with the specific primers 515F (5′-GTGCCAGCMGCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) labeled in a 12 bp barcode. PCR amplification was conducted using primers with a barcode and Premix Taq under the following conditions: 5 min at 94 °C for initialization, 30 s denaturation at 94 °C for 30 cycles, 30 s annealing at 52 °C, and 30 s extension at 72 °C, and a final 10 min elongation at 72 °C. The fragment length and concentration of PCR products were detected by 1% agarose gel electrophoresis, and samples with a main band length in the range of 290-310 bp were selected. After comparing the concentrations of PCR products using GeneTools Analysis Software (Version 4.03.05.0, SynGene), the required volume of each sample was calculated according to the principle of equal quality, and each PCR product was mixed and recovered by the An E.Z.N.A.^® Gel Extraction Kit (OMEGA, USA, Cat. D2500). TE buffer was used to elute and recover the target DNA fragments.

Sequencing and data processing

Libraries were built according to the NEBNext^® Ultra™ DNA Library Prep Kit for Illumina^® and sequenced on an Illumina HiSeq2500 platform (Caporaso et al., 2011), and then 250 bp paired-end reads were generated. Trimmomatic software (V0.33, http://www.usadellab.org/cms/?page=trimmomatic) (Bolger, Lohse & Usadel, 2014) was used to filter the quality of the raw reads data at both ends. At the same time, with reference to the barcode and primer information at both ends of the sequence, Mothur software (V1.35.1, http://www.mothur.org) (Schloss et al., 2009) was used to distribute the sequence to corresponding samples; the allowed mismatch number of barcodes was 2, and the maximum mismatch number of primers was 3. Then, after quality control, barcode and primers were removed to obtain paired-end clean reads. For double-ended sequencing data, it was necessary to splice each pair of paired-end reads using FLASH software (V1.2.11, https://ccb.jhu.edu/software/FLASH/) (Magoc & Salzberg, 2011) in terms of the overlap between paired-end reads to splice the paired reads into a sequence and to filter out nonconforming tags in order to collect the original spliced sequence (raw tags). The minimum overlap length was set to 10 bp and the maximum allowable mismatch ratio in the overlap region of the spliced sequence was set to 0.1. Mothur software was used to carry out quality control and filter the spliced sequences to obtain effective spliced fragments.

OTU and species community analysis

The USEARCH software (V8.0.1517, http://www.drive5.com/usearch/) (Edgar, 2010) was used to cluster all clean tags for all samples. By default, the sequence was clustered into an operational taxonomic unit (OTU) with 97% identity. Singleton OTUs were removed with usearch (http://www.drive5.com/usearch/manual/chimera_formation.html), and chimeric sequences were removed with UCHIME (http://www.drive5.com/usearch/manual/uchime_algo.html) (Edgar et al., 2011). The assign_taxonomy.py script and Ribosomal Database Project (RDP) Classifier method (Cole et al., 2007) in Quantitative Insights Into Microbial Ecology (QIIME, version 1.9.1) were used to obtain species annotation information. The number of valid tag sequences (No. of seqs) and the OTU taxonomic comprehensive information table (otu_table) were obtained by removing chloroplast and mitochondrial sequences as well as OTUs and the tags that could not be annotated at the set limit. Based on otu_table, R software (V2.15.3) (R Core Team, 2013) was used to calculate the annotation proportion of OTUs at each classification level, and the sequence of each sample at each classification level was obtained to form a column chart. All values greater than 0 in each column of values were recorded as 1 and summed, which was the total OTUs of each sample. On the basis of the normalized otu_table, the common or unique OTU analysis was conducted with the ggplot2 package in R software, and meanwhile the OTU triangulation was drawn using the ggtern package to show the common or unique OTUs and their abundance between the two groups. The pheatmap package in R software was used to carry out the cluster analysis between samples and species. With the phylogenetic relationship and relative abundance information of each OTU in the samples, the species annotation results of a single sample was visualized by KRONA software (http://sourceforge.net/projects/krona/) (Ondov, Bergman & Phillippy, 2011). Using GraPhlAn software (http://huttenhower.sph.harvard.edu/graphlan) (Asnicar et al., 2015), a single sample OTU annotation circle graph based on GraPhlAn was obtained. Representative OTU sequences with the top 50 overall relative abundance and genus classification information were selected, and Mafft software (Katoh et al., 2002) was used to carry out a multisequence comparison. FastTree software (http://microbesonline.org/fasttree/) (Price, Dehal & Arkin, 2009) was used to build trees simultaneously combining the relative abundance of each OTU and species annotation confidence information of representative sequences, and the ggtree software package (Yu et al., 2017) was used to carry out the visual display. OTU abundance information were normalized using a standard of sequence number corresponding to the sample with the least sequences. Based on the normalized data, alpha diversity and beta diversity were all performed in the following paragraphs.

Alpha diversity analysis

According to the normalized OTU abundance table, the alpha_diversity.py script (http://huttenhower.sph.harvard.edu/graphlan) in the QIIME software package (version 1.9.1) (Caporaso et al., 2010; Kuczynski et al., 2011) was used to calculate four diversity indices (Chao et al., 2010). According to the OTU abundance table, the alpha_rarefaction.py script (http://qiime.org/scripts/alpha_rarefaction.html) in the QIIME software package was used to calculate dilution curve data of four diversity indices, and the vegan package (Dixon, 2003) was used to draw the dilution curve. Based on the OTU abundance table, the plot_rank_abundance_graph.py script (http://qiime.org/scripts/plot_rank_abundance_graph.html) in the QIIME software package was selected to calculate the rank debt curve index. The data from the specaccum species accumulation curve analysis was calculated and mapped in line with the normalized OTU abundance table.

Beta diversity analysis

According to the normalized OTU abundance table, the beta diversity distance was calculated using jackknifed_beta_diversity.py script (http://qiime.org/scripts/jackknifed_beta_diversity.html) in QIIME software. Principal Coordinate Analysis (PCoA) was performed to get principal coordinates and visualize from complex data (Lozupone et al., 2011). A distance matrix of unweighted Unifrac among samples obtained before was transformed to a new set of orthogonal axes, the maximum variation factor is demonstrated by first principal coordinate, the second maximum one by the second principal coordinate, and so on. PCoA analysis was displayed by QIIME2 and ggplot2 package (Ginestet, 2011). Upgma_cluster.py script (http://qiime.org/scripts/upgma_cluster.html) in QIIME software was applied to build the cluster tree of samples using the UPGMA cluster analysis method. Based on the normalized OUT table (otu_table_subsampled), Nonmetric Multidimensional Scaling (NMDS) analysis was performed with the vegan package and displayed with the ggplot2 package in R software.

LEfSe and PICRUSt analysis

The linear discriminant analysis (LDA) effect size (LEfSe) analyses were performed on the website http://huttenhower.sph.harvard.edu/galaxy (Segata et al., 2011). For OTUs with an average abundance in all samples that was greater than 0.1%, abundances were normalized to the sum of the values per sample in 1 million and then subjected to LDA. The cut off value was the absolute LDA score (log10) >2.0. Functional capacity of gut microbiota was predicted using the Phylogenetic Investigation of Communities by Reconstruction of Unobserved States (PICRUSt) (Galaxy Version 1.0.0) (Langille et al., 2013). The closed reference OTU table was generated from quality control reads in QIIME against the Greengenes reference sequence database. Closed OTU-table drawn by QIIME was compared with KEGG (Kyoto Encyclopedia of Genes and Genomes) database to obtain functional predictions. PICRUSt predictions were categorized as levels 1-3 into KEGG pathways.

Statistical analysis

IBM SPSS 23.0 software were applied in all statistical analyses. Depending on the normality of the underlying data drawn from the Shapiro-Wilk test, differences between groups were examined by two-tailed t-test or Wilcoxon rank-sum test. Test results at an alpha of P < 0.05 were considered statistically significant.

Study population

In this study, we categorized 28 obese children into the obese group and 23 healthy children into the normal weight group to analyze the two group’s gut microbial composition. The two groups revealed no statistical difference regarding gender ratio and age (Table 2).

Sequencing data

A total of 2,349,074 raw sequence reads were obtained from the 51 subjects. After a series of quality filtering, 2,004,646 classifiable 16S rRNA gene sequences were obtained, and the average number of sequences for each individual was 39,307 (ranging from 25,813 to 55,847). All sequences were clustered with representative sequences, and a 97% sequence identity cut-off was used.

The Venn diagram (Fig. 1A) demonstrated the shared and exclusive communities between the groups. The total number of OTUs obtained was 1,213, among which 449 OTUs were shared by both groups, 211 genera were specific to the normal weight group and 104 genera were specific to the obese group, showing a decrease in the gut microbial diversity as the body weight increased.

Estimation of the alpha diversity and beta diversity

Rarefaction curves drawn on the observed species and Shannon indices (Figs. 1B and 1C) indicated that although deeper sequencing may reveal rare OTUs, the majority of microbial diversity had been captured. These curves also revealed that the alpha diversity in samples from the normal-weight group was the highest and that in obesity group was the lowest. To confirm its validity, we calculated Chao1 indices by Wilcoxon rank-sum test (Fig. 2), finding that the mean microbial abundance between two groups decreased significantly (P = 0.006, Wilcoxon rank-sum test). Shown in the rank-abundance curve, the species abundance and uniformity of gut microbiota in the obese group were lower than those in the normal weight group (Fig. 1D).

The alpha diversity indices (the Chao1, observed species, PD whole tree and Shannon index) were used to describe alpha diversity (Fig. 2). Significant differences were found between the normal weight group and the obese group (P < 0.001; P <  0.001; P < 0.001; and P = 0.008, respectively, Wilcoxon rank-sum test), showing a higher abundance and diversity in the normal weight group than in the obese group. In addition to the alpha diversity evaluation, the unweighted UniFrac analysis was applied to compare similarities among gut microbial communities (beta diversity). Looking into the bacterial composition profiles, we observed significant differences between the normal weight group and the obese group (R ² = 0.054, P = 0.001, adonis analysis). NMDS and PCoA based on the abundance of OTUs revealed differences in the microbial composition (Fig. 3). Specifically, an evident clustering was identified for subjects of normal weight and obesity. The observation revealed significant differences in gut microbial community structure between the normal weight group and the obese group, with two principal component scores accounted respectively for 11.59% and 5.57% of the total variations. Moreover, separation between the two groups was particularly obvious. Representing the intestinal microbial composition, data points for subjects of normal weight clustered at the right and obese ones at the left.

The relative abundance of fecal bacterial communities

Statistics of the OTUs suggested the relative abundance of the bacteria at the categorization of phylum, class, order, family and genus. The results showed that fecal bacterial composition differed between the two groups.

Bacteroidetes was the most predominant phylum, contributing 55.48% and 51.35% of the gut microbiota in the normal weight group and the obese group, respectively, followed by Firmicutes, contributing 37.93% and 36.18%, respectively (Figs. 4A and 4B). Proteobacteria, Fusobacteria, Verrucomicrobia and Actinobacteria constituted the next most dominant phyla. Microbial compositions showed high interindividual variability, among which Bacteroidetes accounted for 21.01-73.78%, and Firmicutes accounted for 18.72-59.49% among all individuals. In regard to phyla, proportion of Bacteroidetes (51.35%) in the obese group was lower than that in the normal weight group (55.48%) (P = 0.030, Wilcoxon rank-sum test). No statistical differences were revealed in Firmicutes (P = 0.436, Wilcoxon rank-sum test) and the Firmicutes/Bacteroidetes (P = 0.983, Wilcoxon rank-sum test) between the two groups.

Analysis of the relative abundance of bacterial taxonomic groups showed that the fecal microbiota was dominated by fourteen genera: Bacteroides (mean relative abundance, N, 41.58%; F, 40.06%) (In the following analysis, ”N” is short for “Normal weight” and “F” for “Fat”, namely “Obesity”), Faecalibacterium (N, 7.03%; F, 10.76%), Prevotella (N, 4.72%; F, 6.19%), Oscillospira (N, 4.50%; F, 1.61%), Dialister (N, 4.03%; F, 0.91%), Parabacteroides (N, 4.01%; F, 2.39%), Sutterella (N, 2.68%; F, 5.60%), Roseburia (N, 1.94%; F, 2.32%), Ruminococcus (N, 1.91%; F, 1.13%), Phascolarctobacterium (N, 1.64%; F, 3.43%), Lachnospira (N, 1.55%; F, 3.23%), Escherichia (N, 1.39%; F, 1.32%), Megamonas (N, 0.54%; F, 1.79%), Haemophilus (N, 0.32%; F, 1.96%). Faecalibacterium, Phascolarctobacterium, Lachnospira, Megamonas, and Haemophilus in the obese group were significantly more abundant than those in the normal weight group (P = 0.048, P = 0.018, P <0.001, P = 0.040 and P = 0.003, respectively, Wilcoxon rank-sum test). As for Oscillospira and Dialister, they took lower proportions in the obese group when compared to the normal weight group (P = 0.002, two-tailed t-test; and P = 0.002, Wilcoxon rank-sum test). For the remaining genera, the gut microbiome sequencing data did not indicate significant differences in the abundance for the obese group compared to the normal weight group.

Phylogenetic and taxonomic profiles of gut microbiota and PICRUSt analysis

To analyze the statistical differences in microbial communities between the normal weight group and the obese group, we compared OTUs with the LEfSe analysis. The histogram reflected the LDA scores that were computed for the features at the OTU level (Fig. 5A). Cladograms of the taxa with LDA values >2.0 were depicted in Fig. 5B.

Upon analyzing the fecal samples, a total of 55 species of bacteria were abundant in the normal weight group and 45 species were abundant in the obese group. As shown in Fig. 5A, the relative abundance of taxonomic groups showing an LDA score greater than 10⁴ was summed for the normal weight group (Oscillospira and Bacteroides uniformis) and the obese group (Proteobacteria, Prevotella copri, Alcaligenaceae, Sutterella, Betaproteobacteria, and Burkholderiales). Multiple genera were found in significantly high abundances in the normal weight group. These included Oscillospira, Ruminococcus, Prevotella, Adlercreutzia, Sporobacter, Bifidobacterium, Clostridium, Desulfovibrio, Bilophila, Christensenella, Alistipes, Anaerotruncus, Eubacterium, Holdemania, Oxalobacter, Defluviitalea, and Collinsella. The genera that were enriched in the obese group were Turicibacter, Campylobacter, Actinobacillus, Aggregatibacter, SMB53, Rothia, Granulicatella, Streptococcus, Veillonella, Megamonas, Fusobacterium, Phascolarctobacterium, Haemophilus, Lachnospira, and Sutterella.

Cladograms (Fig. 5B) were generated from the LEfSe analysis, which showed the most differentially abundant taxa enriched in microbiota with green for the normal weight group and red for the obese group. The diameter of each circle is proportional to its abundance. The obese group showed a significant decrease in the phylum Actinobacteria, phylum Tenericutes, class Deltaproteobacteria, class Erysipelotrichia and some taxonomic groups belonging to the order Bacteroidales, such as family Paraprevotellaceae, Barnesiellaceae, S24-7, and Rikenellaceae and a greater abundance of the phylum Proteobacteria, phylum Fusobacteria and class Bacilli when compared with the normal weight group (Fig. 5B).

To study the changes of gut microbial function in obese children, we adopted PICRUSt. Based on KEGG database, PICRUSt revealed a total of six biological metabolism pathways at Level 1 pathways: metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, human diseases. Among them, metabolism, genetic information processing and environmental information processing dominated, accounting for 46.89%-51.41%, 18.33%-22.12% and 10.22%-14.87%, respectively. These six pathways in the obese group were lower than those in the normal weight group. Meanwhile, the secondary function of the predicted gene was analyzed, finding it consisted of 39 sub-functions including membrane transport, carbohydrate metabolism, amino acid metabolism, replication and repair, energy metabolism, translation, cellular processes and signaling. Within the 39 predicted functional categories at KEGG pathway hierarchy level 2, except immune system diseases, neurodegenerative diseases and signal transduction, the remaining 36 sub-functions within the predicted gene all decreased in obese children.