The Cervicovaginal Microbiota-Host Interaction Modulates Chlamydia trachomatis Infection

The vaginal microbiota is believed to protect women against Chlamydia trachomatis, the etiologic agent of the most prevalent sexually transmitted infection (STI) in developed countries. The mechanism underlying this protection has remained elusive. Here, we reveal the comprehensive strategy by which the cervicovaginal microbiota modulates host functions to protect against chlamydial infection, thereby providing a novel conceptual mechanistic understanding. Major implications of this work are that (i) the impact of the vaginal microbiota on the epithelium should be considered in future studies of chlamydial infection and other STIs and (ii) a fundamental understanding of the cervicovaginal microbiota’s role in protection against STIs may enable the development of novel microbiome-based therapeutic strategies to protect women from infection and improve vaginal and cervical health.

reported procedures (1). Briefly, mid-vaginal Eswabs, stored in Amies transport medium (Copan) were thawed on ice and vortexed briefly. A 0.5 ml aliquot of the cell suspension was transferred to a FastPrep Lysing Matrix B (MP Biomedicals) tube containing 0.5 ml of PBS (Invitrogen). A cell lysis solution containing 5 µl lysozyme (10 mg/ml, EMD Chemicals), 13 µl mutanolysin (11,700 U/ml, Sigma-Aldrich), and 3.2 µl lysostaphin (1 mg/ml, Ambi Products) was added and samples were incubated at 37°C for 30 min. Then, 10 µl Proteinase K (20mg/ml, Invitrogen), 50 µl 10% SDS (Sigma-Aldrich), and 2 µl RNase A (10 mg/ml, Invitrogen) were added and samples were incubated at 55°C for 45 min. Cell lysis by mechanical disruption was performed on a FastPrep homogenizer at 6 m/s for 40 s, and the lysate was centrifuged on a Zymo Spin IV column at 7000 x g for 1 min. (Zymo Research). Lysates were further processed on the QIAsymphony platform using the QS DSP Virus/Pathogen Midi Kit (Qiagen) according to the manufacturer's recommendation. DNA quantification was carried out using the Quant-iT PicoGreen dsDNA assay (Invitrogen). Sequencing libraries were constructed using either 1 step or 2 step PCR (details in supplementary methods) (3,4). Composition of the vaginal microbiota was assessed by metataxonomics and amplicon sequencing of the 16S rRNA gene V3-V4 hypervariable regions (5,6). Libraries were sequenced on an Illumina MiSeq instrument using 600 cycles generating paired-end 300 bp sequence reads. The sequences were de-multiplexed using a dual-indexing approach (3).

1-
Step PCR library construction. Primer sequences ranged from 90-97 bp and amplification was performed using Phusion Taq Master Mix (1X, ThermoFisher) with 3% DMSO, 0.4 mM of each primer, and 5 mL of genomic DNA. A standard volume of genomic DNA was used for each library. Cycling conditions were as follows: initial denaturation at 98°C for 30s, 30 cycles of denaturation at 98°C for 15s, annealing at 58°C for 15s, and elongation at 72°C for 15s, followed by a final elongation step at 72°C for 60s (3).

