Impact of oral probiotic Lactobacillus acidophilus vaccine strains on the immune response and gut microbiome of mice

Summary

The potential role of probiotic bacteria as adjuvants in vaccine trials led to their use as nonparenteral live mucosal vaccine vectors. Yet, interactions between these vectors, the host and the microbiome are poorly understood. This study evaluates impact of three probiotic, Lactobacillus acidophilus, vector strains, and their interactions with the host’s immune response, on the gut microbiome. One strain expressed the membrane proximal external region from HIV-1 (MPER). The other two expressed MPER and either secreted interleukin-1ß (IL-1ß) or expressed the surface flagellin subunit C (FliC) as adjuvants. We also used MPER with rice bran as prebiotic supplement. We observed a strain dependent, differential effect suggesting that MPER and IL-1β induced a shift of the microbiome while FliC had minimal impact. Joint probiotic and prebiotic use resulted in a compound effect, highlighting a potential synbiotic approach to impact efficacy of vaccination. Careful consideration of constitutive adjuvants and use of prebiotics is needed depending on whether or not to target microbiome modulation to improve vaccine efficacy. No clear associations were observed between total or MPER-specific IgA and the microbiome suggesting a role for other immune mechanisms or a need to focus on IgA-bound, resident microbiota, most affected by an immune response.

Keywords:

Rice – Immune response – Vaccination and immunization – Vaccines – Microbiome – Immunologic adjuvants – Lactobacillus – probiotics

Introduction

The relationship between the microbiome and its host has been rigorously studied during the last decade resulting in evidence supporting its role in health and disease[1–5]. It was also established that the microbiome is dependent on diet and on its environment and could be modulated, for better or worse, by use of antibiotics, probiotics and/or prebiotics[6–11]. Use of probiotics has been shown to impact the gut microbiome due to direct competition or cooperation with the resident microbiota, and was shown to have a direct effect on the host through immune modulation and enhanced barrier function[12–15].

It is clear that the microbiome greatly influences mucosal health[3,16], but how vaccines, and any subsequent mucosal immune response, influence the microbiome is poorly understood[17]. There is increasing evidence that the immunogenicity and efficacy of current vaccines are related to the intestinal microbiome[18–20]. Multiple studies have also shown that probiotic administration prior to or concurrent with vaccination enhances B cell and antibody responses and provides the mucosa with direct protection from infection through interactions with the innate immune system[20–23]. Trials of both parenteral and nonparenteral vaccines, in conjunction with probiotic administration, also point to probiotic bacteria as adjuvants[17]. The mechanisms behind this phenomenon are incompletely understood but are likely due to probiotic surface structures[24]. This inherent adjuvanticity has led to the use of probiotics, typically lactic acid bacteria, as nonparenteral live mucosal vaccine vectors[25–27].

Species and subspecies of the genus Lactobacillus are an important and heavily studied group of Gram-positive lactic acid bacteria used for food preservation, food bioprocessing, and as probiotics. Most lactobacilli possess acid and bile salt tolerance, allowing them to survive the hostile environment of the stomach and proximal duodenum[28–30]. Additionally, several cell surface components of lactobacilli are recognized by immune cells via pattern recognition receptors (PRR)[31]. In particular, lipoteichoic acid (LTA), peptidoglycan (PG), and muramyl dipeptide (the subcomponent of PG) are the major immune stimulators recognized by the heterodimeric Toll-like receptor (TLR) 2/6 and nucleotide-binding oligomerization domain 2 (NOD2), respectively[32–34]. This capacity to interact with the innate immune system helps explain why some species of lactobacilli are effective inducers of mucosal antibodies, especially IgA[35]. The probiotic strain Lactobacillus acidophilus NCFM is particularly promising as an oral vaccine vector for several reasons: (1) immune stimulation via PRRs as was just described, as well as binding to dendritic cells (DCs) through DC-specific intercellular adhesion molecule 3 (ICAM-3)-grabbing nonintegrin (DC-SIGN)[36], (2) acid and bile tolerance[29,30], and (3) expression of mucus-binding proteins and association with the mucosal epithelium[37,38].

In this study we evaluated the impact of Lactobacillus acidophilus, as a probiotic vaccine vector, and its interactions with the host’s immune response, on the gut microbiome. We utilized three modified L. acidophilus strains with different constitutive adjuvants. All three strains expressed the membrane proximal external region (MPER) from Human Immunodeficiency Virus 1 (HIV-1) within the context of the major Surface-layer protein A (SlpA) that was developed in previous work[39]. The MPER epitope alone is a very weak B-cell immunogen, so to increase immunogenicity the two additional L. acidophilus vaccine strains were modified to either secrete soluble interleukin-1ß (IL-1ß, an inflammatory cytokine) or surface-expressed flagellin protein C (FliC, a TLR5 agonist). Both of these adjuvant strains were previously identified for increasing immunogenicity against MPER[39–41]. In addition, we used the MPER-expressing L. acidophilus (no IL-1β or FliC) along with rice bran as a prebiotic supplement. Rice bran has previously shown adjuvant properties to rotavirus vaccination in pigs and to enhance growth of probiotics[42]. To our knowledge no probiotic vaccine has been tested for gut microbiome alterations, and prior evidence with other oral vectors suggests that oral vaccines do not cause significant perturbations to the host microbiome, unlike what we have observed.

