Quantitative real-time PCR analysis of bacterial biomarkers enable fast and accurate monitoring in inflammatory bowel disease

Sample collection

Fecal and blood samples were collected from 18 IBD patients (six diagnosed with ulcerative colitis and 12 diagnosed with Crohn’s disease) and four healthy (control group) individuals in the Gastroenterology Department of Dokuz Eylul University Hospital (Izmir, Turkey). IBD patients (ulcerative colitis and Crohn’s disease) were diagnosed according to international guidelines based on clinical, endoscopic, histopathological, and radiological examinations (Stange et al., 2008; Bemelman et al., 2018). The healthy controls were candidates without any history of IBD or mucosal lesions in colonoscopy. Fecal samples gathered in sterile and airtight containers, and blood samples collected into EDTA tubes were transported to laboratory within six hours after collection. All study samples were kept at −80 °C until processing. Ethical approval was obtained from the Ethics Committee of the Dokuz Eylül University (2017/08-03). All participants provided informed consent in the format required by the Dokuz Eylül University ethics committee.

Genomic DNA extraction from blood samples

DNA was isolated from blood samples with the Genomic DNA Mini Kit (blood/cultured cell) (Geneaid Biotech Ltd., Taiwan) following the manufacturer’s protocol. The quality of extracted DNAs (A260/A280 and A260/A230 ratios) was checked with a Nanodrop 8000c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

Primer design for candidate variants in IBD associated candidate genes

We targeted rs2066844 (Arg702Trp), rs2066845 (Gly908Arg), rs2066847 (Leu1007insC) SNPs for the NOD2 gene; rs11209026 (Arg381Gln) for the IL-23R gene, and rs2241880 (Thr300Ala) for the ATG16L1 gene. DNA sequences for each gene were obtained from NCBI (https://www.ncbi.nlm.nih.gov/), ENSEMBL (https://www.ensembl.org/index.html), and UCSC Genome Browser (https://genome.ucsc.edu/cgi-bin/hgGateway) databases. These gene sequences were used to design polymerase chain reaction (PCR) primers using the IDT SciTools PCR algorithm (Integrated DNA Technologies) (Table S1). Oligonucleotide properties, melting temperature, hairpins, dimers, and mismatches were identified by IDT SciTools OligoAnalyzer 3.1 (https://www.idtdna.com/calc/analyzer) software (Owczarzy et al., 2008), and specificity of primers were confirmed with Primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast/).

Amplification and sequencing of IL-23R, ATG16L1 and NOD2 variants

With primers designed for the specified SNPs, PCR analysis was performed using the FastStart High Fidelity PCR System, dNTPack kit (Roche Applied Science, Penzburg, Germany). Reaction mixes were made separately for IL-23R, NOD2, and ATG16L1 genes in a final volume 25 µl of using 17.25 µl PCR-grade water, 0.5 µl forward and reverse primers, 0.5 µl PCR Grade Nucleotide Mix, 2.5 µl FastStart High Fidelity Reaction Buffer, 0.5 µl dimethyl sulfoxide (DMSO), 0.25 µl FastStart High Fidelity Enzyme Blend, and 3 µl DNA. The thermal cycling was subjected to the following conditions: denaturation at 94 °C for 10 min followed by 35 cycles of 94 °C for 2 min, annealing at 57 °C for 30 s, and elongation at 72 °C for 1 min using SimpliAmp Thermal Cycler (ThermoFisher Scientific, Waltham, MA, USA).

PCR products were verified by agarose gel electrophoresis. Briefly, 5 µL of PCR product was mixed 1 µL 6X DNA loading buffer, and run on 1.4% agarose gel in 0.5X TBE buffer under a steady voltage of 100 V for 60 min at room temperature.

PCR products were purified using ExoSAP-IT™ PCR Product Cleanup Reagent (Thermo Fisher Scientific, Waltham, MA, USA). 10 µl sequencing mixture contained 4 µl ddH₂O, 1 µl 5X ABI Buffer, 1 µl of primer (3.2 pmol/µl), 2 µl BigDye Ready Reaction Mix, and 2 µl of the PCR product. Samples were sequenced on ABI Prism 3130xl Genetic Analyzer (Applied Biosystems, Waltham, MA, USA). Sequencing results (ABI chromatograms) were analyzed in Unipro UGENE v.33 (Okonechnikov et al., 2012) program. Multiple sequence alignments using ClustalW algorithm were performed in Unipro UGENE v.33 program (Okonechnikov et al., 2012).

Bacterial DNA isolation from stool samples

DNA was extracted from stool samples using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer’s protocol. The concentration (ng/µL) and purity (A260/A280 and A260/A230 ratios) of the DNA samples were determined by Nanodrop 8000c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA).

