Association of Flavonifractor plautii, a Flavonoid-Degrading Bacterium, with the Gut Microbiome of Colorectal Cancer Patients in India

Cohort design and subject enrollment.A considerable sample size consisting of 60 samples (30 cases and 30 controls) was recruited from two different locations (Bhopal and Kerala) in India. For constructing a balanced data set, 15 cases and 15 controls were selected from both the locations. The two selected locations represent different geographies (2,000 km apart) and lifestyles in order to remove the confounding effect of diet and making the observations generalizable for the Indian cohort. Bhopal is a city located in central India and is populated with people from all over the country; hence, samples from here can act as a proxy to represent the diversity of the country. Samples from Kerala were specifically chosen because, among all the other states of India, Kerala has the highest rate of colorectal cancer incidence. The fecal samples were collected only from CRC cases, and those from healthy subjects were taken from a previous study (12). Each fecal sample was collected and immediately transported to the lab at 4°C for further processing. Diagnosis of all the cases was carried out by experienced oncologists at the hospitals through biopsy and colonoscopy. The study exclusion criteria for patients were any previously diagnosed serious medical conditions and recent use of antibiotics, to avoid the effect of confounding factors. Patients with incomplete medical information were also removed from the selected set. Fecal samples were collected prior to colonoscopy in sterile containers.

Fecal metagenomic DNA extraction.Metagenomic DNA was isolated from all the fecal samples using the QIAamp stool minikit (Qiagen, CA, USA) according to the manufacturer’s instructions. DNA concentration was estimated by the Qubit HS double-stranded DNA (dsDNA) assay kit (Invitrogen, CA, USA), and quality was estimated by agarose gel electrophoresis. All the DNA samples were stored at −80°C until sequencing.

Shotgun metagenome sequencing.The extracted metagenomic DNA was used to prepare the sequencing libraries using the Illumina Nextera XT sample preparation kit (Illumina Inc., USA) by following the manufacturer’s protocol. The sizes of all the libraries were assessed on the Agilent 2100 Bioanalyzer using the Agilent high-sensitivity DNA kit (Agilent Technologies, Santa Clara, CA, USA) and were quantified on a Qubit 2.0 fluorometer using the Qubit HS dsDNA kit (Life Technologies, USA) and by quantitative PCR (qPCR) using Kapa SYBR Fast qPCR master mix and Illumina standards and primer premix (Kapa Biosystems, MA, USA) according to the Illumina suggested protocol. The shotgun metagenomic libraries were loaded on an Illumina NextSeq 500 platform using the NextSeq 500/550 v2 sequencing reagent kit (Illumina Inc., USA), and 150-bp paired-end sequencing was performed at the Next-Generation Sequencing (NGS) Facility, IISER, Bhopal, India.

Preprocessing of the metagenomic reads.A total of 150 Gbp of metagenomic sequence data (mean = 1.36 Gb) was generated from 60 fecal samples. The metagenomic reads were filtered using the NGSQC (v2.3.3) toolkit with a cutoff of q of ≥20 (47). The high-quality reads were further filtered to remove the host-origin reads (human contamination) from bacterial metagenomic reads, which resulted in the removal of an average of 1% of reads. The reads from each sample were assembled into contigs at a k-mer size of 63 bp using SOAPdenovo (v2.0) (48). The singletons resulting from each sample were pooled, and de novo assembly was repeated on the combined set of singleton reads from all samples. The open reading frames (ORFs) from each contig (length of ≥300 bp) were predicted using MetaGeneMark (v3.38) (49). Pairwise alignment of genes was performed using BLAT (v2.7.6), and the genes which had an identity of ≥95% and alignment coverage of ≥90% were clustered into a single set of nonredundant genes, from which the longest gene was selected as the representative ORF to construct the nonredundant gene catalogue.

The Integrated Gene Catalogue (IGC), which represents 1,297 human gut metagenomic samples comprising HMP, MetaHIT and Chinese data sets, was retrieved (50). The gene catalogue constructed from Indian samples was combined with the IGC to construct a nonredundant gene catalogue (using ≥95% identity and ≥90% alignment coverage) and is referred to as “updated IGC” in the subsequent analysis.