Results

Vaccination-induced differences in alpha diversity

Results of model fitting highlighted differences in Chao1 and Shannon diversity measures between the vaccine treatments over time (Fig 1 and S1 Table). Vaccine treatments used included Lactobacillus acidophilus (LA) expressing MPER as vaccine vector (MPER), LA expressing MPER and secreting interleukin-1ß as adjuvant (IL1b), LA expressing MPER and also expressing surface flagellin subunit C as adjuvant (FliC), LA expressing MPER and a custom rice bran diet as prebiotic (RB), LA wild-type strain as positive control (WT), and a negative control group (NG) (Materials and methods). The model describing Chao1 richness (S1 Table) highlights variation between treatments and between mice within the same treatment—it also includes a nonlinear temporal trend driven by the FliC and RB treatments (Fig 1A). FliC witnessed an increase in richness for weeks 2 through 6 that was mostly similar to that of NG and WT but significantly different from that of IL1b and RB for all three time points and from that of MPER for time points 2 and 4. FliC richness declined during weeks 8 and 10 and was similar to that observed for RB for week 10, which was qualitatively lower than all other treatments. Richness of the RB treatment was lower than all other treatments starting at week 2 and continuing to the end of the experiment. In addition to being significantly lower than that of FliC, RB richness was also significantly lower than that of NG and WT during weeks 2 through 6 highlighting the combined impact of a change in the microbial environment induced by both vaccine application and the presence of the rice bran prebiotic in the modified diet. The upward trend after the 4^th week could be interpreted as a possible recovery in richness as the mice adjusted to the change in diet and imposition of the vaccine. In general, IL1b and MPER had a similar trend and so did NG, FliC and WT, with NG having the lowest slope of increase in richness over time. IL1b had the second lowest Chao1 richness besides RB, also indicating a putative compound effect of an immune response due to the expression of IL-1ß and the impact of the vaccine vector. No significant differences were observed between the treatment levels at time 0 (highlighted by overlap between the 95% credibility intervals at that time point). Qualitatively, though, the negative control (NG) treatment showed somewhat higher richness than all other treatments followed by both FliC and the positive control (WT) at that point in time. All treatments seemed to converge to the same richness level by the end of the experiment.

Shannon diversity index accounts for both richness and abundance of the different OTUs within the observed samples. A simple linear regression model, with intercepts that varied per-treatment and per-mouse, was selected to describe the trends in this index (Fig 1B and S1 Table). Similar to Chao1, this model highlights differences due to the treatments and variation between mice within each treatment. Fig 1B highlights an increased diversity of the FliC treatment, compared to all treatments other than NG, that was significantly different from that under the MPER and the IL1b treatments, similar to what has been observed for Chao1. Shannon’s diversity for FliC seems to be similar to that of the negative control group. Both the MPER and the IL1b treatments were very similar and have the lowest diversity compared to all other treatment levels though not significantly different from the WT and RB treatments. Fig 1B indicates a linear increase in diversity as time progresses in the experiment. This increase might be a result of the mouse acclimation to the stress imposed by the gavage process used for vaccination.

S1 Fig compares the alpha diversity indexes between the different treatment levels for the cecal samples. The model identified to best fit these data did not show significant differences between the observed microbiomes, reflected in the overlap between the 95% credibility intervals observed in the figure. However, it is worth noting that, qualitatively, some of the trends observed in the fecal samples were also found in the cecal samples including the increased richness and diversity under the FliC treatment and the reduced richness under the RB treatment.

Vaccination-induced trends in Beta diversity

Separation between microbiomes associated with each treatment level is clear in the NMDS presented in Fig 2. An exception was MPER and IL1b, which seemed to overlap, corroborating the trends observed in the alpha diversity analyses above. Maximum separation was observed between the microbiome associated with the RB treatment indicating a possible compound effect of the prebiotic on the response of the microbiome to vaccination, also similar to above. Similar trends were observed for the cecal samples (S2 Fig) where MPER and IL1b, and the NG, seemed to be closer together than those of the other treatment levels.

