Measuring the absolute abundance of the microbiome by adding yeast containing 16S rRNA gene from a hyperthermophile.

High-throughput sequencing (HTS) of 16S rRNA gene amplicons provides compositional information regarding the microbial community, but not the absolute abundance of the bacteria. We aimed to develop a standardized method for quantifying the absolute abundance of bacteria in microbiome studies. To demonstrate the utility of our approach, we quantified the number of bacteria from the compositional data of the fecal and cecal microbiomes. The 16S rRNA gene of a hyperthermophile, Thermus aquaticus, was cloned into Pichia pastoris (yeast) genome, and an equivalent amount of the yeast was added to the stool and cecal samples of mice before DNA extraction. 16S rRNA gene library construction and HTS were performed after DNA extraction. The absolute abundances of bacteria were calculated using T. aquaticus reads. The average relative abundances of T. aquaticus in the five stool and five cecal samples were 0.95% and 0.33%, respectively, indicating that the number of bacteria in a cecum sample is 2.9 times higher than that in a stool sample. The method proposed for quantifying the absolute abundance of the bacterial population in this study is expected to overcome the limitation of showing only compositional data in most microbiome studies.

INTRODUCTION

High‐throughput sequencing (HTS) of 16S rRNA gene amplicon is a reliable and widely used method used in microbiome studies (Poretsky et al., 2014). However, this technique only describes the bacterial composition (relative abundance) and does not quantify the absolute abundances in the samples because the DNA concentration of the 16S rRNA gene is compensated during the construction of the 16S rRNA gene library (Gloor et al., 2017).

Researchers investigating the microbiome have the desire to quantify not only the percentage of each bacterium in the sample but also the absolute quantity of each bacterial species. Bacterial quantification is essential in numerous studies, that is, comparison of the stool microbiome between patients with diarrhea and healthy controls, and comparison of the microbiome in different environments, such as soil and water. In a study of patients with Crohn's disease, absolute abundance of the fecal microbiome was associated with disease phenotype (Vandeputte et al., 2017).

Recently, there have been studies to estimate the absolute abundance of the microbiome by adding internal standards (Harrison et al., 2021). Cellular internal standards use a specific bacterial species mixed into the target sample during the HTS procedure, whereas DNA internal standards use genomic DNA or synthetic DNA. Cellular internal standards may have similar traits to those of microbial organisms in the sample, but they might be more difficult for quantitation and less reproducible than DNA internal standards (Venkataraman et al., 2018). An ideal internal standard should have similar traits to focal organisms, be easily cultured, not occur in the biological samples, not have copy number variation, be easily quantified, and be accurate and reproducible (Harrison et al., 2021).

Thermus aquaticus is a popular hyperthermophile bacterium that was first discovered in the hot spring of Yellowstone National Park (Brock & Freeze, 1969). This bacterium survives at a temperature of 50–80℃ and therefore is not found in the normal environment or animal gut. T. aquaticus is the source of Taq DNA polymerase, an important component of the polymerase chain reaction (PCR) (Chien et al., 1976; Saiki et al., 1988).

This study is aimed at devising a new standardized method for quantifying bacteria in samples when analyzing 16S rRNA gene amplicons using the Illumina MiSeq platform. To fulfill the study objectives, the T. aquaticus 16S rRNA gene was cloned into the Pichia pastoris (yeast) genome for stable amplification. An equivalent amount of this yeast was added to the initial stool and cecum contents of the mice. DNA extraction, 16S rRNA gene library construction, and HTS were performed thereafter for the microbiome analysis. The number of bacteria present in the initial stool and cecum samples was quantified based on the relative abundance of T. aquaticus in the samples.

MATERIALS AND METHODS

Insertion of the T. aquaticus YT 1 16S rRNA gene into the P. pastoris GS115 genome

