Twice-Daily Application of HIV Microbicides Alters the Vaginal Microbiota

ABSTRACT Vaginal HIV microbicides offer great promise in preventing HIV transmission, but failures of phase 3 clinical trials, in which microbicide-treated subjects had an increased risk of HIV transmission, raised concerns about endpoints used to evaluate microbicide safety. A possible explanation for the increased transmission risk is that the agents shifted the vaginal bacterial community, resulting in loss of natural protection and enhanced HIV transmission susceptibility. We characterized vaginal microbiota, using pyrosequencing of bar-coded 16S rRNA gene fragments, in samples from 35 healthy, sexually abstinent female volunteer subjects (ages 18 to 50 years) with regular menses in a repeat phase 1 study of twice-daily application over 13.5 days of 1 of 3 gel products: a hydroxyethylcellulose (HEC)-based “universal” placebo (10 subjects), 6% cellulose sulfate (CS; 13 subjects), and 4% nonoxynol-9 (N-9; 12 subjects). We used mixed effects models inferred using Bayesian Markov chain Monte Carlo methods, which showed that treatment with active agents shifted the microbiota toward a community type lacking significant numbers of Lactobacillus spp. and dominated by strict anaerobes. This state of the vaginal microbiota was associated with a low or intermediate Nugent score and was not identical to bacterial vaginosis, an HIV transmission risk factor. The placebo arm contained a higher proportion of communities dominated by Lactobacillus spp., particularly L. crispatus, throughout treatment. The data suggest that molecular evaluation of microbicide effects on vaginal microbiota may be a critical endpoint that should be incorporated in early clinical assessment of microbicide candidates.

(3). Even more recently, a different clinical trial (VOICE) of the same tenofovir gel but with daily application was ended prematurely by the Data Safety Monitoring Board because of a lack of evidence of a beneficial effect (i.e., futility).
Among the first vaginal microbicide agents studied were nonoxynol-9 (N-9), a nonionic surfactant still widely used as an FDA-approved spermicide, and cellulose sulfate (CS), a highmolecular-weight sulfated carboxymethylcellulose polymer. While initial in vitro and early-phase clinical studies of these microbicides were promising, the results of larger phase 2B and 3 trials showed no protection against HIV compared to placebo treatment and an increased risk of infection when the microbicides were used very frequently (4)(5)(6). The reasons for the failure of these agents remain unknown. One hypothesis that may account for the failure of the microbicides holds that microbicide application alters the vaginal microbiota so as to yield a vaginal environment that has lost its natural protective abilities, either directly enhancing HIV transmission or indirectly acting to activate potential host cells, which would facilitate HIV transmission. This hypothesis was further supported by recent work using an in vitro vaginal microbiota colonization model system (7).