16S rRNA gene amplicon sequencing by next generation sequencing (NGS)

Nine stool samples (three samples from each UC, CD patients, and healthy volunteers) were chosen for amplicon analysis. The variable V3-V4 region of 16S rRNA gene was targeted and amplified with the following PCR primers: 341F (5′-CCTAYGGGRBGCASCAG-3′) and 806R (5′-GGACTACNNGGGTATCTAAT-3′). After the purification of PCR products, sequencing libraries were generated with Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA, USA). The concentration of sequencing libraries are standardized to 4nM each. Normalized samples were pooled and sequenced by Illumina NovaSeq 6000 as paired-end (2 × 250 bp) using the manufacturer’s standard procedure. Raw data quality control check was performed by FastQC, and quality control of the reads was checked by QIIME2 (Caporaso et al., 2010). Effective tags were obtained after removing primer and barcode sequences, chimeric reads, and reads with Phred Score less than 20 by DADA2 (Callahan et al., 2016). By utilizing the effective tags, representative sequence for each Operational Taxonomic Units (OTUs) were acquired with ≥97% similarity against the Greengenes and SILVA databases. QIIME2 was used for taxonomic determination of each OTU. Rarefaction curves plotting sequencing depth vs. number of taxa identified were used to judge the appropriateness of sequencing depth for each sample (Pereira-Marques et al., 2019; Zaheer et al., 2018).