Fig 3 provides a two-dimensional view of the trends of the beta diversity per-treatment. The figure indicates shifts in the microbial community structure occurring over time. This shift is clearly observed for the WT, MPER, and IL1b treatments where the microbial community seems to progress from its state at week 0 to a different state at week 10. This shift occurs in one step for the MPER treatment, where weeks 4–10 cluster at a new equilibrium, but gradually for the WT and the IL1b treatments. Although we observed separation between week 0 and week 10 for the FliC treatment, all time points clustered toward time zero. This observation is in concordance with the alpha diversity, indicating that vaccination supplemented with FliC has the least impact on the microbial community of the host. The RB treatment showed a profound jump in the microbial community structure from its state at time -1, one week prior to application of the treatment and the rice bran diet, and then to a new equilibrium state after progressing through the first treatment application at week 0. The negative control group did not show any significant trends as expected.

Taxonomic level trends and association with treatment groups

There were ten phyla observed within the microbial community structure within all mice under the different treatments for all time points. Five of these phyla were found to have relative abundance larger than 1% per-sample and these are: Firmicutes, Bacteroidetes, Tenericutes, Actinobacteria and Proteobacteria—presented in Fig 4 for the fecal samples and in S3 Fig for the cecal samples. Actinobacteria was only detected in the negative control while Proteobacteria was observed in the MPER treatment. As expected, Firmicutes and Bacteroidetes dominated the mouse gut microbiome in all samples conforming to previous observations[4].

Using Random Forests (RF) our goal was to identify drivers of the putative differences in response to the treatments. We added time as a feature to account for the possible temporal effect in discriminating between treatments. The optimal number of features (OTUs and time) used in constructing the decision trees in the random forest was identified to be 53 for the fecal samples and 9 for the cecal samples. This number was observed to minimize the out-of-bag (OOB) error rate (S4 Fig and Materials and methods).

For the fecal samples the overall OOB estimate of error rate was found to be about 20%. That is, 80% of the time samples were assigned to their correct treatment. The OOB error rate ranged from 14.3% to 24.1% per treatment, as shown in the confusion matrix presented in S2 Table, with assignment to FliC having the lowest error rate and assignment to WT having the highest. The table indicates closeness between MPER and IL1b; 6 out of 35 samples belonging to IL1b were misclassified as MPER and 6 out of 32 MPER samples were misclassified as IL1b corroborating some of the results observed in association with the alpha and beta diversity above. The table also highlights non-reciprocal association between WT and IL1b where 5 out of 29 WT samples were misclassified as IL1b with only one IL1b sample misclassified as WT. S3 Table indicates that much of the misclassification occurred at time zero, the time when the microbial richness was observed to be similar as described above.

Fig 5 shows the relative Mean Decreasing Gini importance plot of the 30 most impactful features, all of which are OTUs. S4 Table associates these OTUs with their taxonomic assignment; the higher the Gini value the more impactful is the feature. It is interesting that time was not identified as highly important for classifying the observed samples. It is also worth noting that an OTU part of the genus Bifidobacterium (phylum Actinobacter) was the most impactful discriminator. Other observed OTUs belonged to the families Lachnospiraceae (16), Ruminococcaceae (5), Clostridales_vadinBB60_group (6), Streptococcaceae (1) and Bacteroidales_S24-7_group (1) (S4 Table). Ten of the Lachnospiraceae OTUs were not classified; the remaining six belonged to the genera Tyzzerella, Acetatifactor, Lachnospiraceae_UCG-001, Marvinbryantia, Lachnoclostridium and Lachnospiraceae_NK4A136_group. OTUs identified part of the Ruminococcaceae family belonged to the Ruminiclostridium, Anaerotruncus, Ruminococcaceae_UCG-014 and Oscillibacter genera. All of the taxa belonging to Clostridiales_vadinBB60_group belonged to the Clostridiales_vadinBB60_ge genus. All of these three families were part of the Firmicutes phylum. Streptococcaceae, phylum Firmicutes as well, was the only Bacilli; Lactococcus was the genus associated with the OTU observed within this family and was ranked sixth in its impact on the correct classification. The observed OTU within family Bacteroidales_S24-7_group belonged to the genus Bacteroidales_S24-7_group_ge. This was the only OTU part of the phylum Bacteroidetes that was identified to impact the correct classification.

Fig 6 shows the trends observed for the five most impactful taxa identified within the fecal samples. Bifidobacterium (OTU0088) was mostly observed in the NG treatment and at week 10 of the WT treatment. Lachnospiraceae_unclassified (OTU0059) was observed for the most part in association with FliC and WT. Ruminiclostridium (OTU0149), on the other hand, was mostly present in NG, MPER, IL1b and WT but not FliC or RB. Clostridiales_vadinBB60_group_ge (OTU0180) was present mostly in MPER and FliC and to a lesser degree in IL1b and NG but not in RB and WT. Tyzzerella (OTU0097) was present for the most part in RB times 0–12 but to a much lower degree in RB time 0 and in the other treatments. Figures describing the other 25 taxa can be found in S5 Fig.