The complete 16S rRNA gene of T. aquaticus YT 1 (1470 bp, NCBI Reference Sequence: NR_025900.1) was synthesized and cloned into pUC57 (Bio Basic Inc.). To clone this gene into pPIC9, primers were designed to make the PCR amplicon have Xho I and Not I restriction sites on its 5′ and 3′ end, respectively. The primers are as follows: forward, 5′‐AAAAGTCGAGGAATTCGAGCTCGGTACC‐3′3 and reverse, 5′‐AAAAGCGGCCGCATGATTACGCCAAGCTTGC‐3′. Then, the 16S rRNA gene of T. aquaticus was cloned into pPIC9 using Xho I and Not I restriction enzymes. The plasmid containing the 16S rRNA gene of T. aquaticus was linearized using the Sac I restriction enzyme and integrated into the P. pastoris GS115 genome using the Pichia Expression Kit (Invitrogen).

Sanger sequencing was used for sequencing the insert in P. pastoris using the 5′hAOX1 primer 5′‐GACTGGTTCCAATTGACAAGC‐3′ and 3′aAOX1 primer 5′‐GCAAATGGCATTCTGACATCC‐3′. To detect the integration location of the 16S rRNA gene of T. aquaticus in the yeast chromosome, Sanger sequencing was performed using the forward primer 5′‐CTGCTGCCTCCCGTAGG‐3′ and reverse primer 5′‐CTACTGGCTTGGCCATAATT‐ 3′. The forward primer is the reverse complement sequence of the forward primer that we used in the HTS of the 16S rRNA gene V3‐V4 region. The reverse primer is a part of the polyamine oxidase gene in the yeast chromosome 4 that is located close to the alcohol oxidase 1 gene, the integration target region.

P. pastoris culture and aliquot

Pichia pastoris was cultured in yeast extract–peptone–dextrose (YPD) medium. One unit of P. pastoris harboring T. aquaticus 16S rRNA gene was defined as 2 million cells of P. pastoris. One unit of P. pastoris was added to fecal and cecal samples of the mice just before DNA extraction. One stool was used as the fecal sample, whereas the entire content of one cecum was used as the cecal sample.

Stool and cecum contents of mice

Five female C57BL/6 mice (7 weeks of age) were obtained from Orient Bio Inc. Mice were euthanized, and the cecum and stool were sampled. All animal studies were approved by the Department of Laboratory Animal Resources Committee of Yonsei University College of Medicine (No. 2020–0056). Animal experiments were performed in pathogen‐free conditions and following standard management practices.

DNA extraction

DNA was extracted from stool and cecum samples of the mice using FastDNA™ SPIN Kit for Soil (MP Biomedicals), as per the manufacturer's instructions.

Serial dilution of stool suspension

Twenty fecal pellets from 10 mice were placed in a 15 ml conical tube and suspended using 5 ml phosphate‐buffered saline, and twofold serial dilutions were performed. DNA extraction was performed using 500 μL of the diluted samples (x1, x1/2, x1/4, x1/8) containing one unit each of P. pastoris.

HTS of 16S rRNA gene

The V3–V4 region of the 16S rRNA gene was amplified by PCR using the bacterial universal primer pair (forward primer: 5'‐ TCGTCGGCAGCGTCAGATGTGTATAAGAGACAGCCTACGGGNGGCWGCAG‐3’; reverse primer: 5'‐ GTCTCGTGGGCTCGGAGATGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC‐3’). A limited cycle amplification step was included to add multiplexing indices and Illumina sequencing adapters. The libraries were normalized, pooled and sequenced using the Illumina MiSeq platform (Illumina MiSeq V3 cartridge [600 cycles]; Illumina), according to the manufacturers’ instructions.

Bioinformatics and statistics