Quantification of gene content.The quantification of gene content was carried out using the strategy performed by Qin et al. (51) in which the high-quality reads were aligned against the updated IGC using SOAP2 (v2.21) in the SOAP aligner with an identity cutoff of ≥90% (52). Two types of alignments were considered for sequence-based profiling: (i) the entire paired-end read mapped to the gene and (ii) one end of the paired-end read mapped to a gene and other end remained unmapped. In both cases, the mapped read was counted as one copy. Further, the read count was normalized based on length of the gene as b_i = x_i/L_i.

The relative abundance of a gene within the sample was calculated as follows:ai=bi∑ jbj=xiLi∑ jxjLj where a_i is relative abundance of gene in sample S, x_i is number of times that gene i was detected in sample S (the number of mapped reads), L_i is length of gene i, j is all the genes, and b_i is copy number of gene i in sequenced data from sample S.

Construction of updated gene catalogue for gut profiling.To construct the gene catalogue for gut microbiome profiling, the high-quality sequencing reads were subjected to a de Bruijn graph-based assembly which resulted in 2,143,541 contigs of >300 bp in length with a total contig length of 1.52 Gb. To capture low-coverage genomic regions or low-abundance genomes, all unassembled reads were extracted and combined with the singletons from each sample to further assemble into an additional 1.2 million contigs (>300 bp) with a total assembled length of 0.76 Gb. The gene prediction on all assembled contigs resulted in 4,591,809 genes, out of which 2,36,4248 genes were nonredundant and which represents the gene catalogue of the Indian population. We incorporated these genes to update the currently available Integrated Gene Catalogue (IGC), which contained 9.8 million genes from 1,267 gut metagenomes from three continents (Europe [53–55], United States [56], and China [51]), as it lacked information on genes specific to the Indian population. The Updated Gene Catalogue (UGC) comprised 11,118,467 genes (an addition of 12.5% genes in the current IGC) with 1,238,571 genes unique to the Indian population.

On this updated gene catalogue, reads from each sample were mapped and the genes present in the Indian population were identified. On average, 54.47% ± 7.84% (average ± SD) of high-quality reads mapped from each sample to UGC and resulted in the identification of 3.8 million genes present in the Indian cohort. Taxonomic assignment and functional annotation were performed for these 3.8 million genes present using 4,097 reference genomes (HMP and NCBI species) and KEGG and eggNOG databases. A total of 2.41 million genes (62.9%) could be successfully assigned a taxonomy at genus level. The remaining genes are expected to be from currently unidentified microbial species. At the functional level, 8,312 KEGG orthologues and 59,303 eggNOG orthologue groups were identified in the updated gene catalogue. Additionally, 24% of the genes which were not mapped to the orthologue groups could be clustered into 649 novel gene families, which did not have any assigned function but were still included in the analysis as novel eggNOG groups.

Diversity and rarefaction analysis.Estimation of total gene richness, α-diversity (within-sample diversity), and β-diversity (between-sample diversity) in the set of 60 metagenomic samples was performed by randomized sampling and replacement, and estimates were compared to a different group of samples. Rarefied matrices were obtained by rarefying at 6 million reads per sample. In total, we performed 10 repetitions, and in each of these, we measured richness, α-diversity (by using the Shannon index), and β-diversity (by using Bray-Curtis distance) for each sample. The median values were taken as the respective measurement for each sample. Intersample distances were calculated using the Bray-Curtis distance matrices. The significance of the association with health status was performed using the Wilcoxon rank sum test.

Phylogenetic assignment of reads.A total of 4,097 reference microbial genomes were obtained from the Human Microbiome Project (HMP) and National Center for Biotechnology Information (NCBI) on 5 December 2015. The databases were independently indexed into two Bowtie indexes using Bowtie 2 (v2.3.4.1) (57). The metagenomic reads were aligned to the reference microbial genomes using Bowtie 2. The mapped reads from the two indexes were merged by selecting the alignment having the higher identity (≥90% identity). The percent identity was calculated using the formula: % identity = 100 × (matches/total aligned length). The normalized abundance of a microbial genome was calculated by summing the total number of reads aligned to its reference genome, normalized by the genome length and the total number of reads in the data set. For reads showing hits to the two indexed databases with equal identity, each genome was assigned an 0.5 read count. The relative abundance of each genome was calculated by adding the normalized abundance of each genome divided by the total abundance. The Calinski-Harabasz index (CHI) was used to calculate the variance between the clusters compared to the variance within clusters (58). A clade-specific-marker-based taxonomic assignment of reads was also done using the mOTUs (v2) (59) approach and Metaphlan (v2.0) (60).