RESULTS AND DISCUSSION
Initial studies found that the microbicides CS and N-9 had limited effects on conventionally cultured organisms (8,9) (mainly Lactobacillus sp.), which was not surprising, since cultivationdependent methods provide biased quantitative and incomplete qualitative information on the composition of bacterial communities (10,11). In this study, we undertook a comprehensive analysis of the effects of CS and N-9 on the vaginal microbiota in a repeat phase 1 study of those agents, using culture-independent molecular methods and Bayesian statistical modeling. Vaginal swabs were collected according to the study design outlined in Fig. 1A. To comprehensively evaluate the effects of the microbicides on the vaginal microbiota, we characterized the vaginal microbial community taxon composition and relative abundance using pyrosequencing of bar-coded 16S rRNA gene fragments (12). A total of 146 longitudinal samples from 35 subjects were successfully collected and analyzed ( Table 1). The microbicides N-9 and CS, and to some extend the hydroxyethylcellulose (HEC) placebo, are major inhibitors of PCR amplification. Using the modified whole genomic DNA extraction procedure developed for this study, we generated DNA from which variable regions 1 and 2 (V1 and V2) of 16S rRNA genes were successfully amplified. Pyrosequencing of these bar-coded 16S rRNA gene amplicons produced a data set consisting of 791,295 high-quality sequence reads with an average length of 359 bp and 5,420 reads per sample. Overall, a total of 296 taxa were observed in the vaginal microbiota of these women. The depth of coverage for each community was sufficient to detect taxa that constituted~0.1% of the community. Complete linkage hierarchical clustering methods were applied and revealed five major bacterial community state types (CSTs) (Fig. 1B). Three CSTs were often dominated by different Lactobacillus species: L. crispatus (CST I), L. iners (CST III), or L. gasseri (CST II) ( Fig. 1B; see also Table S1 in the supplemental material). CST IV-A and CST IV-B were heterogeneous in composition, without significant numbers of Lactobacillus spp., but differed from each other in composition. The frequencies of each CST in this cohort were similar to those previously published (12,13), except that CST V was not observed in the subjects enrolled in this study (Table 1), as its frequency in the general population is less than 2% (12). Vaginal bacterial communities from samples that clustered in CST IV-B were characterized by higher proportions and types of anaerobic bacteria ( Fig. 1B; see also Table S1 in the supplemental material) such as Atopobium, Prevotella, Megaspherae, Sneathia, and Mobiluncus, as well as Gardnerella (the latter in proportions ranging from 0.2% to 9.9%), while CST IV-A comprises members of the genera Streptococcus, Enterococcus, and Escherichia, as well as small proportions of Lactobacillus spp. (Fig. 1B; see also Table S1 in the supplemental material). CST IV-B microbial composition is consistent with vaginal communities found in women with bacterial vaginosis (14,15) and was associated with higher Nugent scores, while CST IV-A was associated with low and intermediate Nugent scores, as shown by modeling this interaction using a log-linear model ( Fig. 2A). Further, using a mixed effect logistic regression model in which the presence or absence of G. vaginalis was the outcome variable and CST was the predictor variable, the log odds ratio of the presence of G. vaginalis in CST IV-B was shown to be significantly (P Ͻ 0.001) higher than in CST IV-A. Other CSTs did not show log odds ratios of the presence of G. vaginalis significantly different from those seen with CST IV-A. This holds true if the abundance of G. vaginalis (ranging from 0.2% to 9.9%) was modeled similarly (p Ͻ 1eϪ16), supporting the correlation observed between high Nugent score and CST IV-B (16). Because high Nugent scores and bacterial vaginosis have been associated with increased transmission and acquisition of HIV (17,18), we evaluated the possibility that CST IV-B was associated with the application of N-9 or CS. However, the analysis revealed that 30.8% of the samples were assigned to CST IV-A (Table 2) and that CST IV-A, not CST IV-B, showed a statistically significant association with N-9 and CS use at visit V4, during the middle of the product application period, compared to placebo and visits V2 and V3 (baseline) ( Fig. 2B and 3A). This result, while surprising, was supported by an analysis of the frequencies of Nugent score categories and product applications over the study period that showed no differences between the N-9, CS, or placebo arms (Fig. 3B). Application of the placebo, HEC, does not appear to affect the frequency of community state types, as shown by a consistent CST frequency distribution at each visit before (V3) and during (V4) the application period and immediately after placebo use ceased (V5) (Fig. 3A). Interestingly, a shift was observed between visits V5 and V6, when an increase in the proportion of CST IV-A was observed (Fig. 3A). However, that shift was not significantly associated with use of the placebo but certainly represented the normal temporal dynamics of the vaginal microbiota, often observed with the L. iners-dominated CST-III (13).
To better characterize the effect of N-9 and CS on shifting the vaginal microbial communities toward CST IV-A, we evaluated the ability of the communities to not shift to CST IV-A by modeling the distances of each community state from the center of CST IV-A (see Fig. S1 in the supplemental material). The model took into account that a subsequent serial microbial community sample in a particular subject depends on the previous microbial community (13) and was designed using mixed effects models inferred using Bayesian Markov chain Monte Carlo methods (see Materials and Methods). Figure 4 shows the interaction plot of the mean values and error bars of each sample's community state with respect to the distance from the center of CST IV-A for each treatment arm and visit, computed using the mixed effects model.
Community state distances to the center of CST IV-A for the N-9 and CS groups were significantly different from those for the placebo group at visit V4 and showed borderline significance at visit 5 ( Fig. 4). This analysis highlights the effect of twice-daily appli-cations of N-9 or CS microbicides on the composition of vaginal bacterial communities after 7 days of use. Upon application of N-9 or CS, most communities shifted in composition from community states that were dominated by species of Lactobacillus to com-  Table 2. munity states mainly dominated by anaerobes and by members of the genera Streptococcus, Enterococcus, and Escherichia (Fig. 1B).
Interestingly, it appears that the vaginal bacterial communities have the ability to rebound rapidly, as no statistically significant difference was observed between placebo and N-9 or CS use at visit 6, 3 days after the last use of the products (Fig. 4). Because CST IV-A contains low number of Lactobacillus spp., it is anticipated that, while not necessarily characterized by high Nugent scores, this state, like CST IV-B, would contribute to an increased risk of transmission or acquisition of sexually transmitted infections (STIs), including HIV. Our results are consistent with previous culture-based evaluations of N-9 (8,9), which found that daily application of N-9 promoted the loss of non-H 2 O 2producing lactobacilli and did not affect vaginal colonization by H 2 O 2 -producing lactobacilli (8); however, that earlier study did not evaluate the effect on other bacteria present in the vagina. The earlier study also found that N-9 application did not produce increases in the Nugent score. A previous phase I safety study of CS did not show increases in Nugent scores either; however, it did show a concomitant reduction in levels of H 2 O 2 -producing lactobacilli and an increase in Escherichia coli numbers (19). Furthermore, in that study CS significantly reduced Lactobacillus sp. colonization of human cervicovaginal epithelial cells and tissues and increased their proinflammatory reaction to bacteria (7). Other polyanions similar to CS have been shown to interfere with Tolllike receptor (TLR)-mediated responses in human cervicovaginal cells (20). Effects of the candidate microbicides on vaginal cell signaling pathways may contribute to observed changes in CSTs that accompanied application of the microbicide candidates in our study.
If the N-9 and CS microbicide agents are associated with both an increased risk of HIV transmission and an alteration in the vaginal microbiota, it may be helpful to consider what may link these two phenomena. One potential plausible explanation is that a distortion in the microbiota may, for example, through interactions with pattern recognition receptors, lead to an increase in levels of inflammatory cytokines, which would activate potential HIV host cells. Activation of TLRs in human vaginal cells by microbial antigens induces innate and proinflammatory responses (21,22). Since HIV replicates preferentially in activated cells and since the risk for HIV transmission at each episode of intercourse  Table S4 in the supplemental material) stratified by community state type and Nugent category, with error bars indicating 95% credible intervals derived from model 2 (described in Materials and Methods). (B) Interaction plot of the mean log(odds) values and the corresponding 95% credible interval of logistic regression model 1, evaluating the association between CSTs (community state types) at visit V4 compared to visits V2 and V3 for each product use arm.
In the N-9 (V4) and CS (cellulose sulfate) (V4 and V5) arms, no communities were assigned to CST II. In addition, in the N-9 arm (V4), no communities were assigned to CST III (Fig. 4). Ravel et al. is low, increasing potential host cell activation may reasonably be expected to affect overall transmission rates. The present study, despite the relatively small number of subjects, demonstrated using culture-independent methods that application of microbicides associated with higher rates of HIV transmission can affect the composition of the vaginal microbiota, although it is impossible to say whether these changes are specific to the active agent or to one or more other components of the microbicide vehicle gel, given that the universal placebo has a composition different from those of the vehicles of both N-9 and CS gels. These changes in the vaginal microbiota could plausibly contribute to higher rates of HIV transmission. While our observation involves candidate microbicides that have been shown to be ineffective at preventing HIV transmission in large phase 3 clinical trials, an analogous effect could produce a decrease in the effectiveness of microbicides that do help prevent HIV transmission. The data suggest that molecular, culture-independent evaluation of microbicide effects on vaginal microbiota, in addition to or combined with culture-based methods, may be an important component of the early clinical assessment of candidate microbicides. The data also suggest that the development of a maximally effective HIV microbicide delivery system may require, in addition to agents focused solely on the virus, components that help maintain or promote a healthy vaginal microbial community.