All analyses were performed with EzBioCloud, a commercially available ChunLab bioinformatic cloud platform for microbiome research (https://www.ezbiocloud.net/) (Kim et al., 2019). Raw reads were processed through a quality check, and low quality (

Whole‐genome sequencing (WGS) of the yeast

WGS was performed to confirm the copy number of the integrated 16S rRNA gene in the yeast. WGS and bioinformatics analysis were performed at Macrogen Inc. The sequencing libraries were prepared according to the manufacturer's instructions of the TruSeq DNA Nano Sample Preparation Kit (Illumina). Briefly, fragmentation of 100 ng of the yeast genomic DNA was performed using adaptive focused acoustics (Covaris Inc.), and the fragmented DNA was end‐repaired to create 5′‐phosphorylated, blunt‐ended dsDNA molecules. Following end‐repair, DNA was size‐selected using the bead‐based method. These DNA fragments undergo the addition of a single “A” base and ligation of TruSeq indexing adapters. Next, we performed sequencing using the HiSeq™ platform (Illumina).

After sequencing, Trimmomatic v0.36 was used to remove adapter sequences and low‐quality reads for reducing biases in the analysis (Bolger et al., 2014). In the sample, filtered data were mapped using BWA v0.7.17 with the mem algorithm to the reference genome (Li & Durbin, 2010). The reference genome was chromosome 4 of P. pastoris GS115 (Komagataella phaffii) sequence (NCBI GenBank: CP014718.1) modified to insert the pPIC9 plasmid containing the 16S rRNA gene into the AOX1 promoter region in the yeast genome (CP014718.1:1592211–1593034). After read mapping, duplicated reads were removed with Sambamba v0.6.7 (Tarasov et al., 2015). The genome coverage and mapping ratio of mapped reads on the reference genome were calculated. The average read depth of a certain region was defined as the sum of the depths of all positions in that region divided by the number of bases in that region.

Quantitative PCR

Quantitative PCR was performed to check the change in the amount of bacterial 16S rRNA gene in the yeast over time. The yeast was fully cultivated in 10 ml of YPD medium for 24 h, and 5 μl was inoculated into 10 ml of fresh YPD medium (n = 5). We collected the yeasts, calculated colony‐forming units (CFU) and extracted DNA for three consecutive days. Then, quantitative PCR was performed using primers for the arginosuccinase gene of the yeast (forward primer 5′‐ACCCGTGAACATGCTTTGCT‐3′ and reverse primer 5′‐CACCATTCTCTCAAGCTCGT‐3′) (Krainer et al., 2012) and primers for the bacterial 16S rRNA gene (BACT1369 5′‐CGGTGAATACGTTCYCGG‐3′ and PROK1492R 5′‐GGWTACCTTGTTACGACTT‐3′) (Buchan et al., 2009) with AMPIGENE® qPCR Green Mixes (ENZO Life Sciences).

RESULTS

The 16S rRNA gene of T. aquaticus inserted into the P. pastoris genome was 100% identical to the NCBI reference sequence, NR_025900.1 (Figure 1a). Then, PCR was used to confirm where the integrated 16S rRNA gene was inserted into the yeast chromosome. PCR was performed using the primer set capable of amplifying the inserted 16S rRNA gene and the polyamine oxidase gene of yeast chromosome 4 because this gene is located next to the AOX1 promoter, which was the suspected integration site (Figure 1b). The PCR product showed the expected band size (1887 bps) in gel electrophoresis, and DNA sequencing (Sanger method) confirmed that it contained the 16S rRNA gene fragment of T. aquaticus, an AOX1 promoter, and the polyamine oxidase gene fragment from the yeast chromosome (Figure 1b). WGS was performed to confirm the copy number of the integrated 16S rRNA gene in the yeast. The average read depth of the total genome of the yeast was 237.88 and that of the integrated 16S rRNA gene of T. aquaticus was 256.21, which means that a single copy of the 16S rRNA gene was inserted into the yeast genome because P. pastoris GS 115 is in a haploid state. This result is in line with the explanation of the kit's protocol that more than 90% of integration cases have one copy of an insert.

FIGURE 1

DNA sequence and the location of the Thermus aquaticus 16S rRNA gene in Pichia pastoris. (a) The inserted T. aquaticus 16S rRNA gene was 100% identical to the NCBI reference sequence NR_025900.1. (b) Integration of the 16S rRNA gene into the AOX1 promoter in the yeast genome was confirmed by PCR and Sanger sequencing using a primer set encompassing the 16S rRNA gene of the pPIC9 and the polyamine oxidase gene from yeast chromosome 4

For microbiome analysis, total reads corresponding to 10 samples (five each for stool and cecum) ranged from 57,252 to 85,958 (Table 1). Total read numbers corresponding to identified bacteria were not different between the stool and cecal samples (p = 0.345), but a greater number of T. aquaticus reads were observed in the stool samples than in cecal samples (p = 0.008). The average relative abundances of T. aquaticus in the stool and cecal samples were 0.95% and 0.33%, respectively indicating significant differences between the two sample groups (Table 2 and Figure 2).

TABLE 1

Reads of Thermus aquaticus in the samples

	Total reads	Reads of Thermus aquaticus	Relative abundance of Thermus aquaticus
Stool 1	80,157	768	0.96%
Stool 2	76,305	728	0.95%
Stool 3	67,976	554	0.81%
Stool 4	85,958	1,068	1.28%
Stool 5	65,856	499	0.76%
Cecum 1	57,252	105	0.18%
Cecum 2	58,837	59	0.10%
Cecum 3	77,891	468	0.60%
Cecum 4	78,261	324	0.41%
Cecum 5	67,136	246	0.37%

TABLE 2

Comparison of relative abundance of Thermus aquaticus in stool and cecum

	Stools	Cecum
Average	0.95%	0.33%
Standard deviation	0.18%	0.18%
95% confidence interval	[0.79%, 1.11%]	[0.18%, 0.49%]
Median	0.95%	0.37%
p value	0.009**

Wilcoxon rank‐sum test was performed.

FIGURE 2

Relative abundance of Thermus aquaticus composition in fecal (n = 5) and cecal (n = 5) samples. The boxplots indicate the minimum, first quartile, second quartile (median), third quartile, and maximum values. Wilcoxon's rank‐sum test was used (** indicates p‐value <0.01)

As the same amount of yeast harboring the T. aquaticus 16S rRNA gene was added to each of the samples before DNA extraction, the total amount of bacteria in the original samples can be estimated. When X is the relative abundance of T. aquaticus in the stool and Y is the relative abundance of T. aquaticus in the cecum, the formula to estimate the ratio of the number of bacteria in the cecum to that in the stool is as follows:

In other words, on average a cecum sample has 2.9 times more bacteria than a stool sample.

For individual samples, the relative abundances of bacteria belonging to the Bacteroidaceae family were very similar between stool 1 and cecum 2 samples, 37.50%, and 36.86%, respectively (Figure 3). However, the relative abundances of bacteria belonging to the Thermaceae family (the family name of T. aquaticus) were 0.96% and 0.1% in the stool 1 and cecum 2 samples, respectively. In the following formula, the cecum 2 sample was found to have 9.68 times more bacteria than stool 1. Therefore, it is estimated that the number of bacteria belonging to the Bacteroidaceae family in cecum 2 is 9.52 times greater than that in stool 1 (9.68 × 36.86%/37.50%). (X = the relative abundance of Thermaceae in stool 1, Y = the relative abundance of T. aquaticus in cecum 2).

FIGURE 3

The composition of stool and cecal microbiome at the family level. The arrows are the values used as examples in the result

Likewise, the ratios of absolute bacterial abundance from stool 1 to caecum 5 were calculated as 1, 1.00, 1.18, 0.77, 1.27, 5.26, 9.68, 1.60, 2.33, and 2.63 using the relative abundance of Thermaceae in stool 1 as the reference.

The number of operational taxonomic units (OTUs) and the Shannon index were significantly higher in cecum samples (Figure 4a, b), suggesting that the cecal microbiome is more rich and diverse than the stool microbiome. Principal coordinates analysis (PCoA) showed that the samples from the same group were clustered closely in the plot. Permutational multivariate analysis of variance (PERMANOVA) confirmed the significant differences between the microbiome composition of the stool and cecum samples (p = 0.01) (Figure 5).

FIGURE 4

Alpha diversity of the microbiome in the stool and cecum of mice. (a) The number of OTUs (richness) and (b) Shannon index. The boxplots indicate the minimum, first quartile, second quartile (median), third quartile, and maximum values. Wilcoxon's rank‐sum test was used (* indicates p‐value <0.05)

FIGURE 5

Principal coordinates analysis (PCoA) and permutational multivariate analysis of variance (PERMANOVA) of the microbiome in mice stool and cecum. Each dot represents a different sample. Red dots are stool samples and green dots are caecum samples

LEfSe analysis was performed to identify the bacterial species that are differentially enriched in either of the two sample groups. Mucispirillum schaedleri had the highest linear discriminant analysis (LDA) score, 4.63, and its average relative abundance was 8.97 in the stool and 0.43 in the cecum (Table A1). The LDA score of T. aquaticus (3.49) was the third‐highest among all bacterial species (Table A1).

TABLE A1

Linear discriminant analysis effect size (LEfSe) analysis of differentially abundant bacterial taxa between stool and cecum samples. Only taxa (species) meeting an LDA significant threshold of >2 are shown

Taxon name	Taxonomy	LDA effect size	Relative abundance in Stool	Relative abundance in Cecum
Mucispirillum schaedleri	Deferribacteres: Deferribacteres_c: Deferribacterales: Deferribacteraceae: Mucispirillum	4.63034	8.96960	0.43160
KE159600_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159600_g	3.59561	0.11520	0.90320
Thermus aquaticus	Deinococcus‐Thermus: Deinococci: Thermales: Thermaceae: Thermus	3.48987	0.95320	0.33560
KE159605_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159605_g	3.39010	0.04400	0.53480
AB626924_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001228_g	3.35775	0.47080	0.01520
PAC001360_g_uc	Firmicutes: Clostridia: Clostridiales: Christensenellaceae: PAC001360_g	3.28489	0.21920	0.60440
PAC001125_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Kineothrix	3.19065	0.34520	0.03520
PAC001757_s group	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g21	3.14647	0.03840	0.31840
PAC001782_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: PAC000661_g	3.13769	0.30400	0.02960
PAC001104_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC000664_g	3.09174	0.04800	0.29480
PAC002512_s group	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g21	3.00307	0.02440	0.22560
PAC001294_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g21	2.99883	0.06840	0.26760
PAC001706_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Oscillibacter	2.99829	0.22080	0.02200
PAC001092_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001092_g	2.95863	0.01600	0.19760
PAC001727_s group	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001385_g	2.92810	0.01320	0.18240
PAC001684_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Oscillibacter	2.82737	0.16000	0.02600
PAC001638_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Oscillibacter	2.81236	0.17840	0.04880
PAC001557_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001043_g	2.75242	0.01120	0.12400
PAC001374_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g24	2.72112	0.01400	0.11880
PAC001361_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001103_g	2.71064	0.00600	0.10840
KE159781_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g21	2.68962	0.04400	0.14160
PAC001382_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001092_g	2.66791	0.01640	0.10920
PAC001772_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g21	2.65644	0.01600	0.10640
KE159714_s group	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Oscillibacter	2.64527	0.11760	0.02960
PAC001542_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Pseudoflavonifractor	2.63569	0.11080	0.02480
PAC002401_s group	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159605_g	2.62727	0.01320	0.09760
PAC001668_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001092_g	2.61322	0.01800	0.09960
PAC000668_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Monoglobus	2.60557	0.00120	0.08160
PAC001744_s group	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Frisingicoccus	2.60481	0.05280	0.13280
AB606336_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Agathobaculum	2.60178	0.02840	0.10800
Adlercreutzia equolifaciens group	ActinoCoriobacteriia: Coriobacteriales: Coriobacteriaceae: Adlercreutzia	2.59700	0.03080	0.10960
Flintibacter butyricus group	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Pseudoflavonifractor	2.59454	0.08120	0.15960
Clostridium_g24_uc	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g24	2.58135	0.00320	0.07920
PAC001088_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Pseudoflavonifractor	2.58021	0.09440	0.01880
PAC002042_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC002042_g	2.51002	0.01200	0.07640
PAC001476_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Agathobaculum	2.50193	0.01120	0.07440
PAC002397_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Pseudoflavonifractor	2.45536	0.06200	0.00520
PAC000692_g_uc	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC000692_g	2.42353	0.00600	0.05880
PAC001535_s group	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC000664_g	2.42156	0.07000	0.01760
JX095379_g_uc	Firmicutes: Clostridia: Clostridiales: Christensenellaceae: JX095379_g	2.41006	0.05840	0.00720
PAC001360_s	Firmicutes: Clostridia: Clostridiales: Christensenellaceae: PAC001360_g	2.38939	0.06560	0.01680
PAC001087_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Pseudoflavonifractor	2.35129	0.08960	0.04520
PAC002147_s	Firmicutes: Clostridia: Clostridiales: Christensenellaceae: PAC002147_g	2.32864	0.04320	0.00080
EF098562_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159600_g	2.32697	0.00360	0.04560
PAC002354_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Acutalibacter	2.32197	0.01360	0.05480
PAC002159_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: PAC000661_g	2.32060	0.04680	0.00520
KE159600_g_uc	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159600_g	2.30840	0.00000	0.04040
PAC001096_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Sporobacter	2.29992	0.04720	0.00760
PAC001092_g_uc	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001092_g	2.25716	0.00920	0.04440
PAC001043_g_uc	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001043_g	2.25524	0.00080	0.03640
PAC002555_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: PAC002555_g	2.23590	0.03480	0.00080
PAC001536_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: PAC000661_g	2.23329	0.04240	0.00840
EU511112_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159600_g	2.22109	0.01200	0.04400
PAC001767_s	Firmicutes: Clostridia: Clostridiales: Mogibacterium_f: PAC001236_g	2.20526	0.01760	0.04920
PAC001396_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159605_g	2.18935	0.01000	0.04040
PAC000184_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g24	2.18879	0.00160	0.03160
PAC002471_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC002471_g	2.17289	0.03320	0.00400
AB626948_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Oscillibacter	2.12901	0.03760	0.01200
PAC001551_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Oscillibacter	2.12092	0.04480	0.01880
PAC001473_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC000671_g	2.10215	0.00360	0.02800
AF371672_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: LLKB_g	2.07844	0.00000	0.02360
PAC001222_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: KE159797_g	2.07604	0.00200	0.02560
PAC001543_s	Firmicutes: Clostridia: Clostridiales: Ruminococcaceae: Pseudoflavonifractor	2.05036	0.01680	0.03840
PAC001085_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Clostridium_g21	2.04864	0.00560	0.02720
PAC001681_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001681_g	2.04810	0.01120	0.03320
PAC002350_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001372_g	2.04045	0.00440	0.02520
PAC001722_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: Coprococcus_g2	2.04004	0.00400	0.02560
PAC001296_s	Firmicutes: Clostridia: Clostridiales: Lachnospiraceae: PAC001296_g	2.00101	0.00000	0.00120

To validate this technique, we analyzed twofold serial dilutions of the stool suspension. An equal amount of yeast was added to the diluted samples (x1, x1/2, x1/4 and x1/8); the DNA was extracted and HTS was performed. The results showed that the relative abundances of the added T. aquaticus sequence in the twofold serially diluted samples were 0.91%, 1.77%, 3.72%, and 6.97%, respectively (Figure 6). The ratio of the estimated absolute abundances of total bacteria in the diluted samples (x1, x1/2, x1/4 and x1/8) was calculated as 8.15: 4.20: 1.96: 1.

FIGURE 6

The microbial composition of the serially diluted (twofold) stool suspension at the family level

DISCUSSION

HTS of 16S rRNA gene amplicon provides information on the relative abundances of bacteria (compositional data), but not the absolute abundances. In compositional data, as one taxon increases within a sample, other taxa must decrease. The compositional data alone may not reflect the true microbiome of the sample thus resulting in biased statistical results when analyzed in an inappropriate manner (Gloor et al., 2017). To address this problem, we developed a method to measure the absolute abundances of bacterial populations in stool and cecal samples by adding Pichia pastoris that had been genetically modified to contain the 16S rRNA gene from a hyperthermophile.

We selected the hyperthermophile T. aquaticus as an internal standard because it is rarely present in clinical and environmental samples. In addition, given that T. aquaticus is commonly used in research, it was anticipated that every bacterial database will have information on this bacterium. Indeed, the 16S rRNA gene databases, such as EzBioCloud, SILVA, Greengenes and NCBI had information regarding the nucleotide sequence of the T. aquaticus 16S rRNA gene. In the beginning, we had tried to use E. coli transformed with pUC57 harboring the T. aquaticus 16S rRNA gene directly, but the ribosomal DNA of E. coli, the host of plasmid, was constantly detected in the tests. Therefore, we used P. pastoris (yeast) as the host for stable amplification of the T. aquaticus 16S rRNA gene, without any E. coli genome contamination.

There have been numerous studies to estimate absolute abundances by quantitative real‐time PCR or cell counting using microscopy or flow cytometry. Compared to using an internal standard, these methods are considered to be costly, labor‐intensive, time‐consuming and require specialized equipment and skills (Harrison et al., 2021). In addition, the quantitative PCR method is less accurate in calculating the converted absolute abundance than the internal standard method (Stämmler et al., 2016).

Recently, studies have been conducted to estimate absolute abundances in microbiota by incorporating internal standards into samples. There are two types of internal standards for metabarcoding: cellular and DNA. Various cellular internal standards have been studied such as Escherichia coli, Salinibacter ruber, Rhizobium radiobacter, Alicyclobacillus acidiphilus, Sporosarcina pasteurii, and Shewanella oneidensis (Ji et al., 2019; Jones et al., 2015; Piwosz et al., 2018; Stämmler et al., 2016). The ideal cellular internal standard should have similar traits to focal organisms, be easily cultured, not occur in the biological samples and not have copy number variation (Harrison et al., 2021). Genomic DNA can be used as internal standards (Deagle et al., 2018; Ji et al., 2020; Lin et al., 2019; Smets et al., 2016; Venkataraman et al., 2018). Synthetic DNA has also been utilized as an internal standard (Hardwick et al., 2018; Jiang et al., 2011; Tkacz et al., 2018; Tourlousse et al., 2017; Zemb et al., 2020). Synthetic sequences do not occur naturally in samples and the problem with copy number variation can be eliminated.

The method we developed is a powerful approach that maximizes the strengths of both cellular and DNA internal standards, as it uses cells (yeasts) that contain the DNA sequence of the 16S rRNA gene of T. aquaticus. The 16S rRNA gene of T. aquaticus, a hyperthermophile, is a sequence suitable for metabarcoding because it is registered in all databases but is not usually detected in animals or most environments. The main benefit of our approach is that yeast (eukaryote) lacks a 16S rRNA gene and will not interfere with the analysis. Yeasts are easier to cultivate than hyperthermophiles or halophiles. Furthermore, the copy number of the T. aquaticus 16S rRNA gene in all yeast cells is the same over time because all yeast cells used in this technique originated from a single recombination clone. We confirmed that a single copy of the 16S rRNA gene was inserted into the yeast genome and the 16S rRNA gene amount per yeast arginosuccinase gene was unchanged over three consecutive rounds of yeast cultivation (Figure A1a).

FIGURE A1

(a) The relative quantity (RQ) of the 16S rRNA gene and (b) colony‐forming units (CFU) of the recombinant yeast over three consecutive rounds of yeast cultivation (n = 5 for each day)

Using an existing 16S rRNA gene rather than a synthetic sequence probably makes it easier to analyze metagenomic bioinformatics in ordinary laboratories. Since the bacterial 16S rRNA gene (such as the 16S rRNA gene of T. aquaticus) is registered in most databases, it can be used as an internal standard in any preexisting pipeline without methodological modification or even be directly applied to a web‐based metagenomics pipeline. Although hyperthermophiles are believed to exist only in extreme environments, a few studies have reported the presence of Thermus spp. in some gut and lung microbiomes (Jones et al., 2018; Yu et al., 2016). However, we believe that T. aquaticus is very suitable if the bacterial 16S rRNA gene should be used as an internal standard. Other studies have also studied the use of the 16S rRNA gene of Thermus spp. as an internal standard (Janes et al., 2020; Lin et al., 2019; Smets et al., 2016). However, when Thermus sp. is present in the sample to be studied, using a synthetic sequence as an internal standard can be a solution. In this study, no Thermus spp. other than the internal standard T. aquaticus were found. The data in Figures 3 and 6 are presented at the level of the family for visual convenience. In these figures, Thermaceae is composed exclusively of T. aquaticus.

In addition, it may be difficult to add the same number of newly cultured yeast cells to a sample, as in previous experiments. The yeast we developed did not show the same cell count per culture batch. (Figure A1b). Therefore, we recommend cultivating large quantities of yeast, making aliquots, storing them and using them for each experiment that requires an internal standard.

In this study, the average relative abundances of T. aquaticus 16S rRNA gene in five stool samples were 0.95% and at 95% confidence interval, [0.79%, 1.11%] the fluctuating abundances centered on 1% (Table 2). In a pilot study using quantitative PCR, we defined that one unit of yeast to prepare the T. aquaticus 16S rRNA gene could represent 1% of the total microbial 16S rRNA gene in one stool sample. When 50,000 total bacterial reads are obtained in one sample, 500 reads of T. aquaticus are required to reach 1% of the total reads. To calculate the differences between the samples while minimizing the effects on the remaining 99% of bacteria, 500 reads (1%) were believed to be sufficient. However, because one unit is set as 1% in one stool sample, it may be necessary to adjust the unit when studying samples other than stools. When the DNA density of the sample is too low, it might be effective to use artificial synthetic sequences or qPCR to measure the absolute abundances of the samples.

The study provided the necessary formula for comparing the absolute abundances of bacteria between the two sample groups, stool, and cecum, that is, five stool samples and five cecal samples. The average relative abundances of T. aquaticus in stool and cecal samples were calculated to 0.95% and 0.33%, respectively. The value obtained using the proposed formula revealed that a cecum sample contains 2.9 times more bacteria than a stool sample. This is likely due to the difference in the original amount of samples between the two groups.

In addition, the yeast we developed and the calculation used showed very accurate estimation ability in validation experiments using serial dilutions (Figure 6). It calculated the ratio of the total bacteria in twofold diluted samples (x1, x1/2, x1/4, and x1/8) as 8.15:4.20:1.96:1. A stool sample without internal standards was not included as a control because T. aquaticus is not expected to be present in fecal samples.

When comparing the absolute abundances of a specific bacterial taxon between two samples, the relative abundance of the target bacteria can be multiplied by the value derived using the formula. In the two sample groups, the relative abundances of bacteria belonging to the Bacteroidaceae family were similar, but the absolute abundance was 9.52 times higher in the cecum 2 sample than that in the stool 1 sample. However, it should be noted that the resulting values from our research method cannot be expressed in an exact unit (e.g., copies per g) and the efficiency of cell lysis for yeast may be different than that for bacteria, so the values have meaning mostly in a comparative sense.

In addition, the number of OTUs (richness) and Shannon index (alpha diversity analysis) revealed that cecal samples exhibit greater microbial diversity than the stool samples. PCoA and PERMANOVA confirmed the differences in microbial composition between the two sample groups. The observations of the current study are concordant with previous studies (Gu et al., 2013; Tanca et al., 2017).

However, the effect of adding the T. aquaticus 16S rRNA gene on the microbiome analyses, such as alpha and beta diversity cannot be excluded. This limitation can be overcome in various ways. One way is to add a trace amount of the T. aquaticus 16S rRNA gene to minimize its effect on the analysis of other bacterial taxa. The next method is to subtract the T. aquaticus reads from the taxonomy output before proceeding with the microbiome analysis. Lastly, HTS can be performed using the same samples with and without adding the T. aquaticus 16S rRNA gene.

CONCLUSION

In this study, we developed a method to measure the absolute abundances of the microbial community when performing HTS and 16S rRNA gene amplicon analysis by adding P. pastoris harboring T. aquaticus 16S rRNA gene to the stool and cecal samples of mice.

CONFLICT OF INTEREST

None declared.

AUTHOR CONTRIBUTIONS

Ju Yeong Kim: Conceptualization (lead); Writing‐original draft (lead). Myung‐hee Yi: Investigation (lead); Methodology (lead). Myungjun Kim: Validation (lead). Seogwon Lee: Data curation (lead). Hye Su Moon: Resources (lead). Dongeun Yong: Formal analysis (lead). Tai‐Soon Yong: Funding acquisition (lead); Supervision (lead).

ETHICS STATEMENT

None required.