Construction of metagenomic species for MGWAS.The gene cohort and its abundance from 291 samples belonging to India (60), Austria (103), and China (128) were combined and used for determining MGS/CAGs. The Pearson's correlation coefficient (PCC) cutoff of ≥0.9 was used for considering association between genes, and only genes having an abundance of >0 in at least 30 samples were considered for association analysis. Furthermore, the genes for which ≥90% abundance was obtained from a single sample were discarded. To determine the taxonomic origin of each MGS/CAG (metagenomic cluster), all the genes were aligned against reference microbial genomes of 4,097 genomes from HMP and NCBI at nucleotide level using BLASTN. The alignment hits were filtered using an E value of ≤10⁻⁶ and alignment coverage of ≥80% of the gene length, and 2,687,688 genes showed alignments against the reference genomes. The remaining genes were aligned against the UNIREF database (UniRef50) at protein sequences (61). The multiple best hits with equal identity and scores were further assigned taxonomy based on the lowest common ancestor (LCA) method. The genes were finally assigned to taxa based on comprehensive parameters of sequence similarity across phylogenetic ranks as described earlier (62). The identity threshold of ≥95% was used for assignment up to species level, ≥85% identity threshold was used for assignment up to genus level, and ≥65% identity was used for phylum-level assignment using BLASTN. The taxonomic assignments of MGS/CAGs were performed with the criterion that ≥50% of genes in each MGS should map to the same lowest phylogenetic group. So, if a particular species is assigned ≥50% genes out of the total, the assignment will be carried out at species level rather than at the level of genus or higher orders. The relative abundance of MGS/CAGs in each sample was estimated by using relative abundance values of all genes from that MGS/CAG. A Poisson distribution was fitted to the relative abundance values of the data. The mean estimated from Poisson distribution was assigned as the relative abundance of that MGS. The profiles of MGS/CAGs were generated and used for further analysis. The MGS/CAGs associated with CRC in the Indian population were scored using log odds ratio, and P values were calculated using the Wilcoxon rank sum test between CRC and healthy individuals. The Wilcoxon rank sum test was adjusted for multiple comparisons using false-discovery rate (FDR) adjustment. The MGS having P values of <0.05 and log odds ratio of >2 (CRC) or ≤2 (healthy) were considered enriched in CRC or healthy groups, respectively.

Fecal metabolomic sample preparation and derivatization.In order to identify the metabolic potential of microbes, metabolomics profiling of a subset of individuals (n = 18; CRC patients = 9, healthy = 9) was performed. Lyophilized fecal samples were used to achieve better metabolite coverage, as described previously (12). Metabolites were extracted from 80 mg of lyophilized samples in 1 ml of ice-cold methanol-water (8:2) by bead beating for 30 cycles (each cycle included 30 s of beating at 2,500 rpm and 1 min of standing at 4°C). The samples were then sonicated for 30 min in a probe-based sonicator (Branson digital Sonifier, model 102 C with double-step microtip) at 4°C followed by 2 min of vortexing. The supernatant was extracted by centrifugation at 18,000 × g for 15 min at 4°C and dried at 50°C under a gentle stream of nitrogen gas. To remove the residual water molecules from the samples, 100 μl of toluene was added to the dry residue and evaporated completely at 50°C under nitrogen gas. The extracted metabolites were first derivatized with 50 μl of methoxyamine hydrochloride (MOX) in pyridine (20 mg/ml) at 60°C for 2 h, and the second derivatization was performed with 100 μl of N-methyl-N-(trimethylsilyl)trifluoroacetamide (MSTFA) in 1% trimethylchlorosilane (TMCS) at 60°C for 45 min to form trimethylsilyl (TMS) derivatives. Finally, 150 μl of the TMS derivatives was transferred into a GC glass vial insert and subjected to GC-time of flight MS (TOFMS) analysis.