MATERIALS AND METHODS
Clinical sampling and study design. The samples were obtained in the course of a carefully monitored repeat phase 1 study of CS and N-9, conducted by CONRAD (protocol no. A04-097). Sixty healthy, sexually abstinent female (tested using Rapid Stain Identification Series [RSID]-Semen) volunteer subjects (ages 18 to 50 years) with regular menstrual cycles (27 to 35 days) not using systemic hormonal contraception, systemic corticosteroid antibiotics, antifungals, or antivirals and without vaginal candidiasis, bacterial vaginosis, sexually transmitted infections, or urinary tract infections were recruited for the study and randomly assigned to use one of the following three gel products in a blinded fashion: a hydroxyethylcellulose (HEC)-based "universal" placebo (96% purified water [USP], 2.7% hydroxyethyl cellulose [NF], 0.85% sodium chloride [NF], 0.1% sorbic acid [NF], sodium hydroxide [NF] to pH 4.4); 6% cellulose sulfate in a gel vehicle that included humectants glycerin and sorbitol, a thickener (carbomer), water, and benzyl alcohol as a preservative; or 4% nonoxynol-9 (Conceptrol) gel, which includes lactic acid, methyl paraben, povidone, propylene glycol, purified water, sodium carboxymethyl cellulose, sorbic acid, and sorbitol. The subjects were instructed to insert their assigned gel twice daily into the vagina for 13.5 days (27 applications of about 3.5 ml each). Tests for seminal fluid were conducted at each visit to confirm that subjects remained abstinent during the course of the study. Vaginal swabs were obtained 14 days prior (V2) and just prior (V3) to the first product application, after the 13th (V4) and last (27th) (V5) applications, and 3 days after (V6) the end of product use (Fig. 1A). V1 was the designation for the initial screening visit, and vaginal swabs were not obtained at that time. At least 8 h separated the insertion of product and sampling. A physician inserted a speculum lubricated with sterile saline water and sampled the midvaginal walls using an Elution Swab system (Copan). The swab was immediately placed into 1 ml of Amies transport medium, frozen upright, and stored at Ϫ80°C. Additional examinations included naked-eye examinations, colposcopies, cervicovaginal lavages, and vaginal smears for Nugent scores (16). Participants were queried about adverse events at each visit.
A total of 60 women were enrolled in the parent study; however, because Institutional Review Board (IRB) approval for this aspect of the study was not received until after the start of subject enrollment, samples for vaginal microbiota analysis were not collected from 10 subjects and only one sample was collected from each of several subjects (4 CS, 3 N-9, and 4 placebo subjects). Another 11 subjects were dropped from the analysis (2 CS, 3 N-9, and 4 placebo subjects deviated from the protocol, and 2 subjects dropped out before randomization). In addition, 35 processed samples were not included in the analysis because the total high-quality read count was less than 500. A total of 146 samples from 35 subjects (13 CS, 12 N-9, and 10 placebo) were included in the analysis ( Table 1). The study was double blinded.
The study was approved by the institutional review boards at Eastern Virginia Medical School, Children's National Medical Center, the University of Maryland School of Medicine, and the University of Pittsburgh.
Nugent Gram stain analysis. The vaginal smears were heat fixed and Gram stained and then blinded and read in random order. A microscopy score of 0 to 10 was assigned by an experienced microbiologist using the standardized method described by Nugent et al. (16). Nugent scores reflect composite scores based on the cellular morphology of the bacteria present in a sample. A score of 0 to 3 is designated normal, 4 to 6 intermediate, and 7 to 10 abnormal and indicative of bacterial vaginosis.
Microbial DNA isolation in the presence of microbicides. Because the nonionic surfactant N-9 and the large polyanion CS were found to interfere with post-DNA-extraction PCR amplification, cell suspensions were first washed prior to lysis to remove the water-soluble microbicide. The use of the flocked nylon swabs facilitated the release of microbial cells in Amies transport media though gentle swirling of the tube (bacteria do not strongly adhere to the nylon swab's axially arrayed fibers and are easily released without the need of a vortex procedure). A total of 1 ml of cell suspension in transport medium was pelleted by centrifugation at 10,000 ϫ g for 10 min and washed twice in 2 ml of 5 mM phosphate-buffered saline (PBS), with bacteria recovered by centrifugation at 10,000 ϫ g for 10 min. After the washes, the cells were resuspended in 2 ml of 5 mM PBS and the suspension was treated with 250 U of cellulase (Sigma-Aldrich) overnight at 37°C. The cells were then washed one more time in 2 ml of 5 mM PBS, centrifuged at 10,000 ϫ g for 10 min, and resuspended in 250 l of 5 mM PBS. Cell lysis and DNA extraction were performed according to Ravel et al. (12), using enzymatic and mechanical lysis. This procedure yielded between 2.5 and 5 g of high-quality, PCR inhibitorfree whole genomic DNA per vaginal swab. The sequence data are available in the NCBI Sequence Read Archive (SRA) (SRA058693) under study accession no. SRP015721.
DNA amplification and pyrosequencing of bar-coded 16S rRNA genes. PCR amplification and 454 pyrosequencing of the V1-V2 hypervariable regions of 16S rRNA genes were performed as previously described (12) using primers 27F and 338R (23).
Sequence analysis. The QIIME software package was used for quality control of sequence reads using the following criteria. Sequences were required to (i) have minimum and maximum lengths of 220 bp and 400 bp; (ii) have an average quality score of q25 over a sliding window of 50 bp (if quality dropped below q25, the read was trimmed at the first base pair of the window and then reassessed for length); (iii) have a perfect match to a bar code sequence; and (iv) include the presence of the 16S primer used for amplification (338r). Sequences were binned by samples using the sample-specific bar code sequences and trimmed by removal of High-quality sequence reads were first dereplicated using 99% similarity and the UCLUST software package, and detection of potential chimeric sequences was performed using the UCHIME component. Chimeric sequences were removed prior to taxonomic assignments. Each processed 16S rRNA gene sequence was classified using the RDP naive Bayesian Classifier (24). RDP classifier quality score filtering was not used, and all reads were classified to the genus or species level as described previously (12). All sequence reads taxonomically assigned to the genus Lactobacillus were identified to the species level using the SpeciateIT software package (http://speciateIT.sourceforge.net). Taxonomic assignments are shown in Table S1 in the supplemental material. Jensen-Shannon divergence. Jensen-Shannon divergence is a measure of dissimilarity between probability distributions that was introduced by Lin (25) to alleviate the limitations associated with relative entropy (13). The Jensen-Shannon divergence between the two discrete probability distributions p ϭ (p 1 , . . ., p n ) and q ϭ (q 1 , . . ., q n ) is defined as follows: where a ϭ (p ϩ q)/2 is the average of p and q, and D KL (x,y) is the relative entropy between x and y. The values of the Jensen-Shannon divergence are normalized to lie between 0 and 1, and its square root is a metric that we refer to as the Jensen-Shannon metric.
The Jensen-Shannon divergence can be expressed using entropy with the following formula: Statistical analyses. Using the methods described previously by Ravel et al. (12), each sample of a vaginal bacterial community, which represents a state of the community at a given time point (community state), was assigned to one of five CSTs (groups of community states with similar microbial species compositions and abundances) using complete linkage hierarchical clustering methods (Fig. 1B).
Modeling the association between CST frequency and product use at V4. CST IV-B, characterized by a lack of Lactobacillus spp. as well as the presence of a variety of strict anaerobes, is often associated with high Nugent scores (14,16). Because high Nugent scores and bacterial vaginosis have been associated with increased risks of transmission and acquisition of several STIs (17,(26)(27)(28), including HIV (17,18), we tested for an association between the frequency of each community state type and product use at visits V2 and V3 and at visit V4. More precisely, the following mixed effects logistic regression model was fitted to the data: y i ϳ Bernoulli͑p i ͒, logit͑p i ͒ ϭ a ϩ b trmt͑i͒ ϩ c stType͑i͒ ϩ d trmt͑i͒,stType͑i͒ ϩ e subjId͑i͒, e s ϳ Norm͑0, b ͒ where Bernoulli(p) is the Bernoulli distribution with mean p, y~Bernoulli(p) means that y is sampled from the distribution Bernoulli(p), and y i is the indicator variable. It is set equal to 0 if the "i"th sample was collected at visits V2 and V3 and set equal to 1 if it was collected at visit V4. trmt(i), stType(i), and subjId(i) are the treatment arm, community state type, and subject identification number (id), respectively, of the community state corresponding to the "i"th sample. The reference class for the treatment variable is the placebo group. The reference class of the community state type variable is CST IV-A. Thus, the baseline group for the model described above consists of placebo samples at CST IV-A. The coefficient a is the mean log odds ratio of samples at visit V4 versus visit V2 or V3 for the baseline group of samples, and b and c are the main effects of treatment and community state type, respectively. The coefficient d captures the interaction between treatment and CST. In order to take into account interactions between samples from the same subject, the model includes a random intercept term (e) that depends on subject id. Estimation of the coefficients of the model was done using Just another Gibbs sampler and jags (29) and the rjags R package (30,31). A noninformative gamma prior with the shape parameter 1 and a scale parameter of 1,000 was used for the precision of the normal distribution, and for the other parameters, the prior was set to the normal distribution with a mean of 0 and a standard deviation of 1,000. Convergence of Markov chains in the model described above was verified using the Gelman-Rubin potential scale reduction factor test (32) at the 1.1 level and the coda R package (33). All Markov chain models were run using 100,000 iterations (with a 100,000-iteration burnin) and a thinning value of 100. The coefficients and their 95% credible intervals are presented in Table S2 in the supplemental material. This analysis revealed that CST IV-B was not associated with product use but that CST IV-A was associated with both N-9 and CS application. To fully characterize this association, we evaluated the ability of the community to not enter CST IV-A by modeling the distances of each community from the center of CST IV-A (see Fig. S1 in the supplemental material). This was accomplished using mixed effects models inferred using Bayesian Markov chain Monte Carlo methods.
Modeling the effect of product application on the community distance from the center of CST IV-A. In order to measure the effect of microbicide treatment on community distance from the center of CST IV-A (see Fig. S1 in the supplemental material), we have used the following double exponential mixed effects model. Community distances were calculated using the Jensen-Shannon divergence: where DExp(, ) is the double exponential distribution with mean and standard deviation , y~DExp(, ) means that y is sampled from the distribution DExp(, ), and y i , trmt(i), visit(i), and subjId(i) are the distances from the center of community state type IV-A, treatment arm, visit, and subject id, respectively, of the community state corresponding to the "i"th sample. The distance from the center of the CST IV-A was computed in the Principal Coordinate Analysis (PCoA) three-dimensional (3D) space, where PCoA was done based on the Jensen-Shannon metric between community states. The center of the CST IV-A was defined as the median position in the PCoA 3D space of CST IV-A community states. The distances to the center of IV-A were computed using the Manhattan metric. The reference class for the treatment variable is the placebo group. The reference class of the visit variable consists of visits V2 and V3. Thus, the baseline group for the above model consists of placebo samples at visits V2 and V3. In the above-described model, the within-sample variance depends on the treatment arm. The coefficient a is the mean distance to CST IV-A for the baseline group of samples, and b and c are the main effects of treatment and visit, respectively. The coefficient d captures the interaction between treatment and visit. In order to take into account interactions between samples from the same subject, the model includes a random intercept term, e, that depends on subject id. Estimation of the coefficients of the model was done using Just Another Gibbs Sampler and jags (29) and the rjags R package (30,31). A noninformative gamma prior with the shape parameter 1 and a scale parameter of 1,000 was used for the precision of the normal distribution, and for the other parameters, the prior was set to the normal distribution with a mean of 0 and a standard deviation of 1,000. Convergence of Markov chains in the above-described The higher frequency of CST IV-A at visits V4 and V5 for treatments with CS and N-9 is notable.
model was verified using the Gelman-Rubin potential scale reduction factor test (32) at the 1.1 level using the coda R package (33). All Markov chain models were run using 100,000 iterations (with a 100,000-iteration burn-in) and a thinning value of 500. The coefficients and their 95% credible intervals are presented in Table S3 in the supplemental material.
Modeling the association between Nugent score category and community state type. The structure of the log-linear model for the Nugent category versus community state type contingency table is as follows: y i ϳ Poisson͑ i ͒, log͑ i ͒ ϭ a ϩ b nugCat͑i͒ ϩ c stType͑i͒ ϩ d nugCat͑i͒,stType͑i͒ where y i , nugCat(i), and stType(i) are the counts, Nugent score category (low ϭ 0 to 3, intermediate ϭ 4 to 6, and high ϭ 7 to 10), and community state type of the "i"th cell of the contingency table (see Table S4 in the supplemental material). Estimation of the coefficients of the model was done using Just Another Gibbs Sampler and jags (29) and the rjags R package (30,31). A noninformative normal distribution with a mean of 0 and a standard deviation of 1,000 was used for all coefficients of the model. Convergence of Markov chains in this model was verified using the Gelman-Rubin potential scale reduction factor test (32) at the 1.1 level and the coda R package (33). The model was run for 50,000 iterations (with a 50,000-iteration burn-in) and a thinning value of 100. Coefficients and their 95% credible intervals for this model are presented in Table S5 in the supplemental material.