Microbial Diversity and Composition is Associated with Patient-Reported Toxicity during Chemoradiation Therapy for Cervical Cancer

Study Design:

A cohort of 35 patients were included in this study with newly diagnosed locally advanced cervical cancer (clinical Stage IB1, IB2, IIA, IIB, IIIB, and IVA) with visible, exophytic tumor on speculum exam) in a multi institutional IRB approved prospective clinical trial between the University of Texas M D Anderson Cancer Center and Harris Health System Lyndon B. Johnson clinic. Patients with any previous pelvic radiation or treatment for cervical cancer were excluded. Enrolled patients followed planned treatment of intact cervical cancer with definitive radiation therapy, including external beam and brachytherapy with cisplatin.

Epic Bowel Function Scores

The Expanded Prostate Cancer Index Composite (EPIC) questionnaire is a patient reported toxicity instrument to evaluate bowel and urinary toxicity before, during and after pelvic irradiation. It has been previously validated in women undergoing pelvic irradiation for gynecologic cancers. The bowel domain of the EPIC questionnaire (EPIC-GI) specifically evaluates bowel function with both “function” and “bother” domains [13]. Patients were asked to complete the Bowel Domain of the EPIC survey, at baseline, week 1, 3, and 5 of treatment (T1, T2, T3, and T4) to evaluate frequency, urgency, and discomfort associated with bowel function with higher scores representing better GI function. Responses to individual questions were compiled into a composite score and bother subscale as have been reported previously. A toxicity score was calculated for each patient by subtracting the patient’s EPIC score at T4 (week5) from their score reported at T1 (Baseline). The median of those toxicity scores was calculated and patients with scores >16 were considered to be reporting high toxicity and those with scores ≤16 as reporting low toxicity.

16S Ribosomal DNA Sequencing

Stool samples were obtained from 35 patients before and during RT (week 1, 3 and 5) at the time of a rectal exam performed in clinic with routine pelvic examination. DNA was isolated from patient swabs and analyzed with 16S rRNA sequencing performed in collaboration with the Alkek Center for Metagenomics and Microbiome Research at Baylor College of Medicine. Samples were obtained using an isohelix swab (SK-2S) and placed into 20uL of Protease K and 400uL of Lysis buffer, and stored at −80C within one hour of sample collection.

16S rRNA sequencing was performed through the Alkek Center for Metagenomics and Microbiome Research at Baylor College of Medicine. 16S rRNA gene sequencing methods were adapted from the methods developed for the NIH-Human Microbiome Project [14]. Briefly, bacterial genomic DNA was extracted using the MagAttract Power Soil DNA Kit (Qiagen). The 16S rDNA V4 region was amplified by PCR and sequenced on the MiSeq platform (Illumina) using the 2×250 bp paired-end protocol yielding pair-end reads that overlap almost completely. The primers used for amplification contain adapters for MiSeq sequencing and single-end barcodes allowing pooling and direct sequencing of PCR products. The 16S rDNA gene pipeline data incorporates phylogenetic and alignment based approaches to maximize data resolution. The read pairs were de-multiplexed based on the unique molecular barcodes, and reads were merged using USEARCH v7.0.1090.

The 16S rRNA gene sequences were clustered into Operational Taxonomic Units (OTUs) at a similarity cutoff value of 97% using the UPARSE algorithm. OTUs were Mapped to an optimized version of the SILVA Database containing only the 16S v4 region to determine OTUs. A custom script constructed an OTU table from the output files generated in the previous two steps for downstream analyses of alpha-diversity, beta-diversity, and phylogenetic trends. ATIMA (R package consisting of APE and VEGAN) was used and samples were rarefied to a depth of 6000 (suppl.1). Inverted Simpson diversity (ISD) and Shannon Diversity index (SDI) was used to characterize alpha (within sample) diversity with comparison of changes over time and toxicity scores performed using linear regression. Principle coordinates analysis (PCOa) was used to compare beta (between sample) diversity between samples using the PERMANOVA test.

Association Network

The association networks (Anets) of samples were generated using OTU table with baseline samples (time point 0) and week 5 samples as described before [15, 16].Briefly, the OTU table was transformed to produce a new table with rows and columns comprised of samples and cells representing the shared richness (number of shared OTUs) for each pair of samples. Two samples were associated if their profiles in the produced table had Pearson correlation coefficient more than 0.4. The same threshold was used for both OTU tables to generate networks visualized in the Cytoscape v3.7.1. [17] and then clustered using the community structure analysis of biological networks [18]. The Fisher’s exact test was used to evaluate how significant the separation of high/low toxicity samples into clusters and to find individual OTUs overrepresented in high- versus low-toxicity patients at T1 and at week 5. Relative abundances of overrepresented OTUs at the time points were visualized as heat maps using Excel.

Linear discriminant analysis Effect Size (LefSe) Analysis