The overall OOB estimated error rate for the cecal samples was about 20.59%. That meant that samples were correctly classified to the respective treatment about 79% of the time, similar to the fecal samples. The per-treatment OOB error rate ranged from 0% for RB to 50% MPER (S5 Table). This result highlights separation of the RB treatment from all others indicating possible interaction or direct impact of the rice bran components on the microbial community response to vaccination. MPER samples did not have sufficient signal to discriminate them from other samples. The OOB for all other treatments was 16.7% (1 in 6 erroneous misclassifications). Note that the small sample sizes were not conducive for accurate classification for this experiment. S6 Fig shows the relative Mean Decreasing Gini importance plot of the 30 most impactful features for the cecal samples. S6 Table associates these OTUs with their taxonomic assignments. Unlike the observed fecal samples, Bifidobacterium was not a driver of discrimination for the cecal samples. Family Lachnospiraceae was represented by 22 OTUs, 16 of which were either unclassified or uncultured; the remaining six belonged to Lachnospiraceae_NK4A136_group, Lachnospiraceae_ge, Lachnospiraceae_FCS020_group, Roseburia, Acetatifactor and Lachnoclostridium. Six taxa represented the family Ruminococcaceae and these belonged to genera Oscillibacter (2), Ruminiclostridium_6, Anaerotruncus, Ruminiclostridium_9 and Ruminiclostridium. Only one taxon belonged to the genus Lactobacillus and another to Clostridiales_vadinBB60_group_ge. The phylum Firmicutes was the only phylum observed within the top 30 impactful taxa. However, due to the near-uniformity of the Gini measure we can conclude that many of the observed taxa, and not only the top 30, had an impact on the correct classification. Many of the top 30 taxa in the cecal samples matched those highlighted in the fecal samples above.

S7 Fig shows the trends of the five most impactful taxa in the cecal samples. The top three of these OTUs (OTU0076, OTU0094 and OTU0039) were unclassified or uncultured Lachnospiraceae. OTU0076 was highly abundant in the WT samples compared to all others. OTU0094 and OTU0039 were highly abundant in IL1b, MPER and the NG samples. OTU0094 also had detectable presence in the RB community albeit at a low relative abundance. The fourth impactful OTU (OTU0012) belonged to the Lachnospiraceae_NK4A136_group genus. This genus was highly represented in the RB samples compared to all other samples. The fifth OTU (OTU0167) belonged to the Oscillibacter genus. This OTU was present in all samples with elevated relative abundance in both the IL1b and the WT treatments. Trends for the remaining 25 taxa can be found in S8 Fig.

Vaccination-induced trends in total and MPER-specific IgA

The model best describing the trend in total IgA production was one that included a treatment- and mouse-specific cubic trend and intercept, also highlighting the variability between and within treatments (Fig 7). The cubic trend was mostly driven by the FliC treatment with fluctuation of the total fecal IgA over time from high to low to high again. Worth noting is the behavior of the WT treatment where the total IgA peaked during the 6^th week and was significantly different than that observed during week 0 within that treatment. Total IgA production of the WT was lower than that of all other treatments at week 0. No significant correlations were observed between the total IgA and the alpha diversity measures described above (S9 Fig) indicating that effectors other than total IgA might be responsible for the changes observed in the diversity of the microbiome.

S10 Fig shows an increase in MPER-specific IgA over time in MPER, FliC and IL1b, all of which include the MPER epitope. This is in accordance with our expectations. The figure indicates that the most pronounced induction of this response was observed in treatment IL1b. Still, some of the mice in these experiments didn’t have a detectable level of this specific IgA. No significant correlations were observed between these trends and measures of microbial diversity (S9 Fig). A lack of MPER-specific IgA in WT and NG groups shows induction of this antibody was due to our vaccine and not non-specific antibody production.

Minimal associations between microbial drivers and IgA

There were no IgA measurements observed for the cecal samples or the RB treatment. Samples with missing IgA measurements were dropped. The presented analysis is for total IgA and antigen-specific IgA samples separately, to minimize the missing data removal. Bars in color represent significance at the 0.1 cutoff p-value, without correction for multiple testing, of positive (red) and negative (blue) Spearman’s correlation between normalized taxonomic abundance of the 30 most impactful taxa and total IgA (Fig 8 and S7 Table) and MPER-specific IgA (S11 Fig and S8 Table). There were no particular patterns observed in these correlations between OTUs and total IgA save possibly one with Clostridiales_vadinBB60_group_ge. Clostridiales_vadinBB60_group_ge seemed to be positively correlated with IgA in all types of vaccine applications, significantly at the FliC treatment level. On the other hand, it seemed to either be neutral (WT treatment) or negatively associated with total IgA (NG treatment). Mixed results were also observed for antigen-specific IgA with most of the high correlation observed between taxa in the MPER treatment. These results are not conclusive though might support our observation that factors other than IgA might be at play in driving the microbiome’s response to vaccination or that a more targeted approach, focusing more on the resident microbiota that are most affected by the immune response, is required to better assess these associations.