Species diversity within samples were assessed by five different Alpha diversity estimates including observed-species, Chao1, Shannon, Simpson, and ACE indices. Alpha diversity indices were visualized via boxplots. All downstream analyzes were performed with “phyloseq” (McMurdie & Holmes, 2013) and “ggplot2” packages in R software (Version 4.0.5) (https://www.r-project.org).

Quantification of selected bacterial species levels using Real-Time quantitative PCR (qPCR) analysis

Six bacterial species (Faecalibacterium prausnitzii, Clostridioides difficile, Akkermansia muciniphila, Bacteroides vulgatus, Streptococcus thermophilus, and Shiga toxin-producing Escherichia coli (stx1 gene positive)) were chosen as bacterial biomarkers for further quantification via qPCR in all samples. Primers were designed using the IDT SciTools^® platform (Owczarzy et al., 2008). Primer sequences for each species and total bacterial quantification are shown in Table 1. In each run both positive and negative controls were added to qPCR-plates. For positive control A. muciniphila (ATCC BAA-835), B. vulgatus (ATCC 8482), S. thermophilus (ATCC 19258), F. prausnitzii (ATCC 27766), E. coli (ATCC 43890), and C. difficile (ATCC 9689) were utilized. A total of 2.5 µl of distilled water was used as a negative control in each run.

The reactions were conducted with 2.5 µl DNA, 1.9 µl PCR-grade water, 0.3 µl primer, and 5 µl LightCycler^® 480 SYBR Green I Master enzyme (Roche Applied Science, Penzburg, Germany) in a total volume of 10 µl. All reactions were carried out with triple replicates on a LightCycler^® 480 II (Roche Applied Science, Penzburg, Germany) qPCR machine. Each sample and species is assayed at least twice. The qPCR reaction was 95 °C for 10 min with initial denaturation followed by 50 cycles of 95 °C for 10 s, 58 °C for 15 s, and 72 °C for 15 s. Melting curve analysis for the qPCR products was performed under the following conditions: 95 °C for 5 s, 63 °C for 1 min and a denaturing temperature ramp from 63 to 97 °C with a rate of 0.11 °C/s. Amplification and melting curves for each sample were obtained using Absolute Quantitation/Second Derivative and Tm Calling analysis modes in the LightCycler^® 480 II Software v.1.5.

SYBR Green dye fluorescence intensity was used for quantification. The target bacterial DNA concentration correlated with the threshold cycle number (Ct), the cycle number at which fluorescence signal was first detected. Roche LightCycler^® 480 System melting curve program analysis is used for the confirmation of success of qPCR reactions (Simenc & Potocnik, 2011). Data analysis was conducted using both the 2^−ΔCt, and 2^−ΔΔCt methods (Livak & Schmittgen, 2001). Target microorganisms were considered as a target while total bacteria measurement was used as a reference (Navidshad, Liang & Jahromi, 2012). Relative abundance of bacteria was expressed as log2 transformed fold change values, and calculated according to the following formulas.

Relative abundance of target bacteria species with respect to abundance of total bacteria: ${\displaystyle 2^{-\Delta \text{Ct}}=2^{\text{- (Ct of target bacteria}-\text{Ct of total bacteria)}}.}$ Fold change of relative abundance of target bacteria in IBD patients compared to healthy controls ${\displaystyle 2^{-\Delta \Delta \mathrm{Ct}}=2^{\text{- [(Ct of target bacteria}-\text{Ct of total bacteria) patient - (Ct of target bacteria}-\text{Ct of total bacteria) control]}}.}$ To verify PCR efficiency, standard curves were generated by 10-fold dilutions of bacterial DNA for all primer sets. In all sets, qPCR efficiency was >90% and calculated by E= 10^{(−1/slope)−1} equation. According to the serial dilutions, the limit of detection of qPCR assays was 10–100 copies.

Statistical analysis

Distribution of candidate gene’s genotypes among UC, CD, and control groups was compared with Chi-square categorical analyses. Shapiro–Wilk test was used to test the normality assumption of numeric variables. All microbial diversity and quantity parameter estimates violated normality assumption, so groups were compared using the non-parametric Kruskal–Wallis and pairwise Wilcoxon rank-sum tests. P-value less than 0.05 was assessed as statistically significant. All statistical analyzes were performed with R software (Version 4.0.5) (https://www.r-project.org). We also estimated the statistical power to detect a 10% change with 95% confidence in the numeric abundances of A. muciniphila, B. vulgatus, S. thermophilus, and F. prausnitzii in UC and CD groups based on qPCR relative abundance estimates. Statistical power calculations followed formulations in Cohen (1988) as implemented in the pwr package (https://cran.r-project.org/web/packages/pwr/index.html) of R software. Experimentally determined relative abundance means and standard errors of A. muciniphila, B. vulgatus, S. thermophilus, and F. prausnitzii are used for statistical power calculations, where Type I error probability (alpha) is set to 0.05.

Patient characteristics

UC and CD patients consisted of similar age and sex groups, however, included a variety of disease locations and disease activity phenotypes (Table 2, Table S2). Distribution of genotypic frequencies of ATG16L1, IL23R and NOD2 variants in UC and CD patients was similar (Table S3), suggesting similar IBD genetic risk profile in these patient groups.

Identification of fecal microbiota profile in the IBD patients and controls

An average of 148,545 reads per sample (range 125,105–172,435) were generated for the nine fecal samples (three samples from each of CD, UC, and control groups). Eight phyla, 15 classes, 23 orders, 43 families, 96 genus, and 233 species were represented in all sequences based on a 97% similarity level. The most abundant phyla in CD samples were Firmicutes (36.22%), Proteobacteria (29.22%), Verrucomicrobia (25.79%), Bacteroidetes (8.68%), Actinobacteria (0.07%), Fusobacteria (0.01%), and Synergistetes (0.01%). Similarly in UC samples, Proteobacteria (45.63%), Firmicutes (29.21%), Bacteroidetes (10.14%), Fusobacteria (9.98), Actinobacteria (2.80%), and Verrucomicrobia (2.25%) were observed (Fig. 1A). The most common phyla in the control samples were Bacteroidetes (62.63%), followed by Firmicutes (31.47%), Proteobacteria (4.01%), Actinobacteria (1.81%), Fusobacteria (0.04%), Lentisphaerae (0.04%), and Verrucomicrobia (0.01%) (Fig. 1A). Despite small sample size Kruskal–Wallis tests showed abundance difference for Proteobacteria (p = 0.03), and Firmicutes (p = 0.07) between IBD and control samples. At the family level, significant abundance differences for Bacteroidaceae (p = 0.006), Rikenellaceae (p = 0.03), Acidaminococcaceae (p = 0.05), Victivallaceae (p = 0.03), and Enterobacteriaceae (p = 0.03) were observed between IBD and control samples (Fig. 1B).

Reduction in alpha diversity estimates in the IBD samples (UC and CD groups) compared to control samples was observed (Fig. S1), however only the Shannon diversity estimate was significantly lower in the CD group compared to the control group (p = 0.04). Overall microbial dysbiosis commonly reported in IBD, and greater dysbiosis in CD is captured in the study.

Identification and quantification of biomarker bacterial species

Fecal microbiota analyses identified not only the taxa unique for UC, CD, and control groups, but also common bacterial species among the groups (Table S4). In fact, 29 bacterial species that were found in all three groups cumulatively constituted nearly 60% of all taxon assigned reads identified in the fecal microbiota (Fig. 1C). Among the 29 bacterial species found in all three groups, Faecalibacterium prausnitzii, Akkermansia muciniphila, Bacteroides vulgatus, Streptococcus thermophilus, and Escherichia coli were the most common ones making up 84% of the reads assigned to these 29 bacterial species.

Shiga toxin-producing Escherichia coli is a proinflammatory bacteria, reported to show increasing abundance in IBD patients, whereas harmful or beneficial association of Bacteroides vulgatus in IBD is less certain (da Silva Santos et al., 2015; Palmela et al., 2018; Zafar & Saier Jr, 2021). Faecalibacterium prausnitzii, Streptococcus thermophilus, and Akkermansia muciniphila are beneficial bacteria, reported to be reduced in IBD patients (Sokol et al., 2008; Prosberg et al., 2016; Zafar & Saier Jr, 2021). These five species are chosen as potential biomarker bacteria for further quantification in all UC and CD samples. Because higher abundance of pathogenic Clostridioides difficile is also reported to be associated with IBD (Issa et al., 2007), C. difficile is also chosen to be further quantified in the IBD samples.

Quantitative real-time PCR analyses were conducted with primers specific to each target species (Table 1). Specificity of primers and success of the qPCR reactions were checked by melting curve analyses (Fig. S2). After quantification of six selected species in all twenty two samples by qPCR, the relative abundance of candidate species compared to total bacteria amount is calculated by the 2^−ΔCt method for each sample. Reduction in the relative abundance of beneficial species F. prausnitzii and A. muciniphila, in UC and CD samples compared to total bacteria amount was evident (Fig. 2). The relative abundance of B. vulgatus was higher in the control samples (Fig. 2). Interestingly, the relative abundance of S. thermophilus was highest in the UC samples (Fig. 2). Shiga toxin-producing E. coli and C. difficile was observed in only three of the CD patients. The next analyses compared the fold change of A. muciniphila, B. vulgatus, S. thermophilus, and F. prausnitzii in IBD samples with respect to control samples calculated by the 2^−ΔΔCt method. Significant reduction in A. muciniphila, S. thermophilus, and F. prausnitzii in combined IBD samples compared to controls was observed (Fig. 3). However, fold change values of these four species were similar in the UC and CD samples (Fig. 3).

Effect of disease course and treatment on common bacterial species is also examined. Patients on biologics treatment presented with higher reduction in the relative abundance of beneficial F. prausnitzii compared to patients who are not on biologics (mean log reduction 5.8 vs. 2.1, p = 0.04). Within the CD group, patients who had surgery showed higher reduction in the relative abundance of A. muciniphila compared to patients who did not have surgery (mean log reduction 21 vs. 15, p = 0.02). Moreover, patients with penetrating perianal CD (more severe version of CD) again showed higher reduction in the relative abundance of A. muciniphila (mean log reduction 21 vs. 16, p = 0.05), and increase in the relative abundance of B. vulgatus (mean log increase 10 vs. 5, p = 0.05). Within the UC group, patients with CRP levels higher than five showed higher reduction in the relative abundance of B. vulgatus compared to patients with CRP levels less than five (mean log reduction 14 vs. 4, p = 0.05).

Utility of NGS (Next Generation Sequencing) and qPCR results

None of the correlations between the NGS read and the qPCR results of either individual selected bacteria or total bacteria was statistically significant. Highest correlation was observed for total bacteria results (Adjusted R² = 0.14) but it was not significant (p = 0.18). Lack of congruency between the two techniques was in part due to great variation in NGS read results. The standard deviation of NGS reads ranged from a lowest of 724 for A. muciniphila to 4,352 for S. thermophilus. Contrary, the standard deviation of qPCR results ranged from a lowest of 6.3 for A. muciniphila to 14.4 for B. vulgatus, showing orders of magnitude smaller variance in the qPCR results compared to NGS results. However, the lower variation in qPCR assays is also due to the concomitant characteristics of the qPCR technique and range of the data generated. Moreover, we did not incorporate unique molecular identifier (UMI) analysis in NGS data that enables identification of PCR duplicates into sequencing libraries, which can introduce noise in the read counts.

The much lower variation inherent in qPCR results suggests higher accuracy, precision, and higher clinical utility for the qPCR technique compared to NGS based amplicon sequencing approach especially with smaller sample sizes. Statistical power calculations based on mean and standard deviation estimates from our experimental data showed that qPCR results have over 80% statistical power to detect a minimum 10% difference in candidate species abundance with less than 200 samples (Fig. 4). Such high statistical power is nearly impossible to achieve with NGS methods. Given that IBD clinics have hundreds of patients, and microbiota alterations as a response to therapy/interventions are usually much higher than 10%, qPCR surveillance of candidate bacteria has good clinical potential to monitor microbiota response through disease and therapy course.