LefSe analysis was performed using the Galaxy web platform [19] to identify differentially represented bacterial taxa between groups with the Kruskal–Wallis test (alpha = 0.1, LDA threshold = 2.0, P values ≤0.1 were considered statistically significant). To highlight the differences in microbial composition at each time point, a toxicity score was calculated for each patient by subtracting the patient’s score at T4 from their score reported at T1. The resulting toxicity score was compiled from each patient, and the median of those toxicity scores was calculated (median decline = 16). Then each toxicity score was compared to the median score on an individual basis and thus categorized each patient as having either “HIGH” or “LOW” levels of toxicity; where “HIGH” toxicity patients had a toxicity score greater than the median decline in EPIC score between baseline and week 5. Each patient had EPIC scores for time points T1 and T4, but some patients lacked the EPIC scores for time points T2 and/or T3 (T1 n = 35, T2 n = 26, T3 = 26, T4 n = 35) (supplementary table1). Simple matrixes were generated using gut microbial metadata in accordance with required formatting parameters for LefSe analysis and then submitted to the Galaxy web platform [19].

Patient Characteristics

Thirty five patients with stage IBI-IVA disease undergoing definitive chemo radiation were enrolled. Patients with negative nodes or involvement of iliac nodes received pelvic radiation.

Patients with involved common iliac or paraortic nodes received extended field radiation. All treatment was given per institutional standard-of-care, including administration of antibiotics, chemotherapy and other medications. All patients were given the same dietary counseling considered standard-of-care for diarrhea management in our department. All patients received concurrent chemotherapy with weekly cisplatin at 40 mg/m2 (Table 1).

Patient reported GI toxicity

For all patients, EPIC bowel toxicity score declined over the course of treatment, reflecting an increase in symptom burden (figure 1A). There was a significant decrease in total EPIC score from baseline compared to week three and week five of treatment Baseline (85.6 ± 14.6) vs week 3 (71.8± 21.5); p = .008; Baseline (85.6 ± 14.6) vs week 5 (70.1 ± 18.8); p = .001, Figure 1A).

For the subset of scores that assessed the degree to which these symptoms bothered patients, there was also a significant decrease in the scores (88.9 ± 12.8 vs 72.9 ± 24.7 at week 1 vs week 3, p=0.008, 88.9 ± 12.8 vs 70.4 ± 21.6 at week 1 vs week 5, p = 0.0007, Figure 1B). The median decline in the total EPIC bowel score was 16 points. Patients with greater than or equal to the median decline were considered to have high toxicity (54.3%, N=19) while patients with less than the median decline wereconsidered to have low toxicity (45.7%; nZ16). Declinefrom baseline to week 5 was used as the primary endpointof the NRG 1203 study and so was chosen prospectively inthis trial as the most well-established metric to evaluatedecline in bowel function during radiation therapy.

Changes in Gut Microbiome over Time

Alpha diversity of the gut microbiome diversity decreased over the course of radiation (2.9± 0.5 at baseline to 2.49 ± 0.7; p=0.012 SDI, Figure 2A). Shifts in beta diversity (composition and structure) from baseline through the end of radiation treatment were also observed on principal coordinate analysis (PCOa) plots based on Jaccard distances (Figure 2B, R2 = 0.03, p=0.04). Analysis of change in relative abundance of taxa demonstrated a progressive decline in relative abundance of Clostridiales over time but relative stability of other phyla (Figure 2C).

Association of Gut Microbiome Alpha Diversity and GI toxicity

Greater microbiome alpha diversity was associated with higher EPIC bowel scores at all individual time points (R2 = 0.062, p= 0.017 for Shannon, Figure 3A); however alpha diversity at baseline was not significantly different for patients with subsequent development of high or low toxicity (Figure 3B).

Association of Gut Microbiome Beta Diversity and GI toxicity

Beta diversity refers to composition differences, or between sample differences in presence of specific species. Two methods were undertaken to explore beta diversity: association network (ANET) and linear effect size discriminant analysis (LefSe). Association networks were created based on similarity of profiles of shared species richness between samples. The analysis demonstrated similarity of profiles of shared species richness at baseline among all samples with high and low toxicity (Figure 3C) but distinctly different profiles at Week 5 ( Figure 3D) between patients with high and low toxicity (Fisher’s p-value <0.001). This reveals that radiation therapy imposes different changes on microbiome of patient experiencing high toxicity versus those with low toxicity and results in two divergent profiles associated with each group.

Considering the association of the shared richness with toxicity at the end of the therapy we focused further analysis on rare species and identified those rare individual OTUs that were over- or underrepresented (Fisher test p-value <0.05) in high- versus low-toxicity patients at week 5. We found nine OTUs that were over-represented in high-toxicity patients. Importantly, 5 out of the 9 OTUs annotated as Phascolarctobacterium, Lachnospiraceae, Veillonella, Erysipelotrichaceae, and Faecalitalea were found to be abundant in high-tox patients, when compared with baseline samples (Fig. 4) demonstrating outgrowth of these putative species after RT in high-toxicity patients.To further explore these differences on a higher level, LEFSe analysis was used to identify differences at the species level with differential abundance in patients who experienced higher toxicity over the course of treatment. At baseline (Figure 5A), again Clostridiales and Desulfovibrio had higher abundances in low toxicity patients, while, Sutterela, Finegoldia and Peptococcaceae (clostridia) had the highest abundance in high toxicity patients. This was consistent at week 5 (Figure 5B). At baseline, relative abundance of Clostridiales was significantly higher in patients with low toxicity. The abundance of individual families within the order Clostridiales was also higher in low toxicity patients (Figure 5C). Although certain subfamilies under class clostridia are associated with higher toxicity, as a whole clostridiales are found to be more abundant in low toxicity patients.