Discussion

In this paper we studied the impact of an oral vaccine vector on the murine gut microbiome, utilizing the probiotic L. acidophilus induced immune response. We evaluated this impact using an HIV-specific peptide vaccine (MPER) augmented with either FliC or IL-1β as adjuvants and compared its use against a regimen that included rice bran as a prebiotic, along with positive and negative controls.

Maintenance of alpha and beta diversity was observed in association with the FliC treatment as compared to the negative control. This treatment utilized recombinant L. acidophilus (rLA) allowing for co-expression of the MPER epitope and the Salmonella enterica serovar Typhimurium flagellin (FliC). Maintenance of diversity could result when there is no effect of rLA on the microbiome directly through competition and/or indirectly through immune response. Our previous work has shown that FliC is an effective adjuvant that improves the response to the vaccine by stimulating activity of the Toll-like receptor 5 (TLR5) present in dendritic cells[41]. Chassaing et al.[43] showed that a loss of TLR5 from DCs did not associate with inflammation or change the microbiome, but did result in a complete loss of the IL-22 immune response, which is induced by flagellin. This is in contrast to a loss of TLR5 on mouse intestinal epithelial cells, which resulted in low-grade inflammation along with metabolic syndrome and an inability to clear pathobionts[43]. Hence, They argued that this intestinal epithelial loss of TLR5 could impact the microbiome. Our results seem to align with the first scenario they describe highlighting a possible TLR5 immune response associated with DCs rather than epithelial cells, hence the observed minimal immune-response impact on the microbiome. Though one would have also expected possible competitive effect of the probiotic on the microbiome that was not observed. Further studies are required to shed more light on this result.

In support of the current literature, we found a putative effect of the probiotic, the vaccine, the prebiotic and their different permutations on the microbiome. This study is the first to describe these effects utilizing a probiotic-vectored vaccine. WT, MPER and IL1b treatments behaved similarly when comparing their beta diversity measures. These treatments resulted in a shift in the microbial community over time to a putative new equilibrium. Similarities between MPER and IL1b were also observed in the confusion matrix resulting from the random forest classification approach where the error of assignment from one of these groups to the other was comparable. While Shannon diversity highlighted a comparable reduction in diversity of these three treatments, as compared to the negative control, Chao1 highlighted a gradient decrease in richness with a clear reduction in the IL1b treatment (Fig 1). IL-1ß, part of the IL1 family, is a proinflammatory cytokine with a role in modulating the adaptive immune response [40,44]. Our results suggest that the observed richness decrease might be an outcome of the compound effect of the probiotic and the inflammatory effect of the adjuvant. An existing body of literature supports that probiotics and/or vaccination could result in modulation of the resident microbiome[13,18,45–50]. In addition, a significant decrease in richness and shift in the microbial community structure was also observed in association with the RB treatment. This treatment also stood out with zero error of assignment, utilizing the Random Forest approach, as observed for the cecal samples. These results could be attributed to the compound effect of the probiotic vaccine and the prebiotic rice bran in this case. Modulation of the microbiome utilizing rice bran[11,51] and the effects of diet are well supported[52–54]. Taken together these results highlight the opportunity for synbiotics, combining probiotics and prebiotics to improve efficacy of vaccination. Careful consideration of the constitutive adjuvants and the use of diet supplementation is needed when the target is modulation of the microbiome to improve efficacy of vaccination or maintenance of the resident microbiome.

The murine microbiome was mostly dominated by taxa belonging to the Firmicutes and Bacteroidetes phyla[4]. The discriminating taxa in both fecal and cecal samples were similar, though the small cecal sample sizes blurred the impact of these taxa. In both cases, taxa were mostly dominated by Firmicutes and belonged to the Lachnospiraceae family. This family belongs to the order Clostridiales with all known taxa being strictly anaerobic and found mainly in the mammalian gut[55]. The Actinobacteria phylum was less dominant in these samples and included an OTU that belonged to the Bifidobacterium genus. This OTU had the highest impact on discriminating between the treatment levels in the fecal samples. This was not a surprise though unlike in previous research, it was not its dominance in the rice bran treatment that caused it to have a prominent effect but its presence in the negative control. Bifidobacterium is also commonly found as part of the mammalian microbiomes[56]. In mammals and humans this genus dominates in early life neonates though still comprises 25% of the microbial population within the healthy gut of adult humans[57]. No taxa of the genus Lactobacillus were observed as part of the top 30 impactful taxa in the fecal samples and only one was observed in the cecal samples, albeit all taxa had similar weight in discriminating between the treatment levels based on these samples. A more targeted approach might be required to verify the possible impact of vaccination on these specific taxa (see below).

