Non-significant influence between aerobic and anaerobic sample transport materials on gut (fecal) microbiota in healthy and fat-metabolic disorder Thai adults

Participant’s recruitment, fecal sample collections and metagenomic extraction

Nine healthy and eleven fat-metabolic disorder Thai participants, males and females of the age range 24–43 years, were recruited, and all methods used in this study were in accordance with the guidelines by the ethical approval. The Institutional Review Board, Faculty of Medicine, Chulalongkorn University (no. 735/61) granted the ethical approval for the study. Written informed consent was obtained from all participants in this study. Fecal samples of these twenty total subjects were collected in fecal containers with one aerobic and one anaerobic transport material; therefore, there were 20 aerobic transport samples and 20 anaerobic transport samples (Fig. 1). All forty samples were individually metagenomic extracted, 16S rRNA gene sequenced and qPCR for microbiota and quantitative microbiota analyses. In aspect of sample size (N), the statistically required sample size: N = (p (1-p) z²)/e² was computed, given p at an estimated incidence between aerobic vs. anaerobic microbiota difference of 50%, z score of ±1.44 for 85% confidence interval, and e of 11.5% for margin of error. This yielded an N of 40 (20 aerobic and 20 anaerobic transport samples).

For aerobic transport material, the fecal container was capped, sealed and placed in a plastic bag. For anaerobic transport material, the fecal container was capped, sealed, and placed in a plastic bag with the AnaeroPack-Anaero (Mitsubishi Gas Chemical, Tokyo, Japan) (<0.1% O₂ and >15% CO₂) (van Horn, Warren & Baccaglini, 1997; Wen et al., 2021). The samples were transported on the same day of fecal collection at a cold temperature (≤4 °C) and processed immediately within 24 h for metagenomic extraction using DNeasy PowerSoil Pro Kit (Qiagen, Hilden, Germany) following the manufacturer’s instruction (Wongsaroj et al., 2021; Ondee et al., 2022). The metagenomic DNA was qualified and quantified by agarose gel electrophoresis and nanodrop spectrophotometry (A260 and A260/A280).

16S rRNA gene V3-V5 library preparation and MiSeq sequencing

PCR amplification of the 16S rRNA gene at the V3-V5 region was performed using the universal prokaryotic primers 342F (5′-GGRGGCAGCAGTNGGGAA-3′) and 895R (5′-TGCGDCCGTACTCCCCA-3′) with appended barcode and adaptor sequences (Human Microbiome Project (HMP) Consortium, 2012a; Castelino et al., 2017; Wongsaroj et al., 2021; Dityen et al., 2022). The 342F was used elsewhere and the 895R position was shared with the 909R. The in-silico analysis revealed that the V3-V5 primers could identify bacteria on phylum/class/order/family levels with >77% efficiency, genus 56.6% and species 21.1% (Wang & Qian, 2009; Human Microbiome Project (HMP) Consortium, 2012a; Castelino et al., 2017; Johnson et al., 2019; Darwish et al., 2021; Suwarsa et al., 2021; Wongsaroj et al., 2021; Dityen et al., 2022). Each PCR reaction comprised 1 × EmeraldAmp GT PCR Master Mix (TaKaRa, Shiga, Japan), 0.2 μM of each primer, and 50–100 ng of the genomic DNA in a total volume of 75 µL. The PCR conditions were 94 °C 3 min, and 25 cycles of 94 °C 45 s, 50 °C 1 min and 72 °C 1 min 30 s, followed by 72 °C 10 min. A minimum of two independent PCR reactions were performed and pooled to prevent PCR stochastic bias. Then, the ~640-base pair (bp) amplicon was excised from agarose gel resolution and purified using PureDireX PCR Clean-Up & Gel Extraction Kit (Bio-Helix, Keelung, Taiwan), and quantified using a Qubit 3.0 Fluorometer and Qubit dsDNA HS Assay kit (Invitrogen, Waltham, MA, USA). Finally, 180 ng of each barcoded amplicon product was pooled for sequencing using the Miseq600 platform (Illumina, San Diego, CA, USA), along with the sequencing primers and index sequence (Caporaso et al., 2012; Wongsaroj et al., 2021; Dityen et al., 2022; Ondee et al., 2022), at the Omics Sciences and Bioinformatics Center, Chulalongkorn University (Bangkok, Thailand).

Quantification of total bacteria copy number