2-
Step PCR library construction. This library preparation (4) is modified from an Illumina protocol (https://support.illumina.com/downloads/16s_metagenomic_sequencing_library_preparation.htm l). The 16S rRNA gene V3-V4 region from genomic DNA was targeted using bacterial primers 338F and 806R combined with a heterogeneity spacer of 0-7 bp, and Illumina Sequencing Primers. A single PCR master mix containing an equal ratio of all primers, was used. Each PCR reaction contained 1X Phusion Taq Master Mix (ThermoFisher), Step 1 Forward and Reverse primers (0.4 mM each), 3% DMSO, and 5 mL of genomic DNA. The following cycling conditions were used for PCR amplification: initial denaturation at 94°C for 3 min, 20 cycles of denaturation at 94°C for 30s, annealing at 58°C for 30s, and elongation at 72°C for 60s, and a final elongation step at 72°C for 7 min. The resultant amplicons were diluted 1:20, and 1 mL was used in the second step PCR which introduced an 8 bp dual-index barcode to the 16S rRNA gene amplicons and the flow cell linker adaptors using primers containing a sequence that anneals to the Illumina sequencing primer sequence introduced in step 1. Each primer was added to a final concentration of 0.4 mM in each sample specific reaction, along with Phusion Taq Master Mix (1X) and 3% DMSO. Phusion Taq Polymerase (ThermoFisher) was used with the following cycling conditions: initial denaturation at 94°C for 30s, 10 cycles consisting of denaturation at 94°C for 30s, annealing at 58°C for 30s, and elongation at 72°C for 60s, followed by a final elongation step at 72°C for 5 min.
Microbiota bioinformatic analysis. The resulting forward and reverse sequence fastq files were split by sample using the QIIME-dependent script split_sequence_file_on_sample_ids.py, and primer sequences were removed using TagCleaner (version 0.16) (7). Further processing followed the DADA2 Workflow for Big Data and dada2 (v. 1.5.2) (https://benjjneb.github.io/dada2/bigdata.html) (8,9). Taxonomy was assigned to each read using a novel fifth-order Markov Chain Monte Carlo taxonomic classifier (available at http://ravellab.org/speciateit) and taxa frequencies normalized to total per-sample read counts. Community state types (CSTs) were identified by calculating the Jensen-Shannon divergence among samples, followed by hierarchical clustering with complete linkage (1). Clusters were assigned to one of the five CSTs described previously (1,5).

Antimicrobial resistance test. Minimal inhibitory concentration (MIC) for 22 vaginal bacteria
strains was determined by broth microdilution using azithromycin and doxycycline in TSB and NYCIII medium (Sigma-Aldrich). Using a previously published protocol, bacteria were exposed to 12 dilutions ranging from 0.03-640 mg/L and the inhibitory concentration determined by OD (10).
Lactic acid solutions pH adjusted by increasing lactic acid concentration. D(-), D/L (racemic mixture of each isomer) or L(+) lactic acids (Sigma-Aldrich) were added to A2EN complete medium at concentrations of 15 mM, 22.5 mM and 28 mM. pH was measured and when necessary adjusted to pH values of 7, 5.5 and 4, respectively using 100mM lactic acid.
Hydrochloric acid (HCl) (Sigma) was added at concentrations of 15 mM, 17 mM and 19 mM, and adjusted to pH 7, 5.5 and 4, respectively using 100 mM HCl.
Lactic acid solutions pH adjusted from a 1% lactic acid with NaOH. A 30% (w/v) stock solutions of each D(-), D/L or L(+) lactic acid were made and used to prepare 1% lactic acid in A2EN complete medium (pH ~3) for each experiment. The pH of the media was adjusted to pH 7 and 4 using 1M NaOH. A 1% (w/v) stock solution of HCl was made, followed by a 30X dilution in A2EN complete medium (pH~3) and adjusted to pH 7 and 4 using 1M NaOH.
A2EN and VK2 cell viability assay. A2EN and VK2 cell viability was performed using the Viability/Cytotoxicity Assay Kit for Animal Live and Dead cells (Biotium) as per the manufacturer's instructions. Briefly, epithelial cells were either exposed to media containing lactic acid or to different culture supernatants for 30 min or to 20% diluted culture supernatants for up to 24h, rinsed and incubated with 4 µM Calcein AM and 1 µM EthDIII for 1 h in the dark.
Images were taken using the FITC and RFP filters on a Zeiss Axio Imager Z1(Zeiss). Manual analysis was performed to determine the percentage of live (green) relative to dead (red) cells within each sample.
Mucus penetration imaging. mCherry labeled Chlamydia trachomatis serovar L2 were exposed to the cervicovaginal mucus and analyzed in a manner similar to fluorescent HIV (11). Briefly fluorescent bacteria or beads were added at 5% (vol/vol) to 20µl of CVM placed in a custommade glass chamber and incubated for 1 h at 37°C prior to microscopy. An aliquot of CVM was titrated to pH 6.8 to 7.1 using small volumes (~3% [vol/vol]) of 3 N NaOH and then analyzed.
Using an electron multiplying charge-coupled-device (EMCCD) camera (Evolve 512; Photometrics, Tucson, AZ) mounted on an inverted epifluorescence microscope (AxioObserver D1; Zeiss, Oberkochen, Germany) the translational motions of the particles were recorded. Using MetaMorph imaging software (Molecular Devices, San Jose, CA) videos (512 X 512 pixels, 16-bit image depth) were captured and the tracking resolution was determined by tracking the displacements of particles. Trajectories of ≥40 particles per frame on average were analyzed for each sample, typically corresponding to >100 individual particle traces throughout the video.
Sample selection for small RNA-sequencing.
Samples for small RNA-seq were selected from vaginal microbiota profiles previously characterized by metataxonomics analysis (16S rRNA gene sequencing) (2). Samples were selfcollected daily for 10 weeks by 135 reproductive-age women and dropped off weekly to the study site. A clinical visit was performed at enrollment, week 5 and week 10 of the study.  Table S1B). free water to accommodate 50-100 ng total RNA input. The RNA template was denatured for 2 min at 70°C, then 1µl diluted 3' adaptor, 1µl RNase Inhibitor, 1µl Enzyme 1 and 5µl Buffer 1 were added to the template, mixed, and incubated at 28°C for 1 hour followed by incubation at 65°C for 20 min. Following this step, 4µl nuclease-free water, 1µl Buffer 2, 1µl RNase Inhibitor, and 2µl Enzyme 2, was added to the RNA template and the 3' adaptor mixture. The diluted 5' adaptor was denatured for 2 min at 70°C, then 2µl was added to the mixture and incubated at The beads were magnetized for 4 min, then the supernatant containing the library was transferred to a new tube where 144µl beads were added and incubated for 10 min to bind DNA. The beads were magnetized again for 4 min, the supernatant discarded, and the beads washed twice with 500µl 70% ethanol. After the wash, the beads were left to air-dry before resuspending in 17µl nuclease-free water for 2 min. The solution was re-magnetized and 15µl was transferred as the small RNA-seq library. . Strand-specific genomic feature overlaps were counted using HTSeq version 0.5.3p3 (16) with default parameters (mode=union, minaqual=0, stranded='yes') and the iGenomes annotation as above (Table S3 and  Mix, 4µl 5X miScript HiSpec Buffer and 5µl nuclease-free water. The reaction mixture was incubated at 37°C for 1h, then heated to 95°C for 5 min to stop the reaction per the manufacturer's protocol (17,18). The resulting cDNA was diluted 1:10 in water before use in the subsequent qPCR assay. qPCR reactions were performed using 1µl diluted cDNA template, and the miScript SYBR Green PCR Kit by mixing 5µl 1X SYBR Green Mastermix and 1µl each of Universal and miRNA-specific primers in a 10µl reaction volume (QIAGEN #218073 and #218160, miR-specific primers: MS00031549 (miR-193b-3p, 5'AACUGGCCCUCAAAGUCCCGCU) and MS00033740 (RNU6-2-11) (QIAGEN). qPCR was carried out on a 7900HT thermocycler (ThermoFisher) at the following cycle conditions: 95°C for 15 min, then 40 cycles of: 94°C for 15s, 55°C for 30s and 70°C for 30s. The ΔΔCt method was used to compute the ΔCt between miR-193b and RNU-6 within a sample, and then the ΔCt between Lactobacillus spp. CS or lactic acid and cell culture medium (19). The ΔΔCt standard deviation was computed using σΔΔCt =(σ,ΔCt Lactobacillus spp. 2 + σ,ΔCt non-Lactobacillus spp. 2 ) 1/2 , where σ is the standard deviation. A two-tailed t-test was applied to each ΔΔCt, with p<0.05 considered significant.
Identification of Lactobacillus spp.-associated miRNAs using Random forest. A Random forest model was applied to select miRNAs that best predict Lactobacillus spp. relative abundance by utilizing a combination of the R software packages rfPermute (version 2.0.1, (20)), randomForest (version 4.6-12, (21)) and custom subroutines (available in the R scripts at https://github.com/ravel-lab/smith_thesis_2017/tree/feature/manuscript-2.0/AnalysisPipeline/Scripts). Each sample's Lactobacillus spp. relative abundance (proportion of reads mapping to Lactobacillus crispatus, Lactobacillus jensenii, Lactobacillus iners, and Lactobacillus gasseri, available in SRA under project PRJNA208535) were used as the response variable to determine the most important miRNAs that predicted the community state while the log2 transformed normalized miRNA counts were used as predictors. The training set consisted of 70% of available data while model performance was assessed using the remaining 30% of the held-out data. To increase the confidence of feature calls, miRNAs with zero counts in any sample were excluded, as zero miRNA counts could be due to under sampling. Each model underwent 10-fold cross-validation, with 500 permutations to determine the null distribution for p-value calculation. Default parameters for ntree (500) and mtry (1558 input features/3=519) were used. The algorithm accounted for non-independent samples that originated from the same subject by evenly splitting each cross-fold iteration among subjects using a custom script (available at https://github.com/ravel-lab/smith_thesis_2017/tree/feature/manuscript-2.0/AnalysisPipeline/Scripts). Importance metrics and p-values were calculated based on rfPermute and randomForest R packages (20,21). The increase in mean squared error and increase in node purity were used to rank features (22,23). Statistically significant features were defined as features with p-value<0.05 for any importance metric within a model result.

RNA-seq differential expression and pathway analysis.
All analysis scripts can be found at https://github.com/ravel-lab/smith_thesis_2017/tree/feature/manuscript-  (Table S3). All of the top 10 most abundant unaligned reads from all samples were of human origin.
The R package edgeR, version 3.10.5, was used to compute pairwise differential expression between combinations of each exposure time and BCS vs. cell culture medium (25). Negative binomial dispersion was estimated for samples passing QC by applying the estimateDisp function available through edgeR. Samples were normalized using the calcNormFactors function (26). Reads were fit to a negative binomial generalized linear model using the glmFit function available in edgeR, using the sample's treatment as the design matrix. Differential expression using edgeR's likelihood ratio test was computed for each gene using the glmLRT function.
Genes with an average log counts per million (logCPM) >1 , log2-transformed Fold Change (logFC) > 1 and false discovery rate (FDR) < 0.01 (27,28)(87) were considered differentially expressed between treatments. The mean logCPM as calculated by edgeR is the log2 counts per million reads, averaged over all the libraries, while log2FC is the coefficient of the Generalized Linear Model used by edgeR (25). Gene expression plots were created using ggplot2 (version 2.2.1) and custom scripts (29).
Differentially expressed genes for each comparison were used to generate pathway enrichment scores using Ingenuity Pathway Analysis (IPA) (QIAGEN, build version 439932M, content version 33559992). IPA computes a pathway activation score (z-score) for each pathway comparison based on inferred expression directionality using the logFC of each gene (29).
Determining gene ontology for Lactobacillus spp.-associated miRNAs. Experimentally validated miRNA targets were identified using the "strong evidences" list from miRTarBase, Release 6.0 (Sept 15, 2015) (30). The Gene Ontology DIRECT process terms from DAVID (release 6.8, May 2016, (31)) were mapped to experimentally validated miRNA targets. The proportions of targets for each GO DIRECT term were computed for each miRNA. The miRNA with the highest overall expression and most correlated with Lactobacillus spp. relative abundance was chosen for further experimentation.
Scratch assay using bacterial culture supernatants. VK2 epithelial cells, seeded at 7.5x10 4 cells/well, were grown to confluence, starved for 24h in base medium at which time a scratch made with a 1 ml sterile pipette tip. VK2 cells were then exposed for up to 22h to either lactic acid medium (0.06% of either D(-) or L(+) lactic acid in VK2 cell culture medium) or culture supernatants (created as previously described) diluted to 20% (v/v) using complete VK2 cell culture medium. Phase contrast images were taken at 100X using a Zeiss Primovert microscope (Zeiss) at 0 and 13h post-culture supernatant exposure. The proportion of cells occupying the scratched area relative to time 0h was used to quantify migration (ImageJ software (version 1.50i, (32)). Total RNA from exposed cells were then extracted by first adding 300 µl RNAlater and mechanically detaching cells from the plate, then using the total RNA extraction method described above. To monitor active DNA synthesis for cell proliferation, EdU (5-ethynyl-2'deoxyuridine) assay was carried out per manufacturer's instructions (ThermoFisher). Briefly, cell nuclei were stained using Hoechst 33342 (1:1,000) and imaged at 20X using a Zeiss Axio Imager Z1 fluorescence microscope (Zeiss). The amount of DNA synthesis was calculated using CellProfiler (version 2.2.0 rev 9969f42) (33) by counting the number of green nuclei (EdU stained) relative to blue nuclei (DAPI stained) in five fields per duplicate coverslip.
Protein band intensities were measured using ImageJ Studio software, then CCND1 intensity values were normalized to b-actin loading control value.