GC-MS analysis.GC-MS was performed on an Agilent 7890A gas chromatograph with a 5975 C MS system. An HP-5 (25-m by 320-μm by 0.25-μm-inside-diameter [i.d.]) fused silica capillary column (Agilent J&W Scientific, Folsom, CA) was used with the open split interface. The injector, transfer line, and ion source temperatures were maintained at 220, 220, and 250°C, respectively. Oven temperature was programmed at 70°C for 0.2 min and increased at 10°C/min to 270°C, where it was sustained for 5 min, and further increased at 40°C/min to 310°C, where it was held for 11 min. The MS was operated in the electron impact ionization mode at 70 eV. Mass data were acquired in full scan mode from m/z 40 to 600 with an acquisition rate of 20 spectra per second. To detect retention time shifts and enable Kovats retention index (RI) calculation, a standard alkane series mixture (C₁₀ to C₄₀) was injected periodically during the sample analysis. RIs are relative retention times normalized to n-alkanes eluted adjacently. The injector port temperature was held at 250°C, and the helium gas flow rate was set to 1 ml/min at an initial oven temperature of 50°C. The oven temperature was increased at 10°C/min to 310°C for 11 min, and mass data were acquired in full scan mode from m/z 40 to 600 with an acquisition rate of 20 spectra per second.

Metabolomic analysis.The raw peaks were processed for peak identification and alignment using the XCMS package in R. Initially, prior to alignment the parameters used for peak picking and retention time correction were optimized using the IPO package in R. The package using iterative modes and a range of values optimizes the best parameter settings for the GC-MS experiment. Using the centwave algorithm, the peaks were detected from 18 samples and were further corrected for retention time corrections using a “binwidth” of 1. The parameters optimized for preprocessing of the peaks from the GC-MS experiment are as follows: minimum peak width = 3, maximum peak width = 18.75, parts per million [ppm] = 575, difference in m/z values = 0.004588, signal-to-noise ratio [S/N] threshold = 10, bin width size = 1. The peaks were annotated using NIST/Massbank databases.

The most prominent features of the GC-MS data which could be annotated were used for differential analysis. In total, 234 features that were well annotated and showed prominent peaks were used. The differential analysis of peaks was performed using the MetaboDiff package in R. The concentrations of fecal metabolites in CRC and healthy subjects were quantified from peak intensities and were further normalized.

Meta-analysis.In order to define a non-reference-based approach of assigning genetic markers to a particular condition (CRC in this case), we used the MGWAS approach, which uses correlations among genes to generate clusters of genes. Various methods have been described for clustering genes in previous studies (51, 63). We here used a recently published method of determining clusters of genes using MGS-canopy-based correlations (55). The approach uses a matrix of gene abundance across samples and the correlation coefficient cutoffs used to identify fine-grained clusters of genes belonging to single species or closely related species.

The gene cohort and its abundance from 291 samples belonging to India (60), Austria (103), and China (129) were combined and used for determining MGS/CAGs. The PCC cutoff of ≥0.9 was used for considering association between genes, and only genes having an abundance of >0 in at least 30 samples were considered for association analysis. Furthermore, the genes for which ≥90% abundance was obtained from a single sample were discarded. To determine the taxonomic origin of each MGS/CAG (metagenomic cluster), all the genes were aligned against reference microbial genomes of 4,097 genomes from HMP and NCBI at nucleotide level using BLASTN. The alignment hits were filtered using an E value of ≤10⁻⁶ and alignment coverage of ≥80% of the gene length, and 2,687,688 genes showed alignments against the reference genomes. The remaining genes were aligned against the UNIREF database (UniRef50) at protein sequences (61). The multiple best hits with equal identity and scores were further assigned taxonomy based on the lowest common ancestor (LCA) method. The genes were finally assigned to taxa based on comprehensive parameters of sequence similarity across phylogenetic ranks as described earlier (62). The identity threshold of ≥95% was used for assignment up to species level, an ≥85% identity threshold was used for assignment up to genus level, and ≥65% identity was used for phylum-level assignment using BLASTN. The taxonomic assignments of MGS/CAGs were performed with the criteria that ≥50% of genes in each MGS should map to the same lowest phylogenetic group. So, if a particular species is assigned ≥50% genes out of the total, the assignment will be carried out at species level rather than at the level of genus or higher orders. The relative abundance of MGS/CAGs in each sample was estimated by using relative abundance values of all genes from that MGS/CAG. A Poisson distribution was fitted to the relative abundance values of the data. The mean estimated from Poisson distribution was assigned as the relative abundance of that MGS. The profiles of MGS/CAGs were generated and used for further analysis.

Global comparative analysis.The MGS/CAG abundance table was rarefied for 100 times at a lowest sequencing depth of 220,000 sequences/sample. The average of these rarefied counts from 100 iterations was calculated and used for further analysis. The principal-component analysis was performed on the MGS abundance table, and the components were correlated with covariates (country and status) using polyserial correlation. PCA with polyserial correlations showed countrywide variations as a major factor correlated with PC1.