The 16S rRNA gene qPCR was performed to quantify total bacteria in copy unit, using universal primers 1392F (5′-CGGTGAATACGTTCYCGG-3) and 1492R (5′-GGTTACCTTGTTAC GACTT-3′), and Quantinova SYBR green PCR Master Mix (Qiagen, Hilden, Germany) in a 20 µL total volume and 1 ng metagenomic DNA (or reference DNA), as previously established (Suzuki, Taylor & DeLong, 2000; Oldham & Duncan, 2012; Wongsaroj et al., 2021). The qPCR thermocycling parameters were 95 °C 5 min, followed by 40 cycles of 95 °C 5 s and 60 °C 10 s. They ended with a 50–99 °C melting curve analysis to validate a single proper amplicon peak (i.e., neither primer-dimer nor non-specific amplification). The reference for copy number computation was Escherichia coli, in which the ~120-bp 1392F-1492R amplicon fragments were cloned into pGEM-T-Easy Vector (Promega, Wisconsin, WI, USA) and the recombinant plasmids were transformed into competent E. coli DH5α for expression (Hanahan, Jessee & Bloom, 1991). The inserted fragments were verified by colony PCR using the primers M13F (on vector) and 1492R (inserted fragment). Ten-fold serial dilutions of the extracted plasmids (10⁵–10¹⁰ copies/µL) were used as the reference standard curves in the bacterial copy number computation as following equation (Smith et al., 2006).

$$\mathrm{C}\mathrm{o}\mathrm{p}\mathrm{y}\mathrm{n}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r}\mathrm{p}\mathrm{e}\mathrm{r}\mu \mathrm{L}={\displaystyle \frac{\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{n}\mathrm{t}\mathrm{r}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}(\mathrm{n}\mathrm{g}/\mu \mathrm{L})\times 6.023\times 10^{23}\left(\mathrm{c}\mathrm{o}\mathrm{p}\mathrm{i}\mathrm{e}\mathrm{s}/\mathrm{m}\mathrm{o}\mathrm{l}\right)}{\mathrm{l}\mathrm{e}\mathrm{n}\mathrm{g}\mathrm{t}\mathrm{h}(\mathrm{b}\mathrm{p})\times 6.6\times 10^{11}(\mathrm{n}\mathrm{g}/\mathrm{m}\mathrm{o}\mathrm{l})}}$$

The qPCR experiments were performed using Rotor-GeneQ (Qiagen, Hilden, Germany). Three replicates were conducted per reaction. The bacteria copy number of each sample was quantified against the reference standard curve by Rotor-Gene Q Series Software (Qiagen, Hilden, Germany).

Bioinformatic and statistical analyses for bacterial microbiota diversity and potential metabolisms

