Characterization of the first vaginal Lactobacillus crispatus genomes isolated in Brazil

Genome sequencing, assembly, and annotation

In a previous study, vaginal fluid samples were collected from individuals diagnosed as healthy or with BV, with the approval of the ethics committee in research (COEP) of the Federal University of Minas Gerais (protocols ETIC 062/03 and 212/07) (Branco et al., 2010). L. crispatus strains were identified by cellular morphology (Gram-positive bacilli or coccobacilli), biochemistry test (catalase-negative) and 16S-23S rRNA restriction profiling (Branco et al., 2010). The genomes of strains CRI4, CRI8, CRI10 and CRI17 isolated from four healthy patients were sequenced using HiSeq 2500 (Illumina, San Diego, CA, USA) with paired-end libraries of 2 × 150 bp. The sequencing reads quality was examined using FastQC v0.11.8 (Andrews, 2015). The sequencing reads were mapped to the genomes of 46 L. crispatus vaginal isolates to filter out contaminants (Table 1) using bowtie v2 (Langmead & Salzberg, 2012), and the mapped reads were extracted using Samtools v1.7-2 (Li et al., 2009).

Ten different assemblies were generated for each genome using SPAdes v3.14.0 (Bankevich et al., 2012) (assembly 1), Unicycler v0.4.5 (Wick et al., 2017) with SPAdes v3.14.0 (2), ABySS v2.0 (Simpson et al., 2009) (3), MaSuRCA v3.4.0 (Zimin et al., 2013) (4) and an in house pipeline (https://github.com/engbiopct/assembly-hiseq) that generates six other assemblies (5–10). In the in-house pipeline, six strategies combine different software. The best k-mer values were identified using KmerStream v1.1 (Melsted & Halldórsson, 2014). Adapters were removed using AdapterRemoval v2.3.1 (Schubert, Lindgreen & Orlando, 2016). The genome assemblers were Edena v3.131028 (Hernandez et al., 2008) and SPAdes v3.13.0 (Bankevich et al., 2012). The six assembly strategies were: Edena (5), KmerStream and SPAdes (6), KmerStream and SPAdes, using Edena assembly as trusted contigs (7), AdapterRemoval, KmerStream and SPAdes (8), reads processed by AdapterRemoval and the raw reads as input, KmerStream, and SPAdes, using Edena assembly as trusted contigs (9), and AdapterRemoval, KmerStream, and SPAdes, using Edena assembly as trusted contigs (10).

The best assembly for each genome was determined using QUAST v5.0.2 (Gurevich et al., 2013). Then, the genome’s paired-read sequencing data was used for scaffolding with SSPACE v3.0 (Boetzer et al., 2011) and contig extension and gap filling with GapFiller v1.1.1 (Boetzer & Pirovano, 2012). Finally, more gaps were closed using contigs from the other nine assemblies and the chromosome of L. crispatus strain AB70 (CP026503.1) (Chang et al., 2019) as a reference, using GFinisher (Guizelini et al., 2016). The presence of plasmids was investigated using PlasmidFinder 2.1 (Carattoli et al., 2014). The four genomes were identified from the species L. crispatus using the Type (Strain) Genome Server (Meier-Kolthoff & Göker, 2019). The assemblies completeness was evaluated using BUSCO v4.0.6 (Seppey, Manni & Zdobnov, 2019), based on the presence of 402 single-copy orthologous genes shared within Lactobacillales. The genomes were annotated using Prokka v1.11 (Seemann, 2014). The GenBank accession numbers of the genomes from strains CRI4, CRI8, CRI10 and CRI17, isolated from healthy patients, are JABERN01, JABERO01, JABERP01 and JABERQ01, respectively (Table 2). Type Strain Genome Server (Meier-Kolthoff & Göker, 2019) was used to confirm the taxonomic classification of the 50 samples as L. crispatus strains.

Criteria of public genomes selection

We included in the analysis genomes available in public databases using the following criteria: (i) vaginal isolate, (ii) the metadata explicitly informs the health condition of the individual, and/or the microbiome classification. A total of 46 samples were classified as belonging to “healthy” (26) or “BV” (20) condition groups and were used in previous studies (Hemmerling et al., 2010; Ojala et al., 2014; Abramov et al., 2014; Abdelmaksoud et al., 2016; Eslami et al., 2016; Dols et al., 2016; Chang et al., 2019; Clabaut et al., 2019, 2020; Van der Veer et al., 2019). Table 1 shows the genomes list, including the microbiome classification, when available, and the terms in metadata used to classify the sample.

Probiotic features

Bacteriocins and linear azol(in)e-containing peptides (LAPs) were predicted using BAGEL4 (Van Heel et al., 2018). The enzymes involved in the production of L- and D-lactate and hydrogen peroxide were identified using KEGG Mapper/BlastKOALA, under the pyruvate metabolism pathway. Adhesins were predicted using eggNOG-mapper v2 (Huerta-Cepas et al., 2017). Protein IDs were identified using BLASTp (Camacho et al., 2009) using the GenBank non-redundant (nr) database, selecting hits with 100% identity and coverage. Pathways involved in the biosynthesis of antimicrobial drugs with clinical importance (Chokesajjawatee et al., 2020) were predicted using KEGG Mapper.

Safety assessment of the strains for probiotic applications

Detection of plasmids, insertion sequences, prophages, and CRISPR-Cas elements was performed using PlasmidFinder (Carattoli et al., 2014), ISEScan (Xie & Tang, 2017), PHASTER (Arndt et al., 2016), and CRISPRCasFinder (https://crisprcas.i2bc.paris-saclay.fr/CrisprCasFinder/Index), respectively. The domains of the identified Cas proteins were predicted using InterProScan (Jones et al., 2014) and NCBI’s Conserved Domain Database (Lu et al., 2020). Multiple alignments of sequences of interest were performed using the Clustal Omega web service (Larkin et al., 2007). Local alignment of sequences of interest across genomes was performed using BLASTn (Camacho et al., 2009) implemented in PATRIC (Davis et al., 2020).

Virulence factor genes were detected using the databases VFDB (Liu et al., 2019) and Ecoli_VF (https://github.com/phac-nml/ecoli_vf), while antimicrobial resistance genes were detected using the databases ARG-ANNOT (Gupta et al., 2014), CARD (Alcock et al., 2020), MEGARes (Doster et al., 2019), NCBI AMRFinderPlus (Feldgarden et al., 2019) and ResFinder (Bortolaia et al., 2020). The screening using these databases was performed using ABRicate (https://github.com/tseemann/abricate) with default parameters.

The presence of the toxins hemolysins, enzymes involved in the synthesis of biogenic amines and other undesirable genes, as listed by Chokesajjawatee et al. (2020), were manually screened using KEGG Mapper (Kanehisa & Sato, 2020) implemented in BlastKOALA v2.2 (Kanehisa, Sato & Morishima, 2016).

In vitro detection of phages

One representative strain (L. crispatus CRI4) was subjected to the induction of lysogenic bacteriophages according to previous studies (Kiliç et al., 1996; Raya & H’bert, 2009). A total of 200 µL of the culture was initially inoculated in 5 mL of Man–Rogosa–Sharpe medium. The optical density (OD) was measured until reaching DO₆₀₀ = 0.2. Then, 0.4 µg/mL of Mitomycin C (Sigma, St. Louis, MO, USA) was added in the culture and incubated at 37 °C overnight. Afterward, the supernatant was collected and filtered in sterile 0.22 µm membrane. A 10 µL sample of the filtered lysate was applied to a 200 mesh grid at the UFMG electron microscopy center (CM-UFMG) and visualization was performed in Tecnai G2-12 Transmission Electron Microscope, SpiritBiotwin FEI, 120 kV.

Comparative analyses using public genomes

We performed comparative analysis including the for strains from Brazil and 46 public genomes (Table 1) to reconstruct their phylogeny and to identify adaptations that could be related to the colonization of the vaginal niche by testing for association between gene presence/absence and features of interest, and detection of positive selection in protein-coding genes.

For phylogenomic analysis, L. helveticus DSM20075 genome (GenBank accession ACLM01) was used as an outgroup, the conserved genes across all 51 genomes were estimated by Roary v3.6.0, and their nucleotide sequences were aligned using MAFFT (Katoh et al., 2005) implemented in Roary. The alignment was used as input for IQ-Tree v1.6.12 (Nguyen et al., 2015) for phylogenetic inference using the Maximum Likelihood. The confidence values were estimated using 1,000 rounds of bootstrapping. The tree was edited using iTOL (Letunic & Bork, 2019).

We tested the associations of gene presence/absence with a health condition and geographical locations suggested by the phylogenetic tree. A GWAS based on gene presence/absence was performed using Scoary v1.6.16 (Brynildsrud et al., 2016). Scoary estimates association by pairwise comparisons on a phylogeny (Maddison, 2000) to correct population structure and permutation. The input was a gene presence-absence matrix from the 50 L. crispatus genomes estimated using Roary and a phylogenetic tree generated by IQ-Tree v1.6.12, utilizing the core gene alignment calculated by Roary, and a matrix containing the presence-absence of the features across the samples.

A genome-scale positive selection analysis was performed using POTION v1.1.2 (Hongo et al., 2015). To generate the input, FastOrtho (https://github.com/PATRIC3/FastOrtho) was used to identify ortholog groups across the 50 L. crispatus genomes. The file containing the orthologous group’s information and multifasta files containing nucleotide sequences of protein-coding genes were used as input for POTION v1.1.2. The genome-scale positive selection analysis used site tests with the models M1 and M2, and M7 and M8 (Yang & Nielsen, 2002). The POTION configuration file is available as Data S1. The function of the identified proteins was annotated using eggNOG-mapper (Huerta-Cepas et al., 2017), the subcellular localization using SufG+ v1.2.1 (Barinov et al., 2009), and the GenBank protein ID using BLASTp (Agarwala et al., 2016).

Genome sequencing and taxonomy

The four sequenced genomes were assembled as drafts, and no plasmid was found. Table 2 shows the statistics of genome assembly and annotation. The genomes completeness ranged from 99% to 99.5%, while the reference strain AB70, available as a complete genome sequenced using PacBio RS II platform (Chang et al., 2019), had its completeness estimated as 99%. The four sequenced genomes and the 46 public genomes were classified as L. crispatus by TYGS, with dDDH > 70% and G+C content divergence of less than 1% to the strain JCM 1185^T (Data S2).

Probiotic features

Features associated with probiotic effects in the four strains from Brazil are shown in Table 3. We identified genes involved in the biosynthesis of D-lactate (1), L-lactate (3), hydrogen peroxide (1), bacteriocins (9), LAPs (1) and adhesins (10, five classes) across the four genomes. No pathway for biosynthesis of antimicrobial drugs of clinical importance was found.

Safety assessment of the strains for probiotic applications

Features associated with safety in the four strains from Brazil are shown in Table 4. About mobile elements, no plasmid was predicted, as previously stated. A total of 131 to 184 IS from 14 known families were identified in the four genomes, and one new family was predicted across CRI8 (2 copies), CRI10 (2) and CRI17 (1) (Data S3). The new IS has a size of 2,088 bp and harbors two genes coding a site-specific integrase-resolvase (IS607-like family, GenBank protein ID WP_005728427.1) and a transposase (IS605 family, AZR15009.1) (Data S4). The multiple alignments of the five copies show the IS is fragmented in the contig 44 of strain CRI10, and the difference among the complete sequences is a SNP (T to G) in position 1,350 (Data S5). A BLASTn against the four strains from Brazil and other 123 public L crispatus genomes database identified hits with ≥99% identity and ≥98% coverage with 42 genomes from strains isolated from the female human urogenital tract, 39 of them public genomes (Data S3).

Five or six prophages were predicted across the four genomes, classified as questionable (score 70–90) or incomplete (score < 70), with most of them close to contig ends. In strains CRI8 and CRI17, some of the predicted prophages harbored the new IS (Fig. 1; Data S6). We identified CRISPR-Cas systems subtype II-A in the four strains and subtype II-C only in strain CRI10 (Table 4; Data S7). A closer examination of the subtype II-A system reveals that all four genomes contain CRISPR loci and Cas proteins-encoding genes cns2, cas2, cas1, and cas9. However, a transposase is inserted in csn2 in each strain, resulting in two gene fragments. The gene cas9 is also fragmented in all strains, with a transposase inserted between fragments in strains CRI4 and CRI8. CRI17 has an extra copy of cas9. The subtype II-C system is represented in CRI10 by a single cas9 gene (Fig. 2). All the cas9 CDSs lack one or more domains.

No virulence or antimicrobial resistance gene was predicted using ABRicate. With the search for virulence and undesirable genes using KEGG Mapper/BlastKOALA, we identified two hemolysins (tlyC and hlyIII), one bile salt hydrolase (cbh), and one enzyme involved in the biosynthesis of putrescine (Ornithine decarboxylase, odcI). Manual screening of Prokka annotation showed a third hemolysin, Hemolysin A. The search for toxin biosynthesis identified an “Ornithine decarboxylase” (odcI) (WP_005727730.1) involved in the synthesis of the biogenic amine putrescine, in all for strains.

In vitro detection of phages

Several rod-shaped particles of varying length and thickness were observed by electron microscopy, some of them similar to Myoviridae bacteriophages. This pattern was repeated throughout the slide (Fig. 3).

Comparative analyses using public genomes

The strains did not cluster according to a health condition. However, strains from Brazil were the closest to strains from the Netherlands, even in different clusters (Fig. 4). The GWAS based on gene presence/absence implemented in Scoary did not identify an association with a health condition (p > 0.05). Due to the clustering of strains from Brazil and the Netherlands, we tested the association with these two geographical locations and the result was also negative. A total of 8 protein-coding genes were identified as under positive selection (q < 0.05). After manual curation for false positives caused by alignment artifacts, five genes were obtained: a surface exposed “S-layer protein precursor,” three phage related proteins, and a hypothetical protein (Table 5; Data S8).