Total IgA also showed different trends between the applied treatments. Induction of MPER-specific IgA was observed in association with the treatments that included the vaccine, which was expected. Non-parametric correlation analysis, utilizing the Spearman correlation coefficient, between the diversity measures and between the abundance of the 30 most impactful taxa and total and antigen-specific IgA did not produce highly significant associations. This indicates that either other immunological factors, in addition to the direct impact of the vaccine vector, were at play or that a more targeted approach is required to better assess these associations. A possible more direct approach to assess the correlation between the change in IgA and its possible microbiome targets in the gut is to focus on the study of the IgA-coated microbiota. These taxa represent the putative resident microbiome most affected by changes in the immune response. Changes in abundance and presence/absence of components of this part of the microbiome could highlight possible off-target effects of the vaccine whereby resident taxa that are evolutionarily or functionally similar to the vaccine vector might be targeted by the immune response. We are currently assessing this hypothesis in a proof of principle follow-up study.

Materials and methods

Study design and vaccine strains

6–8 week old BALB/c mice (Jackson Laboratory) were housed together in isolated cages based on vaccine strain to be received. Animal gut microbiomes were normalized using multiple procedures including standardized diet, identical cage locations, cage swapping, and oral gavage of cecum content one week prior to first vaccination[58]. Six animals per group received either Lactobacillus acidophilus wild-type strain (WT), WT expressing the membrane proximal external region from HIV-1 within the context of the major S-layer protein A (MPER), MPER secreting interleukin-1ß (IL1b) or MPER expressing surface flagellin subunit C (FliC). Table 1 provides a full description of the used strains along with associated originating references. Five mice received MPER and a custom 10% rice bran diet (RB) from Envigo, prepared as previously described[59]. The six mice in the negative control group (NG) received only the dosing buffer composed of 8.4 mg/mL NaHCO₃ and 20 mg/mL soybean trypsin inhibitor (SIGMA T9128) in ultrapure water[29,60]. Bacterial cells were prepared from overnight culture and suspended in dosing buffer. Each animal received 5 x10⁹ CFU/day in 200 μL of dosing buffer by oral gavage for three consecutive days at weeks 0, 2, 4, 6, 8, and 10 (six 3-day doses) as previously described[39,61,62]. Samples were collected on the first day of each dose, prior to dose administration. Rice bran diet group samples were also collected one week prior to starting the RB diet (week -1) and at week 12.

The care and use of experimental animals was carried out in accordance with relevant guidelines and regulations, with approval from Colorado State University’s Institutional Animal Care and Use Committee (IACUC 14-5332A). All animals were age matched, maintained in specific pathogen free conditions, individually tracked and monitored daily for clinical signs of stress or illness, including but not limited to changes in skin and hair, eyes and mucous membranes, respiratory system, circulatory system, central nervous system, salivation, diarrhea, or lethargy. Upon arrival, animals were housed socially (n = 2–5) in commercially available, individually ventilated caging systems with a 12 h light/12 h dark cycle. Animals were provided ad libitum commercial irradiated rodent chow (Teklad Global) or the rice bran diet and tap water filtered via reverse osmosis in autoclaved water bottles; all bedding and enrichment materials were autoclaved prior to use and changed regularly.

Sample collection

Fecal pellets were collected from mice in cleaned collection cups and transferred to pre-weighed collection tubes. Pellets for supernatant extraction were weighed and 10 μL of 2x ProteaseArrest (G-Biosciences) was freshly diluted in cold PBS (from 50x stock) per mg of fecal pellet weight within 5 minutes of excretion. Samples were then homogenized via three 20-second cycles (4 m/s) using a FastPrep-24 Homogenizer (MP Biomedicals) and centrifuged at 10,000 RCF, 4°C for 10 minutes. The supernatant was removed and stored at -80°C for later analysis. Fecal pellets for microbiome analysis were collected directly from the anus into collection tubes, immediately placed on ice, and stored at -80°C until extraction. Cecal contents for microbiome analysis were obtained by necropsy at week 10 for treatments NG, WT, MPER, IL1b, and FliC, and at week 12 for RB, squeezed into collection tubes, immediately placed on ice, and stored at -80°C until extraction. Animals were euthanized by CO₂ asphyxiation and subsequent cervical dislocation, according to IACUC protocol.

ELISA