Raw sequences (reads) were processed following Mothur 1.39.5’s standard operation procedures for MiSeq (Schloss et al., 2009) (https://github.com/mothur/mothur/releases/), including removal of (a) reads shorter than 100 nucleotides (nt) excluding primer and barcode sequences, (b) ambiguous bases ≥4, (c) chimera sequences, and (d) homopolymer of >7 homopolymers. The sequences were aligned with the 16S rRNA gene references and taxonomic database SILVA 13.2 (McDonald et al., 2012), and Greengenes 13.8 (Quast et al., 2013) to remove lineages of mitochondria, chloroplasts, eukaryotes, and chimera sequences. Then, the quality sequences were clustered into operational taxonomic units (OTU) with 97% nt similarity (78% for phylum, 88% order, 91% class, 93% family, 95% genus, and 97% species) based on naïve Bayesian taxonomic method with default parameters (Wang et al., 2007; Schloss et al., 2009). Samples were normalized for an equal sequencing depth (7,137 quality sequences per sample). The count of total bacteria copy number from the 16S rRNA gene qPCR data were analyzed along with the percent microbiota composition to yield the quantitative microbiota (the bacterial copy number for each individual OTU) (Vandeputte et al., 2017a, 2017b; Jian et al., 2020; Wongsaroj et al., 2021). Alpha diversity including Good’s coverage index (percent sequence coverage to true estimate), rarefaction curve, Chao1 richness, inverse Simpson and Shannon diversity; and beta diversity including Smith theta (Thetan), Sorenson (Sorabund), Morisita–Horn, Yue and Clayton theta (Thetayc), Bray-Curtis (BC), Jaccard (jclass), and Lennon (Lennon) coefficients, and two-dimension non-metric multidimensional scaling (NMDS), were computed using Mothur 1.39.5 (Schloss et al., 2009; Schloss, 2020). Estimates of the microbial metabolic profiles were determined by PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) based on the reference genome annotations in KEGG (Kyoto Encyclopedia of Genes and Genomes pathways) (Langille et al., 2013), and statistically compared by STAMP (Statistical Analysis of Metagenomic Profiles) (Parks et al., 2014). The differences in microbial metabolic profiles were further analyzed by linear discriminant analysis effect size (LEfSe) method with pairwise Kruskal–Wallis and Wilcoxon tests to identify the microbial metabolic biomarker representing healthy and disease groups. For general statistics, non-parametric multiple t-tests were used and a P-value < 0.05 was considered significant.

Availability of supporting data

The nucleic acid sequences in this study were deposited in the NCBI open access Sequence Read Archive database, accession number PRJNA1020208.

16S rRNA gene sequencing results and percent microbiota compositions

The 16S rRNA gene sequencing yielded 2,365,959 total raw sequences (Table S1: aerobic sample transport 1,517,643 sequences, and anaerobic sample transport 848,316 sequences), and 1,623,517 total quality sequences (aerobic sample transport 1,062,335 sequences, and anaerobic sample transport 561,182 sequences). The average quality sequences per sample were 40,587 ± 24,139 (avg. ± SD), and the numbers of OTUs ranged 5–10 at phylum (Table 1: average 6.80 ± 1.22 OTUs), 55–93 genus, and 77–133 species levels, respectively (Table S1 and 1). The number of OTUs at phylum, genus and species levels were found approximately equal between aerobic and anaerobic sample transports (Table 1: phylum OTUs 6.55 ± 1.19 aerobic, 7.05 ± 1.23 anaerobic; genus OTUs 71.40 ± 10.45 aerobic, 72.70 ± 11.29 anaerobic; and species OTUs 101.15 ± 16.83 aerobic, 101.60 ± 15.67 anaerobic). Following the successfully high number of quality sequences, the Good’s coverage (estimated percent sequence coverage to true diversity) of all samples were above 99.5% at phylum, genus and species level OTUs: avg. 100% phylum, 99.82% genus and 99.72% species (Tables 1 and S1). Once data normalization was performed of all samples, each to the same sequencing depth, the Good’s coverages remained average >99% and the rarefaction curves were relative plateau (Fig. S1). The data disclosed that the further microbiota bioinformatic analyses had no bias from various quality sequencing numbers per sample.

The percent bacterial compositions at phylum, genus, and species levels across all participants were compared between aerobic vs. anaerobic sample transport materials, and no statistical difference in the phylum/genus/species was found (AMOVA, P > 0.05) (Fig. 2). Five major phyla, ranging from Firmicutes as the top abundant (averagely, 52.03 ± 17.30%), Bacteroidetes (24.32 ± 14.11%), Proteobacteria, Actinobacteria, to Fusobacteria, were presented. The latter three phyla accounted for an average <24%. Twenty-two bacterial genera (equating 24 bacterial species OTUs), excluded <1% genus, or species, were revealed and the individual percent genus (or species) was compared between aerobic vs. anaerobic sample transport materials: no statistical difference were found (t-test, P > 0.05) (Table S2). The OTU compositions indicated no statistical difference in microbiota percents and compositions at phylum, genus and species levels, between aerobic and anaerobic sample transport groups.

Quantitative microbiota composition analyses between aerobic and anaerobic sample transport groups

Following the quantification of bacteria by the universal 16S rRNA gene qPCR, the number of bacterial counts and the quantitative microbiota compositions could be analyzed. The quantity of bacterial counts was not significantly different between aerobic and anaerobic sample transport groups, although slightly lower for the aerobic sample transport group (Fig. 3A: P = 0.057). Noted that the relatively low in the aerobic sample transport group was due to ID3a and the relatively high in the anaerobic sample transport group was due to ID1an; if except these two, the average bacterial counts of both groups will even become closer to each other and P value increases (Fig. S2).

Next, individual bacterial species corresponding to obligate (or strictly) anaerobes that consisted of five bacterial species and facultative anaerobes that consisted of three species were quantitatively compared. No statistically significant difference in quantity was pointed in these bacterial species between aerobic and anaerobic sample transport groups (Fig. 3B). In detail, the obligate anaerobic Bacteroides spp. were found most dominated than other obligate anaerobic bacterial genera in both groups and presented in approximately comparable counts, followed by Prevotella, Faecalibacterium, Oscillospira, Bifidobacterium, and the facultative anaerobic Haemophilus, Streptococcus and Enterococcus, respectively. Nonetheless, the slight but non-statistically significant higher counts of obligate anaerobic bacteria were shown. Still, this trend was minute and found inconsistent for facultative anaerobic bacteria genera (Fig. 3B), highlighting the differences in obligate vs. facultative oxygen requirement effect yet at the non-significant statistic. Overall, the percent microbiota composition and the quantitative microbiota did not demonstrate significant difference between aerobic and anaerobic sample transport materials. Subsequently, the alpha diversity by OTU species richness (OTUs and Chao1) and OTU species diversity (inverse Simpson and Shannon) showed very high P values between 0.3827 and 0.9497 (Fig. 4), and the beta diversity among individual samples belonging to aerobic and anaerobic sample transport groups showed no separate clustering pattern (Fig. 5A). Noted that the detail analyses of alpha diversity at OTU phylum and genus levels were also analyzed. No statistic differences were found (P > 0.05) (Fig. S3). Additionally, other beta diversity coefficients, such as Sorabund, Morisita-Horn, Thetayc and Bray-Curtis, were computed and all dissimilarity coefficient indices did not separate the microbiota community differences between aerobic and anaerobic sample transport groups (Table S3: P > 0.05). Meanwhile, we further classified the samples into healthy and unhealthy categories, and the alpha diversities showed relatively no difference between aerobic and anaerobic sample transports (Fig. S3E).

Quantitative microbiota analyses between healthy and fat-metabolic disorder groups

When we analyzed the quantitative microbiota structure differences by different beta diversity coefficients, we found the statistical difference between healthy vs. fat-metabolic disorder (from now on referred as “unhealthy”) groups (Fig. 5B: P = 0.02). The differences were found when considering only aerobic healthy vs. unhealthy, anaerobic healthy vs. unhealthy, and combined aerobic+anaerobic healthy vs. unhealthy. Supportively, the clinical parameters corresponding to fat-metabolic disorders demonstrated statistically (P < 0.05: age, liver stiffness, GGT, BMI, TC, AST, ALT, TG, LDL, and CAP) and non-statistically (P > 0.05, HDL) associated the same direction with the unhealthy microbiota community structure (Fig. 5C). Figure 5D exhibited bacterial species that significantly associated with unhealthy community structure patterns such as Prevotella, Haemophilus and Bacteroides plebeius; and healthy community structure such as Bifidobacterium, Ruminococcus and Clostridium.

Furthermore, the low-abundance OTUs of <1% and non-shared inter-individual microbiota were tested and filtered out (remaining as “core microbiota”) for the NMDS analysis. The result remained consistent, demonstrating no statistical difference in quantitative core microbiota between aerobic and anaerobic sample transport groups (Fig. 5E: P = 0.87), yet the statistical difference between healthy and unhealthy groups (Fig. 5F: P = 0.019). This finding might infer the importance in the core microbiota pattern that aligned the unhealthy microbiota association with the fat-metabolic disorder.

Metabolic function prediction levels via quantitative profiles of prevalent health-associated bacteria, and microbial metabolic function species biomarkers for healthy and fat-metabolic disorder groups

The metabolic potentials of the potentially important bacteria were analyzed. These included Bacteroides, Prevotella, Megamonas, Bifidobacterium, Hemophilus, Clostridium, Ruminococcus and Pasteurellaceae (Wu, Bushmanc & Lewis, 2013; Schirmer et al., 2019; Sun et al., 2020; Sabo & Dumitrascu, 2021). The generally most active microbial-related functions were metabolism pathway (49.92%: primarily amino acid and carbohydrate metabolisms followed by energy, cofactors and vitamins, lipid and xenobiotics biodegradation metabolisms), 19.94% in genetic information processing, 16.22% in environmental information processing, 3.11% cellular process, 0.91% human diseases, 0.65% organismal systems, and 5.09% poorly characterized. The OTUs of Bacteroides and Prevotella copri represented the topmost varying functional metabolisms (Fig. 6A). Meanwhile the functional redundancy among bacterial OTUs, the relative abundances of these health-associated bacteria showed the dynamic functions with some distinguished categories of metabolisms, cellular process, and genetic information processing between healthy and fat-metabolic disorder groups. For instance, the relatively more abundance of amino acids, carbohydrate and energy metabolism functions, cellular processes, genetic information processing, and human diseases were reported in the fat-metabolic disorder group. Prevotella copri, Prevotella stercorea, and Bacteroides plebeius, were estimated to have more diverse and abundant functions in the fat-metabolic disorder group, while Bacteroides and Bifidobacterium longum, were estimated to be more diverse and abundant in the healthy group (Fig. S4). These microbial metabolism differences between groups allowed LEfSe to identify the specific microbial metabolic functions along the bacterial species as the biomarkers to differentiate between healthy vs. fat-metabolic disorder groups, with statistical P values. Prevotella copri and Bacteroides plebeius were the biomarkers for the fat-metabolic disorder. Their microbial metabolic functions included many functions involved in diseases (immune system diseases, metabolic diseases, and neurodegenerative diseases). In contrast, the healthy group showed a greater variety of bacterial species and their associated metabolic functions when compared to the unhealthy group. This supports the existence of diverse microbial-related metabolic functions in healthy human guts. It was noted that the commonly reported functions were related to metabolism and organismal systems pathways, meanwhile the human disease pathway was rare in the healthy group (Fig. 6B).