Since variations between countries were higher than those between CRC patients, we performed multivariate analysis distance-based redundancy analysis (db-RDA). Canberra distances were selected based on the highest rank index for performing db-RDA. The metainformation used was the status (CRC/healthy) and country/study from which the data set came (India, China, and Austria). The Bray-Curtis distances were used for constrained ordinations and were constrained for country and status. The function “capscale” was used for performing ordinations with country and status separately, and the ordinations were plotted on a PCA plot.

Identifying global taxonomic patterns.In order to identify global CRC MGS/CAG, we performed univariate testing for differential abundance of MGS/CAG. Since countries and populations are stratified and have variations between them, we controlled for population variations by applying block. The independence test was performed using ytrafo = rank for Wilcoxon rank sum test and teststat = “scalar” using the COIN package in R. The population variations were controlled, and differential abundances of MGS/CAGs were calculated.

Gene biomarker identification.For gene biomarker identification, the 60 samples were divided into two cohorts: cohort A and cohort B. Cohort A comprised 48 samples and cohort B comprised 12 samples, by randomly selecting from each location and health status. Using the samples of cohort A, a group of genes that were highly correlated with each other (Pearson correlation ρ > 0.9) were identified, and the longest gene from each correlated group was used to construct a statistically nonredundant set. Further, we used the “CfsSubsetEval” method from Weka to identify a subset of genes that are highly correlated with the health status while having low intercorrelation with each other. The genes from this subset were further validated using the Boruta algorithm, which uses Random Forest to perform a top-down search for relevant features by comparing original attribute importance with importance achievable at random and eliminates irrelevant features to stabilize the test. To test the accuracy of these markers, a Random Forest model was constructed using these genes and was used for making the prediction on samples from cohort B.

CRC index.To compare the performances of markers, we computed a CRC index, as defined by Yu et al. (6), for each of the individuals on the basis of 33 gene markers identified using the methodology mentioned above.

Supervised learning.Predictive models were built using supervised machine learning algorithm Random Forest (RF). The models were optimized using 10,000 trees and default settings of mtry (number for variables used to build the model). The mean 3-fold cross-validation error rates were calculated for each of the binary trees and the ensemble of trees. The mean decrease in accuracy, which is the increase in error rates on leaving the variable out, was calculated for each prediction and tree and was used to estimate the importance score. The variables showing a higher mean decrease in accuracy of prediction were considered important for the segregation of the data sets into groups based on the categorical variable.

Network plot.In order to derive associations between microbial markers, a species cooccurrence network was generated from pairwise correlations using sparCC, which takes into account the compositional data and estimates the correlations between species. The species associations with correlation coefficients of >0.3 were considered for construction of networks and inferring associations among species.

Statistical analysis.All the statistical comparisons between groups were performed using a nonparametric Wilcoxon rank sum test with FDR-adjusted P values to control for multiple comparisons. The correlations between two variables and the correlations within the variable were calculated using Spearman’s correlation coefficient with adjusted P values. The correlations between categorical and numeric variables were performed using polyserial correlation/biserial correlations. To identify the enrichment of enzymes/species associated with a host, odds ratio was used as a measure of the enrichment of an enzyme in a host. The odds ratio was calculated as OR (k) = [∑_S = LOC1 A_Sk/∑_S = LOC1 (∑_{i ≠ k} A_Si)]/[∑_S = LOC2 A_Sk/∑_S = LOC2 (∑_{i ≠ k} A_Si)], where A_Sk denotes abundance of enzyme k in sample S. Apart from that, the Reporter features algorithm was used for gene-set analysis of significant pathways associated with different groups of samples. The algorithm takes the adjusted P values and fold changes (log odds ratio) as input for each KO. The gene statistic is calculated based on the significant association of KO and its direction of change through which the pathway is scored by calculating the global P value. All the graphs and plots were generated using the ggplot2 package in R.

Ethics approval and consent to participate.The recruitment of volunteers, sample collection, and other study-related procedures were carried out by following the guidelines and protocols approved by the Institute Ethics Committee of the Indian Institute of Science Education and Research (IISER), Bhopal, India. A written informed consent was obtained from all the subjects prior to any study-related procedures.

Availability of data and material.The data sets generated and/or analyzed during the current study are available in the NCBI BioProject database under project numbers PRJNA531273 and PRJNA397112.