MaxiSorp high-binding 96-well plates (Nunc) were coated overnight at 4°C with 100 μL/well of 10 μg/mL goat anti-mouse IgA antibody (Bethyl Laboratories, Inc.) for total IgA ELISA, or with 1 μg/ml of synthetic 17-mer HIV-MPER peptide (GNEQELLELDKWASLWN, Bio-Synthesis Inc.) for antigen-specific ELISA, both suspended in carbonate coating buffer (15 mM Na₂CO₃, 35 mM NaHCO₃ in ultrapure water). Wells were blocked with 1% BSA in PBS for 1 hour at room temperature. After washing (0.05% Tween 20 in PBS), fecal supernatant was serially diluted 1:2 (from 1:10 to 1:640) in PBS with 1% BSA and 0.05% Tween 20, added to wells, and incubated for 2 hours at room temperature. Following washing, HRP-conjugated goat anti-mouse IgA antibody (50 ng/mL, Bethyl) was added and incubated for 1 hour at room temperature. Color development with 3,3’,5,5’-tetramethylbenzidine (TMB) was terminated after 15 minutes with sulfuric acid and absorbance (570–450 nm) was measured.

To quantify total IgA, 1:10 dilutions of Mouse Reference Serum (Bethyl Laboratories) starting at 1000 ng/mL were run on each plate in duplicate to generate a standard curve. Fecal supernatants from each animal’s week 0 time point were included in the assay to determine endpoint titer for antigen-specific IgA. The optical density cutoff was calculated as the mean value of all negative controls for each vaccination group + 3.365 standard deviations, based on the 99% confidence interval standard deviation cutoff multiplier for n = 6[63].

DNA extraction and sequencing

Microbial genomic DNA was extracted from whole fecal pellets using PowerFecal DNA Isolation Kit (QIAGEN); buffer-only controls were included in each extraction batch. The hyper variable region 4 (V4) of the 16S rDNA gene was amplified by PCR according to the Human Microbiome Project protocols[64,65]. Primers used for multiplexed 16S library generation included an Illumina adapter, an 8-nt index sequence for barcoding, a pad and linker, and the 16S-specific sequence, as described previously[65]. ZymoBIOMICS D6305 Microbial Community Standard (mock community) and no-template controls were included in each PCR plate. Dual-indexed library molecules were then purified using Sera-Mag beads (GE Life Sciences) to select for DNA fragments larger than 150bp. Purified library molecules were quantified using Quant-iT Broad-Range dsDNA kit (Invitrogen) according to manufacturer’s protocols. Libraries were pooled at equimolar concentrations and quality controlled prior to sequencing on an Illumina MiSeq at the Colorado State University’s Next Generation Sequencing Core Facility. Paired-end 2 x 250 reads and index reads for multiplexing were generated as previously described[65]. Samples were sequenced targeting an average read-depth of 4x10⁴ reads per sample. All raw sequence data were deposited and are available on the National Center for Biotechnology Information’s (NCBI) Sequence Read Archive (SRA) repository under accession number PRJNA542488.

Data processing and bioinformatics

Sequence data were processed using mothur[66] version 1.39.5 and utilizing the developers’ standard operating procedure (SOP) for OTU calling and taxonomic classification of MiSeq data first presented in Kozich, et al., 2013[65]. For alignment and classification within this SOP we used the SILVA database[67] version 128. Sequencing error rate was assessed using the ZymoBIOMICS^™ microbial community standard (mock community) described above, which includes 8 bacterial and 2 fungal species. Our sequencing included five mock community samples, nine negative controls with no sample, and four negative controls with no template. Clustering, for OTU identification, was performed using OptiClust[68] utilizing 0.97 sequence similarity. We used maximum counts per OTU of all mock samples to identify a cutoff number of reads per OTU to guard against overestimation of sample diversity. Rarefaction curves were generated using the package vegan[69] as implemented in R version 3.4.4[70] to assess diversity and suitability of depth of coverage per sample. The resulting OTU table was utilized in further data analyses as follows. Results of data processing and the associated bioinformatics can be found in the S1 Appendix.

Diversity and IgA analysis

Alpha diversity

Chao1 and observed richness along with Shannon and inverse Simpson diversity indices were computed utilizing the R package phyloseq[71]. We used a Bayesian framework to fit and compare 30 putative models that evaluate linear and nonlinear associations between time and these measures of diversity. Model fitting was done using the R-package rjags[72] an interface to the Just Another Gibbs Sampler (JAGS) software[73] and model selection was done using the Deviance Information Criterion (DIC)[74] also implemented in rjags. S2 Appendix provides detailed description of the general characteristics of the different models used. The model with minimum DIC was selected to be the best-fit model. Convergence of the Markov Chain Monte Carlo (MCMC) sampling was assessed using the trace diagrams associated with the chosen models and Gelman and Rubin’s convergence diagnostic[75] based on three chains that were run for 500,000 iterations using the first 200,000 as burn-in and sampling every 500 steps to assure sample independence. This resulted in 600 samples per chain (a total of 1,800 samples) that we used to construct the posterior distribution of the parameters. Times -1 and 12 associated with the RB treatment were dropped from this analysis to put all treatments on the same footing.

Beta diversity

Nonmetric Multidimensional Scaling (NMDS)[76] was used on the OTU level to assess possible trends and clustering in the microbial community structure per treatment condition and per time point. NMDS was performed separately per treatment level using the vegan package and utilizing Bray-Curtis dissimilarity[76] and was based on data transformed using Cumulative Sum Scaling[77].

IgA

Total IgA was compared between treatment levels, other than RB, using the same Bayesian framework described above for alpha diversity. Total IgA data were not collected for the RB treatment. MPER-specific IgA data were transformed by first replacing zero measurements with 1 and then computing the log. This was done to overcome the high variation observed in the measurements.

Trends and association on the taxonomic level

Bar graphs

Utilizing relative abundance data based on the resulting OTU table, bar graphs were generated using the ggplot2[78] package in R. These plots were generated for the relative abundance data at the phylum level and meant to describe the microbial community structure per time point under each of the treatments.

Random Forest Classification: We used Random Forests (RF)[79] to identify the major microbial taxa that influenced the accurate prediction of the state of vaccination (treatments) over the period of the experiment. We used an iterative approach to identify the optimal number of features (variables) to incorporate in constructing the regression trees in the RF. In this approach we iteratively fit a random forest including one all through 100 features using the mtryStart option in the function tuneRF in the randomForest[80] package. This was done 100 times and the out-of-bag (OOB) error was averaged over all runs. The number of features to include in sampling to construct the regression trees in the RF was chosen to be that with the minimum median OOB. We used the parameter ntree = 2000 for the total number of trees to grow in the forest and set the parameter importance to TRUE to assess importance of the different features in prediction. The predictive features included all OTUs in addition to time of vaccine administration per treatment level to account for the impact of possible change over time.

Association with immune response

We used Spearman’s correlation coefficient[81] to assess the association between time, total IgA, antigen-specific IgA and the alpha diversity measures. This was also done to study the associations between taxa identified as most predictive using RF and the immune response.

Supporting information

S1 Fig [pdf]
Differences in alpha diversity per treatment level in cecal samples.

S2 Fig [pdf]
Beta diversity per treatment level in cecal samples.

S3 Fig [pdf]
Bar-plots of phylum level diversity per treatment in cecal samples.

S4 Fig [pdf]
Out of Box (OOB) error plots.

S5 Fig [pdf]
Per time point and treatment level differential abundance of influential OTUs 6–30 identified by Random Forest analysis for fecal samples.

S6 Fig [pdf]
Mean decreased Gini importance for cecal samples.

S7 Fig [pdf]
The normalized abundance of the five most impactful OTUs in cecal samples.

S8 Fig [pdf]
Per time point and treatment level differential abundance of influential OTUs 6–30 identified by random forest analysis for cecal samples.

S9 Fig [pdf]
Per treatment correlation analysis between alpha diversity and total IgA, MPER-specific IgA and time.

S10 Fig [pdf]
Change in MPER-specific IgA per treatment over time.

S11 Fig [pdf]
Spearman correlation between MPER-specific IgA and the 30 influential OTUs per treatment.

S12 Fig [pdf]
Rarefaction curves.

S13 Fig [pdf]
Nonmetric multidimensional scaling plot identifying the outlier group at time period 2 of the wild type (WT) treatment.

S1 Table [docx]
Best fit models explaining trends in alpha diversity.

S2 Table [txt]
Confusion matrix of the random forest classification model per treatment for fecal samples.

S3 Table [csv]
True vs predicted assignment per sample using the Random Forests classifaction model for fecal samples.

S4 Table [csv]
Taxonomic assignment of the 30 influential OTUs sorted by their Gini inortance score for fecal samples.

S5 Table [txt]
Confusion matrix of the random forest classification model per treatment for cecal samples.

S6 Table [csv]
Taxonomic assignment of the 30 influential OTUs sorted by their Gini inortance score for cecal samples.

S7 Table [csv]
Spearman correlation tables between total-IgA and normalized abundance of the 30 influential OTUs for fecal samples.

S8 Table [csv]
Spearman correlation tables between MPER-specific IgA and normalized abundance of the 30 influential OTUs for fecal samples.

S1 Appendix [docx]
Results of data processing and bioinformatics.

S2 Appendix [pdf]
Univariate models for analysis of alpha diversity and total-IgA data.