Lateral Gene Transfer Shapes Diversity of Gardnerella spp.

Abstract

Gardnerella spp. are pathognomonic for bacterial vaginosis, which increases the risk of preterm birth and the transmission of sexually transmitted infections. Gardnerella spp. are genetically diverse, comprising what have recently been defined as distinct species with differing functional capacities. Disease associations with Gardnerella spp. are not straightforward: patients with BV are usually infected with multiple species, and Gardnerella spp. are also found in the vaginal microbiome of healthy women. Genome comparisons of Gardnerella spp. show evidence of lateral gene transfer (LGT), but patterns of LGT have not been characterized in detail. Here we sought to define the role of LGT in shaping the genetic structure of Gardnerella spp. We analyzed whole genome sequencing data for 106 Gardnerella strains and used these data for pan genome analysis and to characterize LGT in the core and accessory genomes, over recent and remote timescales. In our diverse sample of Gardnerella strains, we found that both the core and accessory genomes are clearly differentiated in accordance with newly defined species designations. We identified putative competence and pilus assembly genes across most species; we also found them to be differentiated between species. Competence machinery has diverged in parallel with the core genome, with selection against deleterious mutations as a predominant influence on their evolution. By contrast, the virulence factor vaginolysin, which encodes a toxin, appears to be readily exchanged among species. We identified five distinct prophage clusters in Gardnerella genomes, two of which appear to be exchanged between Gardnerella species. Differences among species are apparent in their patterns of LGT, including their exchange with diverse gene pools. Despite frequent LGT and co-localization in the same niche, our results show that Gardnerella spp. are clearly genetically differentiated and yet capable of exchanging specific genetic material. This likely reflects complex interactions within bacterial communities associated with the vaginal microbiome. Our results provide insight into how such interactions evolve and are maintained, allowing these multi-species communities to colonize and invade human tissues and adapt to antibiotics and other stressors.

Introduction

Gardnerella spp. are Gram-variable, facultative anaerobes found in the vaginal microbiome of healthy women (Schellenberg et al., 2017). However, Gardnerella spp. are also associated with vaginal dysbiosis and bacterial vaginosis (BV), which is characterized by symptoms such as abnormal vaginal discharge, malodor, and pain (Hilbert et al., 2017). Additionally, BV can increase the risk of preterm birth and transmission of sexually transmitted diseases, including HIV (Hilbert et al., 2017).

Gardnerella spp. were previously considered a single species (i.e., Gardnerella vaginalis), but more recent research proposes that what was previously designated G. vaginalis in fact comprises several distinct species (Vaneechoutte et al., 2019). Biochemical tests (Piot et al., 1984) and phylogenetic methods (Ahmed et al., 2012; Cornejo et al., 2018) have delineated at least four distinct clades/species. While Gardnerella spp. are known to be associated with BV, clear consensus is lacking for the clades or combinations of clades that are most consequential for disease (Janulaitiene et al., 2017; Hill et al., 2019). Women with BV are often infected with strains from two or more clades (Hilbert et al., 2017). Recently, a study used average nucleotide identity and digital DNA-DNA hybridization to update the description of G. vaginalis, describe 3 new species, and 9 different genomospecies within “Gardnerella spp.” (Vaneechoutte et al., 2019). Researchers have tried to identify which groups of isolates are more pathogenic than others; however, this complex relationship remains to be fully defined (Harwich et al., 2010; Balashov et al., 2014; Janulaitiene et al., 2017; Hill et al., 2019). Our analyses build on previous work that used comparative genomics approaches to predict potential functional differentiation of this multi-species community (Cornejo et al., 2018). Here we investigate how these functional differences have evolved and are maintained in Gardnerella spp. populations.

We previously found evidence suggesting lateral gene transfer (LGT) is frequent among Gardnerella spp. (Devault et al., 2017). We and others found that LGT appeared to be structured by species, and that despite abundant intergenomic recombination, the species have maintained their genetic distinctiveness (Ahmed et al., 2012). Here we sought to characterize LGT in detail in a larger sample and to investigate whether patterns of recombination differ between species. The mechanism(s) of recombination in Gardnerella spp. are not known, but previous investigators found four predicted competence genes in a sample of three isolates (Yeoman et al., 2010). Additionally, prophage genes have been identified in Gardnerella spp. (Malki et al., 2016). As part of our investigation of LGT in Gardnerella spp. we further sought to identify and characterize genes with a potential role in shaping LGT, including competence genes, phage, and restriction modification systems, in a larger and more diverse sample of isolates. Elucidating how DNA is exchanged in bacterial populations helps to illuminate how clinically relevant traits, such as virulence and antibiotic resistance, evolve in these complex communities.

Materials and Methods

Data Set

We obtained whole genome sequencing reads or de novo assemblies (when reads were unavailable) for 97 isolates of Gardnerella spp. from NCBI. We additionally sequenced nine clinical Gardnerella spp. isolates. The accessions for the nine newly sequenced isolates are listed in Table S1. Accession numbers and available clinical data for all 106 isolates are listed in Table S2.

Bacterial Growth and Isolation

Strains were streaked on human blood bilayer with Tween (HBT) agar (Bd Diagnostic Systems) and incubated at 37°C in 5% CO₂ for 3–4 days. Growth was subcultured in 6 mls of brain heart infusion (BHI) broth (Teknova) + 10% fetal bovine serum (FBS) + 5% Fildes enrichment + 1 μg/mL amphotericin for 48 h at 37°C in 5% CO₂.

DNA Extraction

We used the gBac Mini gDNA Bacteria Kit (IBI, Lot No JM14117) for DNA extraction of 9 clinical isolates with the following modifications: the entire 6 ml culture was centrifuged for 5 min at 5,000 g to pellet cells and incubated with lysozyme for 1 h.

Library Preparation and Sequencing

For the 9 new clinical isolates, libraries were prepared using a modified Nextera protocol as described by Baym et al. (2015) with a reconditioning PCR to minimize chimeras with fresh primers and polymerase for an additional 5 cycles and a bead based size selection (650 b). Libraries were sequenced on an Illumina HiSeq 2500 (paired-end, 150 bp).

De novo Assembly and Annotation

We used the iMetAMOS (Koren et al., 2014) pipeline to compare de novo assemblies from SPAdes (Bankevich et al., 2012), MaSurCA (Zimin et al., 2013), and Velvet (Zerbino and Birney, 2008). KmerGenie (Chikhi and Medvedev, 2014) was used to select kmer sizes for assembly. Quality of reads and assemblies were assessed using FastQC (Andrews, 2010), QUAST (Gurevich et al., 2013), REAPR (Hunt et al., 2013), LAP (Ghodsi et al., 2013), ALE (Clark et al., 2013), FreeBayes (Garrison and Marth, 2012), and CGAL (Rahman and Pachter, 2013), and contamination was detected with Kraken (Wood and Salzberg, 2014). De novo assemblies were annotated with Prokka (Seemann, 2014).

Core Genome Identification and Alignment

We used PIRATE (Bayliss et al., 2019) to identify the core and pan genomes for newly sequenced and publicly available Gardnerella spp. genomes (Table 1). Using PIRATE, we clustered orthologous gene families using an amino acid identity threshold ranging from 50 to 100% to obtain a clearer understanding of the breadth of diversity across the pangenomes of Gardnerella spp. We used the pangenome information from PIRATE and concatenated core genome alignments of single copy genes at 100% frequency for the entire dataset as well as individual species. In addition, we used Roary (Page et al., 2015) to identify core and pan genomes at an amino acid percent identity threshold of >75%. This allowed us to compare gene homologs at a consistent threshold across all genes.

Recombination Detection

To identify recombination events between the major clades/species of Gardnerella spp., we used FastGEAR (Mostowy et al., 2017) on a concatenated core genome aligned with MAFFT (Katoh and Standley, 2014). Briefly, FastGEAR uses a Hidden Markov Model approach to cluster isolates into lineages, detect ancestral and recent recombination, and measure the statistical strength of the recombination events. We used Gubbins (Croucher et al., 2015) to identify recombination events within the core genome alignments of clade 1 (Gardnerella vaginalis) and clade 2 (Gardnerella piotii). Briefly, Gubbins identifies recombination by using spatial scanning statistics to identify loci with elevated single nucleotide polymorphism (SNP) densities. To account for differences in sample size in G. vaginalis and G. piotii, we subsampled G. vaginalis to the size of G. piotii and used Gubbins to identify recombination in the subsampled dataset. We calculated the proportion of sites affected by recombination per isolate within G. vaginalis and G. piotii and compared the means using a Mann-Whitney-Wilcoxon test (Mann and Whitney, 1947).

Codon Usage

We calculated codon adaptation index (CAI) for the core genes of G. vaginalis and G. piotii using the EMBOSS CUSP program (Rice, 2000). We compared the means of CAI values using a Student's t-test.

Phylogenetic Network Inference

Using the concatenated core genome alignment constructed from PIRATE output, we inferred a phylogenetic network using SplitsTree 4 (Huson and Bryant, 2006) of all Gardnerella isolates, G. vaginalis, and G. piotii isolates.

Maximum Likelihood Phylogenetic Inference

We performed maximum likelihood phylogenetic inference on the concatenated alignment of core genes using RAxML v 8.2.3 (Stamatakis, 2014) with the GTR model of nucleotide substitution and gamma distribution of rate variation. Twenty trees were estimated for the alignment, and the tree with the maximum likelihood was chosen. We performed bootstrapping using the autoMR convergence criteria.

Competence Machinery and Vaginolysin Identification

We systematically identified competence genes, tad pilus assembly homologs, and vaginolysin using PIRATE and Roary gene homolog output, as well as manually using the Prokka (Seemann, 2014) annotations of the Gardnerella isolates. We included homologs of cinA, recA, comEA, comEC, cpaB, cpaF, tadB, tadC, tadE, and tadG in our selection analyses (Tomich et al., 2007; Yeoman et al., 2010). We aligned each gene with PRANK and constructed individual maximum likelihood gene trees using FastTree (Price et al., 2010). We used the aBSREL method implemented in Hyphy (Smith et al., 2015) to test for selection along branches in the phylogenies of competence genes. To investigate the effects of LGT and selection on the diversity of the competence machinery, vaginolysin, and the core genome, we compared these genomic regions in the two species that are well-sampled: G. vaginalis and G. piotii. We calculated dN/dS using the yn00 implementation (Yang and Nielsen, 2000) in PAML (Yang, 1997).

Prophage Identification

We used ProphET (Reis-Cunha et al., 2019) to detect prophage in our collection of genomes. Briefly, ProphET performs a similarity search of annotated proteins from bacterial genomes against a database of known phage proteins to identify prophage within bacterial genomes. ProphET discards regions with a low density of phage-associated genes. Next, we blasted the known prophage regions against a custom nucleotide database of the de novo assembled contigs from all 106 Gardnerella spp. isolates to identify additional phage. To assess whether we missed prophage that were split between contigs, we blasted the known prophage against the custom database and filtered the results to identify hits found within 50 bp of a contig end and plotted the sequence length distribution of these hits. After identifying prophage regions, we calculated pairwise mash (Ondov et al., 2016) distances, which is based on shared k-mer (sequences of length k) content between prophage nucleotide sequences. Using these pairwise distances, we performed multidimensional scaling (MDS) to identify clusters of similar prophage. We also created a presence/absence matrix for each prophage cluster in our data set. To compare genetic content of each prophage cluster, we annotated the nucleotide sequences of all prophage using Prokka (Seemann, 2014), and cross referenced these results with the Clusters of Orthologous Groups (COG) database (Tatusov et al., 2000).

CRISPR/cas Identification

We identified CRISPR/cas genes in genome annotations produced by Prokka. Additionally, we used the PIRATE output to identify homologous CRISPR/cas genes. To look for an association between the presence of CRISPR/cas and particular prophage cluster we used Fisher's Exact Test (Fisher, 1934) with Bonferroni correction (Bonferroni, 1935). We also performed a Mann-Whitney-Wilcoxon test (Mann and Whitney, 1947) to test for an association between the presence of CRISPR/cas and genome assembly length.

Restriction Modification Identification

We used the Prokka annotations and PIRATE output to identify homologous clusters of genes associated with restriction modification (RM) systems. We created a presence/absence matrix for each RM associated gene in our data set.

Pangenome Diversity Analyses

We calculated pangenome accumulation and rarefaction curves of G. vaginalis and G. piotii isolates. Additionally, we calculated the gene frequency of the accessory genomes of both species. Using Egglib (De Mita and Siol, 2012), we calculated average π per accessory gene within and between G. vaginalis and G. piotii isolates using a subset of genes found at intermediate frequencies (1–99%) in both species. We performed a Kruskal-Wallis test in R to determine differences in average gene π values by group (Kruskal and Wallis, 1952). We then performed pairwise Mann-Whitney-Wilcoxon tests in R with Bonferroni correction to identify which distribution pairs were significantly different (Bonferroni, 1935; Mann and Whitney, 1947).

Results

Pangenome and Phylogenetic Analysis of Gardnerella spp.

Pan genome analysis of 106 clinical isolates of Gardnerella spp. identified 4,653 genes in the pan-genome. Six hundred and eight of these formed the “strict core” genome (i.e., found in 100% of isolates) whereas an additional 195 genes were found in 95% of the sample (Table 1).

For initial categorization of genes in this diverse sample, we used an amino acid identity threshold of 50% or greater, implemented in PIRATE (Bayliss et al., 2019), to identify 4,653 gene families. We found an average of 75% amino acid percent identity among genes found in at least 95% of the isolates. Based on this finding, we performed additional pangenome analysis using Roary with an amino acid percent identity threshold of 75% to identify core and accessory genes (Table 1). At this more restrictive threshold, we found 6,055 total gene families, 343 of which were found in at least 99% of isolates.

We found that G. vaginalis, G. piotii, G. swidsinskii, G. leopoldii, clade 3/D, and four unnamed additional clades were clearly differentiated in their core genomes, as shown by long branches separating them on a core genome phylogeny (Figure 1; Figure S1). The same pattern held in the accessory genome, where Gardnerella spp. could be clearly distinguished on the basis of accessory gene content (Figure 2).

Figure 1

Figure 2

Lateral Gene Transfer in the Gardnerella Core Genome

Our analyses of lateral gene transfer indicated that this species structure is maintained by barriers to recombination. We used FastGEAR to identify sub-groups within the sample and infer patterns of core genome recombination in the recent and remote past. FastGEAR analyses defined eight clades, consistent with previously published genetic classifications. Few recombination events were inferred between these eight clades (Figure 3A). Interestingly, between-clade recombination appeared to have occurred more commonly in the remote past (Figure 3B). More recombination events were inferred between closely related clades/species (i.e., between G. vaginalis and G. piotii) than distantly related clades/species (e.g., G. vaginalis and clade D/3).

Figure 3

Mechanisms of LGT

To identify potential mechanisms of LGT in Gardnerella, we first systematically searched for and identified a collection of competence related genes. A previous study identified four competence genes (cinA, recA, comEA, comEC) in a sample of three Gardnerella isolates (Yeoman et al., 2010). We expanded this finding to our collection of 106 Gardnerella isolates and identified 6 additional genes (cpaB, cpaF, tadB, tadC, tadE, tadG) involved in tad pilus assembly (Tomich et al., 2007). We found competence homologs to be encoded by most isolates (Figure 4). These genes were highly differentiated among species (e.g., Figure S2), mirroring patterns of diversity in the core genome. This suggests that competence genes are not commonly exchanged across species.

Figure 4

It is possible that patterns of diversity in competence genes reflect functional differentiation that could contribute to reproductive isolation among Gardnerella spp. To test this hypothesis, we performed selection analyses: in the event that divergence was driven by functional differentiation we expect to find evidence of positive selection reflected in excess non-synonymous (coding) variation separating species-specific versions of these genes. Pairwise values of dN/dS are consistent with purifying selection (i.e., selection against deleterious mutations) as the predominant influence on both core genes and competence homologs (Figure 5A). Estimates of branch specific omega values also indicate that the competence genes are generally evolving under purifying selection (Figure 5B; Figure S3). Thus, patterns of diversity in competence machinery do not suggest that their differentiation has been driven by selection for advantageous mutations/ functional differentiation.

Figure 5

Patterns of Diversity at a Virulence Locus

Vaginolysin is a pore-forming toxin and acute virulence factor (Gelber et al., 2008) previously identified as a core gene with evidence of between-clade recombination in a sample of 17 isolates (Ahmed et al., 2012). Here we sought to further investigate evolution of this locus: genes that mediate host pathogen interactions are often characterized by high levels of diversity reflecting selection for advantageous mutations (Andrews and Gojobori, 2004; Kennemann et al., 2011; Osório et al., 2013) and we were interested in dynamics of this well-characterized virulence factor in Gardnerella spp. By contrast with Ahmed et al., we did not find vaginolysin to be part of the core genome in this larger sample of isolates. Using a >80% amino acid similarity threshold in the PIRATE gene homolog output, we found vaginolysin to be present in 83% of Gardnerella isolates. A presence/absence matrix for this gene implies that vaginolysin has been gained or lost a limited number of times during evolution of our sample (Figure S4). The vaginolysin gene tree is not consistent with the core phylogeny, indicating that while it is infrequently lost or gained, it appears to be readily exchanged between Gardnerella spp. (Figures S4, S5). Interestingly, we found that vaginolysin appears to be evolving under strong constraint, similar to the core genome and competence machinery in G. vaginalis (Figure 5A). Taken together, these results suggest that the vaginolysin toxin performs an important function in Gardnerella spp., but bacteria occasionally adapt to loss of this function.

Prophage

To further explore mechanisms of LGT, we used ProphET (Reis-Cunha et al., 2019) to identify 130 prophage within our set of 106 genomes. We found Gardnerella spp. to encode between zero and four prophage per genome (median 1). Using mash distances (Ondov et al., 2016), we identified 5 clusters of prophage, with at least one prophage found in 70% of Gardnerella isolates. Prophage clusters 1 and 2 are found across Gardnerella spp. (Figure 6), while prophage clusters 3 and 4 are restricted to G. vaginalis and G. piotii, with the exception of one cluster 4 prophage found in a G. swidsinskii isolate. Cluster 5 is the smallest group, found in 2 isolates from an unnamed Gardnerella spp. clade (Figure 6). Phages in clusters 2 and 3 form well-differentiated sub-clusters that correspond with host species designations (Figures 7B,C), whereas cluster 1 and 4 phages from different hosts co-mingle (Figures 7A,D). This suggests there are barriers to between-species transfer of phage clusters 2 and 3 whereas clusters 1 and 4 are readily transferred among diverse bacterial hosts. We found the genetic content of prophage to vary between clusters (Figure S6); however, the majority of the genes were uncharacterized hypothetical genes (82%). To assess if we missed prophage split over contigs, we blasted all prophage against a custom database of all de novo assembled contigs and filtered the results to find hits within 50 bp of the end of a contig. We plotted the sequence length distribution of these blast results and found the majority were very small in length, which suggests they are unlikely to be identified prophage (Figure S7).

Figure 6

Figure 7

CRISPR/cas

CRISPR/cas protect bacterial genomes from foreign DNA, including phage, and thus could potentially play a role in maintaining genetic barriers among Gardnerella spp. Presence/absence matrix of CRISPR/cas loci in our sample is shown in Figure 6. Patterns of carriage were variable, with some strains carrying up to 10 distinct loci and others without any loci identified. Certain loci tended to co-occur. Loss of CRISPR/cas can lead to proliferation of mobile genetic elements (Hullahalli et al., 2018), in which case we might expect genomes lacking CRISPR/cas to accumulate mobile genetic elements and increase in size. To test this hypothesis, we compared the presence of CRISPR/cas to the total length of de novo assembled contigs. We found that CRISPR/cas is not associated with increased genome length (Mann-Whitney-Wilcoxon test, W = 1,058, p = 0.08) (Figure S8).

Restriction Modification

Restriction modification (RM) systems can also protect bacterial genomes by cleaving foreign DNA (Tock and Dryden, 2005), and thus may shape barriers to lateral gene transfer in Gardnerella spp. To identify genes involved in restriction modification (RM) systems, we used PIRATE homologous gene family output and Prokka genome annotations. We found a wide diversity of RM genes found at varying frequencies across the Gardnerella spp. phylogeny, and Gardnerella spp. are not defined by the presence/absence of particular RM gene families (Figure S9).

Codon Usage in G. vaginalis and G. piotii

Differences in codon usage (Tuller et al., 2011) are another potential mechanism driving differentiation and genetic isolation of Gardnerella species. To investigate this hypothesis, we focused on G. vaginalis and G. piotii, sister species in the Gardnerella phylogeny that are common and well-sampled (50 G. vaginalis, 28 G. piotii). We found codon usage to be similar in G. vaginalis and G. piotii isolates and thus it does not appear to be an explanation for the reproductive isolation of these species (t-test, p > 0.99) (Figure S10).

Lateral Gene Transfer in Core Genomes of G. vaginalis and G. piotii

To further characterize patterns of recombination in G. vaginalis and G. piotii, we identified recombinant tracts in core genomes with Gubbins. LGT appears qualitatively more frequent in G. vaginalis than in G. piotii (Figure 8). The proportion of the total alignment affected by recombination is similar in the two species (0.94 for G. vaginalis and 0.92 for G. piotii), in an uncorrected analysis. To account for differences in sample size in G. vaginalis and G. piotii isolates, we sub-sampled the number of G. vaginalis isolates to the number of G. piotii isolates and used Gubbins to identify recombination in the subsampled dataset. The mean proportion of each isolate's core genome affected by recombination is greater for G. vaginalis (26.1%) than for G. piotii (23.1%) (Mann-Whitney-Wilcoxon test, W = 473, p = 0.0138) (Figure 9). The G. vaginalis network contains more reticulations than G. piotii, which also supports overall higher rates of recombination in G. vaginalis (Figure S11).

Figure 8

Figure 9

Patterns of Lateral Gene Transfer in Accessory Genomes of G. vaginalis and G. piotii

Inter-species barriers to gene exchange were evident in comparisons of the accessory genomes of G. vaginalis and G. piotii. Accessory gene frequencies in G. vaginalis and G. piotii isolates indicated that most accessory genes were unique to an individual species, with few genes found at similar frequencies across species (Figure 10). If accessory genes were freely exchanged between species, we would expect their frequencies to equilibrate.

Figure 10

To further characterize lateral gene transfer in the accessory genome, we quantified pairwise diversity (π) for accessory genes common to G. vaginalis and G. piotii. We performed a Kruskal-Wallis test and found there was a statistically significant difference between average gene π values by group (H = 214.7, p < 2.2e-16). We then performed pairwise Mann-Whitney-Wilcoxon tests with Bonferroni correction and found the distributions of average gene π values of G. vaginalis (W = 373,928, p < 3.32e-13) and G. piotii isolates (W = 399,762, p < 3.37e-13) were lower than between species (Figure 11). In species that regularly exchange accessory gene content, we expect to see similar levels of diversity in between- and within-species comparisons. However, we observed lower diversity within species, indicating that accessory gene variants are transferred more frequently within species than between them.

Figure 11

Differences in the amount of recombination between G. vaginalis and G. piotii appear to affect not only the core genome, but also the accessory genome. G. vaginalis has a larger pangenome than G. piotii, consistent with acquisition of novel gene content from diverse sources (Figure 12).

Figure 12

Discussion

Structure of Gardnerella spp. Pangenome

Gardnerella spp. comprise diverse bacteria with a distinct genetic structure (Ahmed et al., 2012). This diversity has been delineated both phylogenetically and metabolically (Ahmed et al., 2012; Cornejo et al., 2018; Vaneechoutte et al., 2019). Using WGS data from 106 Gardnerella isolates, we identified eight major clades from a core genome alignment that are consistent with published genomic analyses (Ahmed et al., 2012; Devault et al., 2017; Schellenberg et al., 2017; Cornejo et al., 2018; Vaneechoutte et al., 2019). Our data primarily consist of isolates from G. vaginalis and G. piotii, which are found most commonly in clinical samples with higher Nugent scores (Janulaitiene et al., 2017).

The eight major clades in our sample were differentiated with respect to both allelic variation in their core genomes (Figure 1) and gene content in their accessory genomes (Figure 2). Prior research found evidence of frequent homologous recombination within clades/species (Ahmed et al., 2012) and functional differentiation of the accessory genome (Cornejo et al., 2018). Our results extend these observations with the observation that allelic variants in the core genome appear to be more readily exchanged within species than between them (Figures 3, 8). In addition, we found evidence of barriers to recombination in the accessory genome, as shared accessory genes are more similar within species than between species (Figure 11). Accessory genes are not maintained at similar frequencies across species, suggesting that selection pressures for the shared accessory genes are not the same across species (Figure 10).

Mechanisms and Barriers to LGT in Gardnerella

Previous studies identified competence genes in a handful of Gardnerella isolates (Yeoman et al., 2010). We systematically examined our sample for previously identified competence genes as well as genes known to encode tad pili (Tomich et al., 2007; Yeoman et al., 2010). We found these genes to be encoded by most isolates from the eight clades/species in our sample (Figure 4). The competence machinery, was however, highly genetically differentiated among clades/species (Figure S2). This raises the possibility of functional differentiation among competence genes, which could contribute to genetic isolation of clades/species (Porse et al., 2018). We investigated this hypothesis by computing the ratio of non-synonymous to synonymous variation in pairwise comparisons and across phylogenies of competence genes (Figure 5; Figure S3). We did not find any evidence of functional divergence among competence genes. Our results instead suggest that competence genes are under similar evolutionary pressures as the core genome, with both evolving under purifying selection (Figure 5; Figure S3).

A previous study identified predicted prophage genes in Gardnerella spp. (Malki et al., 2016). This does not necessarily indicate the presence of prophage, particularly if phage orthologs are not found within a cluster of phage associated genes. We used ProphET (Reis-Cunha et al., 2019) to identify prophage clusters, which we found in 70% of isolates (Figure 6). One limitation of this method is that prophage found in poorly assembled regions or across multiple contigs in the assembly may remain unidentified thus the true carriage frequency of prophage could be higher than our estimate. However, we blasted the prophage regions against the de novo assembled contigs and found that hits along the end of contigs are very small in length, suggesting they are not unidentified prophage (Figure S7).

CRISPR/cas are adaptive immune systems that can protect bacterial genomes from mobile genetic elements, such as phage, and thus could potentially shape and reinforce the genetic barriers in Gardnerella spp. Loss of CRISPR has been shown to enable the proliferation of mobile elements in Enterococcus (Pleckaityte et al., 2012; Hullahalli et al., 2018). We did not find an association between the presence of CRISPR/cas genes and the absence of any prophage clusters, suggesting these interactions may not be straightforward (Figure 6). Restriction modification (RM) is another mechanism by which foreign DNA is cleaved and thus prevented from integrating into bacterial genomes (Tock and Dryden, 2005). We found a wide diversity of RM genes at varying frequencies among Gardnerella spp. (Figure S9). The role they play in shaping and maintaining LGT barriers is unknown.

Differences in codon usage are another possible barrier to recombination among species (Tuller et al., 2011). To investigate this possibility, we compared codon usage between the two most well-sampled species. We did not, however, find evidence of differential codon usage among G. vaginalis and G. piotii isolates (Figure S10). It's intriguing that codon usage appears to be harmonized amongst these highly genetically differentiated species. By comparison, a study of a single relatively clonal species (S. aureus) found evidence of differences in codon usage among ecotypes (Richardson et al., 2018). We and others have found evidence that the pore-forming toxin vaginolysin is freely exchanged among Gardnerella spp. (Figure S5) (Ahmed et al., 2012). Harmonization of codon usage could facilitate exchange of genes like vaginolysin that are critical to diverse Gardnerella spp. We found evidence that vaginolysin genes are evolving under strong constraint (Figure 5A), which supports the idea that the toxin is important to the fitness of diverse Gardnerella spp.

Our results indicate that clades/species of Gardnerella spp. are reproductively isolated despite being found in the same niche. This differentiation does not appear to be driven by functional differentiation of competence genes, nor by differences in patterns of codon usage. The likelihood of LGT events and compatibility of transferred genes with the recipient genome increases between closely related genomes (Popa and Dagan, 2011). It is probable that general patterns of differentiation in the core and accessory genomes of Gardnerella spp. have a role maintaining species separation (Porse et al., 2018). It is intriguing that between-species LGT appears to have been more common in the remote past (Figure 3B), when the species may have shared more sequence similarity. Other potential mediators of reproductive isolation among species may be restriction-modification systems and CRISPR/cas systems, both of which can target degradation of foreign DNA, and thus shape patterns of LGT (Tock and Dryden, 2005; Marraffini and Sontheimer, 2008; Dupuis et al., 2013; Hullahalli et al., 2018). We did not find evidence of interactions between CRISPR/cas loci and prophage (Figure 6); suggesting interactions among these elements are complex in Gardnerella spp.

Patterns of Within and Between Species LGT in G. vaginalis and G. piotii

We found patterns of LGT to vary among Gardnerella spp. Specifically, we found evidence of distinct patterns of within-species recombination in comparisons of G. vaginalis and G. piotii. G. vaginalis appears to engage more frequently in LGT (Figure 8; Figure S11) as a larger proportion of each isolate's core genome is estimated to have been affected by recombination (Figure 9). Additionally, G. vaginalis has a larger pangenome (Figure 12), consistent with higher levels of gene importation in the accessory genome. These findings raise the possibility that G. vaginalis exchanges DNA with more diverse partners than does G. piotii, potentially with other non-Gardnerella members of the complex polymicrobial BV biofilm.

We do not know the underlying mechanism for differences in LGT between G. vaginalis and G. piotii. Variation in rates of transduction as well as restriction-modification systems could play a role. A previous study found G. vaginalis to be enriched for “phage-associated protein” and three genes involved in a type I RM system (Cornejo et al., 2018). We observed a wide diversity of RM systems, with multiple type I variants found within each species (Figure S9). Thus, even if RM systems explain the differences in LGT patterns within species, it is not a one to one relationship between species and RM systems. Teasing apart those interactions will require more research.

In conclusion, Gardnerella spp. are genetically distinct in both their core and accessory genomes. We found evidence of more within species LGT in the core and accessory genomes, suggesting active maintenance of reproductive barriers despite similar patterns of codon usage. The putative competence machinery is genetically differentiated between clades/species; however, we found no evidence of functional divergence/positive selection driving clade/species separation. We identified a larger pangenome in G. vaginalis than in G. piotii as well as more LGT in the core genome, suggesting G. vaginalis engages more frequently in LGT with more diverse partners. Taken together, our results demonstrate that co-localized bacterial populations can maintain a complex genetic structure in which genetic exchange appears to be restricted to specific sub-populations with exceptions for individual genes (e.g., vaginolysin). The forces maintaining this structure are yet to be fully elucidated but likely include patterns of sequence similarity and possibly phage, CRISPR/cas, RM systems and interactions among them. Defining evolutionary interactions in bacterial populations helps to illuminate how clinically important traits such as antibiotic resistance and virulence emerge and are maintained in these complex communities.

References

• 1

  AhmedA.EarlJ.RetchlessA.HillierS. L.RabeL. K.CherpesT. L.et al. (2012). Comparative genomic analyses of 17 clinical isolates of Gardnerella vaginalis provide evidence of multiple genetically isolated clades consistent with subspeciation into genovars. J. Bacteriol. 194, 3922–3937. 10.1128/JB.00056-12

  • CrossRef
  • Google Scholar

• 2

  AndrewsS. (2010). FastQC: A Quality Control tool for High Throughput Sequence Data. Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/

  • Google Scholar

• 3

  AndrewsT. D.GojoboriT. (2004). Strong positive selection and recombination drive the antigenic variation of the PilE protein of the human pathogen neisseria meningitidis. Genetics166, 25–32. 10.1534/genetics.166.1.25

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 4

  BalashovS. V.MordechaiE.AdelsonM. E.GygaxS. E. (2014). Identification, quantification and subtyping of Gardnerella vaginalis in noncultured clinical vaginal samples by quantitative PCR. J. Med. Microbiol. 63, 162–175. 10.1099/jmm.0.066407-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 5

  BankevichA.NurkS.AntipovD.GurevichA. A.DvorkinM.KulikovA. S.et al. (2012). SPAdes: a new genome assembly algorithm and its applications to single-cell sequencing. J. Comput. Biol. 19, 455–477. 10.1089/cmb.2012.0021

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 6

  BaylissS. C.ThorpeH. A.CoyleN. M.SheppardS. K.FeilE. J. (2019). PIRATE: a fast and scalable pangenomics toolbox for clustering diverged orthologues in bacteria. GigaScience8:598391. 10.1093/gigascience/giz119

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 7

  BaymM.KryazhimskiyS.LiebermanT. D.ChungH.DesaiM. M.KishonyR. (2015). Inexpensive multiplexed library preparation for megabase-sized genomes. PLoS ONE10:e0128036. 10.1371/journal.pone.0128036

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 8

  BonferroniC. E. (1935). Il calcolo delle assicurazioni su gruppi di teste. Studi Onore Profr. Salvatore Ortu Carboni13–60.

  • Google Scholar

• 9

  ChikhiR.MedvedevP. (2014). Informed and automated k-mer size selection for genome assembly. Bioinformatics30, 31–37. 10.1093/bioinformatics/btt310

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 10

  ClarkS. C.EganR.FrazierP. I.WangZ. (2013). ALE: a generic assembly likelihood evaluation framework for assessing the accuracy of genome and metagenome assemblies. Bioinformatics29, 435–443. 10.1093/bioinformatics/bts723

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 11

  CornejoO. E.HickeyR. J.SuzukiH.ForneyL. J. (2018). Focusing the diversity of Gardnerella vaginalis through the lens of ecotypes. Evol. Appl. 11, 312–324. 10.1111/eva.12555

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 12

  CroucherN. J.PageA. J.ConnorT. R.DelaneyA. J.KeaneJ. A.BentleyS. D.et al. (2015). Rapid phylogenetic analysis of large samples of recombinant bacterial whole genome sequences using gubbins. Nucleic Acids Res. 43:e15. 10.1093/nar/gku1196

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 13

  De MitaS.SiolM. (2012). EggLib: processing, analysis and simulation tools for population genetics and genomics. BMC Genet. 13:27. 10.1186/1471-2156-13-27

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 14

  DevaultA. M.MortimerT. D.KitchenA.KiesewetterH.EnkJ. M.GoldingG. B.et al. (2017). A molecular portrait of maternal sepsis from byzantine troy. ELife6:e20983. 10.7554/eLife.20983

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 15

  DupuisM.-È.VillionM.MagadánA. H.MoineauS. (2013). CRISPR-Cas and restriction–modification systems are compatible and increase phage resistance. Nat. Commun. 4:2057. 10.1038/ncomms3087

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 16

  FisherR. A. (1934). Statistical Methods for Research Workers, 5th Edn.

  • Google Scholar

• 17

  GarrisonE.MarthG. (2012). Haplotype-based variant detection from short-read sequencing. arXiv preprint arXiv:1207.3907.

  • Google Scholar

• 18

  GelberS. E.AguilarJ. L.LewisK. L. T.RatnerA. J. (2008). Functional and phylogenetic characterization of Vaginolysin, the human-specific cytolysin from Gardnerella vaginalis. J. Bacteriol. 190, 3896–3903. 10.1128/JB.01965-07

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 19

  GhodsiM.HillC. M.AstrovskayaI.LinH.SommerD. D.KorenS.et al. (2013). De novo likelihood-based measures for comparing genome assemblies. BMC Res. Notes6:334. 10.1186/1756-0500-6-334

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 20

  GurevichA.SavelievV.VyahhiN.TeslerG. (2013). QUAST: quality assessment tool for genome assemblies. Bioinforma. Oxf. Engl. 29, 1072–1075. 10.1093/bioinformatics/btt086

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 21

  HarwichM. D.AlvesJ. M.BuckG. A.StraussJ. F.PattersonJ. L.OkiA. T.et al. (2010). Drawing the line between commensal and pathogenic Gardnerella vaginalis through genome analysis and virulence studies. BMC Genomics11:375. 10.1186/1471-2164-11-375

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 22

  HilbertD. W.SchuylerJ. A.AdelsonM. E.MordechaiE.SobelJ. D.GygaxS. E. (2017). Gardnerella vaginalis population dynamics in bacterial vaginosis. Eur. J. Clin. Microbiol. Infect. Dis.36, 1269–1278. 10.1007/s10096-017-2933-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 23

  HillJ. E.AlbertA. Y. K.VOGUE Research Group. (2019). Resolution and cooccurrence patterns of gardnerella leopoldii, G. swidsinskii, G. piotii, and G. vaginalis within the vaginal microbiome. Infect. Immun. 87:e00532–19. 10.1128/IAI.00532-19

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 24

  HullahalliK.RodriguesM.NguyenU. T.PalmerK. (2018). An attenuated CRISPR-Cas system in enterococcus faecalis permits DNA acquisition. MBio9:e00414–18. 10.1128/mBio.00414-18

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 25

  HuntM.KikuchiT.SandersM.NewboldC.BerrimanM.OttoT. D. (2013). REAPR: a universal tool for genome assembly evaluation. Genome Biol. 14:R47. 10.1186/gb-2013-14-5-r47

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 26

  HusonD. H.BryantD. (2006). Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol.23, 254–267.

  • Pubmed Abstract
  • Google Scholar

• 27

  JanulaitieneM.PaliulyteV.GrincevicieneS.ZakarevicieneJ.VladisauskieneA.MarcinkuteA.et al. (2017). Prevalence and distribution of Gardnerella vaginalis subgroups in women with and without bacterial vaginosis. BMC Infect. Dis. 17:394. 10.1186/s12879-017-2501-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 28

  KatohK.StandleyD. M. (2014). MAFFT: iterative refinement and additional methods. Methods Mol. Biol. Clifton N.J.1079, 131–146. 10.1007/978-1-62703-646-7_8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 29

  KennemannL.DidelotX.AebischerT.KuhnS.DrescherB.DroegeM.et al. (2011). Helicobacter pylori genome evolution during human infection. Proc. Natl. Acad. Sci. U.S.A. 108, 5033–5038. 10.1073/pnas.1018444108

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 30

  KorenS.TreangenT. J.HillC. M.PopM.PhillippyA. M. (2014). Automated ensemble assembly and validation of microbial genomes. BMC Bioinformatics15:126. 10.1186/1471-2105-15-126

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 31

  KruskalW. H.WallisW. A. (1952). Use of ranks in one-criterion variance analysis. J. Am. Stat. Assoc. 47, 583–621. 10.1080/01621459.1952.10483441

  • CrossRef
  • Google Scholar

• 32

  MalkiK.ShapiroJ. W.PriceT. K.HiltE. E.Thomas-WhiteK.SircarT.et al. (2016). Genomes of gardnerella strains reveal an abundance of prophages within the bladder microbiome. PLoS ONE11:e0166757. 10.1371/journal.pone.0166757

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 33

  MannH. B.WhitneyD. R. (1947). On a test of whether one of two random variables is stochastically larger than the other. Ann. Math. Stat. 18, 50–60. 10.1214/aoms/1177730491

  • CrossRef
  • Google Scholar

• 34

  MarraffiniL. A.SontheimerE. J. (2008). CRISPR interference limits horizontal gene transfer in staphylococci by targeting DNA. Science322, 1843–1845. 10.1126/science.1165771

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 35

  MostowyR.CroucherN. J.AndamC. P.CoranderJ.HanageW. P.MarttinenP. (2017). Efficient inference of recent and ancestral recombination within bacterial populations. Mol. Biol. Evol.34, 1167–82.

  • Pubmed Abstract
  • Google Scholar

• 36

  OndovB. D.TreangenT. J.MelstedP.MalloneeA. B.BergmanN. H.KorenS.et al. (2016). Mash: fast genome and metagenome distance estimation using MinHash. Genome Biol. 17:132. 10.1186/s13059-016-0997-x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 37

  OsórioN. S.RodriguesF.GagneuxS.PedrosaJ.Pinto-CarbóM.CastroA. G.et al. (2013). Evidence for diversifying selection in a set of mycobacterium tuberculosis genes in response to antibiotic- and nonantibiotic-related pressure. Mol. Biol. Evol. 30, 1326–1336. 10.1093/molbev/mst038

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 38

  PageA. J.CumminsC. A.HuntM.WongV. K.ReuterS.HoldenM. T. G.et al. (2015). Roary: rapid large-scale prokaryote pan genome analysis. Bioinformatics31, 3691–3693. 10.1093/bioinformatics/btv421

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 39

  PiotP.DyckE. V.PeetersM.HaleJ.TottenP. A.HolmesK. K. (1984). Biotypes of Gardnerella vaginalis. J. Clin. Microbiol. 20, 677–679. 10.1128/JCM.20.4.677-679.1984

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 40

  PleckaityteM.ZilnyteM.ZvirblieneA. (2012). Insights into the CRISPR/Cas system of Gardnerella vaginalis. BMC Microbiol. 12:301. 10.1186/1471-2180-12-301

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 41

  PopaO.DaganT. (2011). Trends and barriers to lateral gene transfer in prokaryotes. Curr. Opin. Microbiol. 14, 615–623. 10.1016/j.mib.2011.07.027

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 42

  PorseA.SchouT. S.MunckC.EllabaanM. M. H.SommerM. O. A. (2018). Biochemical mechanisms determine the functional compatibility of heterologous genes. Nat. Commun. 9:522. 10.1038/s41467-018-02944-3

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 43

  PriceM. N.DehalP. S.ArkinA. P. (2010). FastTree 2 – approximately maximum-likelihood trees for large alignments. PLoS ONE5:e9490. 10.1371/journal.pone.0009490

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 44

  RahmanA.PachterL. (2013). CGAL: computing genome assembly likelihoods. Genome Biol. 14:R8. 10.1186/gb-2013-14-1-r8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 45

  Reis-CunhaJ. L.BartholomeuD. C.MansonA. L.EarlA. M.CerqueiraG. C. (2019). ProphET, prophage estimation tool: a stand-alone prophage sequence prediction tool with self-updating reference database. PLoS ONE. 14:e0223364. 10.1371/journal.pone.0223364

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 46

  RiceP. (2000). EMBOSS: the european molecular biology open software suite. Trends Genet. 16, 276–277. 10.1016/S0168-9525(00)02024-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 47

  RichardsonE. J.BacigalupeR.HarrisonE. M.WeinertL. A.LycettS.VrielingM.et al. (2018). Gene exchange drives the ecological success of a multi-host bacterial pathogen. Nat. Ecol. Evol. 2, 1468–1478. 10.1038/s41559-018-0617-0

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 48

  SchellenbergJ. J.PattersonM. H.HillJ. E. (2017). Gardnerella vaginalis diversity and ecology in relation to vaginal symptoms. Res. Microbiol. 168, 837–844. 10.1016/j.resmic.2017.02.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 49

  SeemannT. (2014). Prokka: rapid prokaryotic genome annotation. Bioinformatics30, 2068–2069. 10.1093/bioinformatics/btu153

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 50

  SmithM. D.WertheimJ. O.WeaverS.MurrellB.SchefflerK.Kosakovsky PondS. L. (2015). Less is more: an adaptive branch-site random effects model for efficient detection of episodic diversifying selection. Mol. Biol. Evol. 32, 1342–1353. 10.1093/molbev/msv022

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 51

  StamatakisA. (2014). RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics30, 1312–1313. 10.1093/bioinformatics/btu033

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 52

  TatusovR. L.GalperinM. Y.NataleD. A.KooninE. V. (2000). The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res. 28, 33–36. 10.1093/nar/28.1.33

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 53

  TockM. R.DrydenD. T. (2005). The biology of restriction and anti-restriction. Curr. Opin. Microbiol. 8, 466–472. 10.1016/j.mib.2005.06.003

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 54

  TomichM.PlanetP. J.FigurskiD. H. (2007). The tad locus: postcards from the widespread colonization island. Nat. Rev. Microbiol. 5, 363–375. 10.1038/nrmicro1636

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 55

  TullerT.GirshovichY.SellaY.KreimerA.FreilichS.KupiecM.et al. (2011). Association between translation efficiency and horizontal gene transfer within microbial communities. Nucleic Acids Res. 39, 4743–4755. 10.1093/nar/gkr054

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 56

  VaneechoutteM.GuschinA.VanL. S.GansemansY.VanF. N.CoolsP. (2019). Emended description of Gardnerella vaginalis and description of Gardnerella leopoldii sp. nov., Gardnerella piotii sp. nov. and Gardnerella swidsinskii sp. nov., with delineation of 13 genomic species within the genus Gardnerella. Int. J. Syst. Evol. Microbiol. 69, 679–687. 10.1099/ijsem.0.003200

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 57

  WoodD. E.SalzbergS. L. (2014). Kraken: ultrafast metagenomic sequence classification using exact alignments. Genome Biol. 15:R46. 10.1186/gb-2014-15-3-r46

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  YangZ. (1997). PAML: a program package for phylogenetic analysis by maximum likelihood. Comput. Appl. Biosci.13, 555–556. 10.1093/bioinformatics/13.5.555

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 59

  YangZ.NielsenR. (2000). Estimating synonymous and nonsynonymous substitution rates under realistic evolutionary models. Mol. Biol. Evol. 17, 32–43. 10.1093/oxfordjournals.molbev.a026236

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 60

  YeomanC. J.YildirimS.ThomasS. M.DurkinA. S.TorralbaM.SuttonG.et al. (2010). Comparative genomics of Gardnerella vaginalis strains reveals substantial differences in metabolic and virulence potential. PLoS ONE5:e12411. 10.1371/journal.pone.0012411

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 61

  ZerbinoD. R.BirneyE. (2008). Velvet: algorithms for de novo short read assembly using de bruijn graphs. Genome Res. 18, 821–829. 10.1101/gr.074492.107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 62

  ZiminA.MarçaisG.PuiuD.RobertsM.SalzbergS. L.YorkeJ. A. (2013). The MaSuRCA genome assembler. Bioinformatics29, 2669–77. 10.1093/bioinformatics/btt476

  • Pubmed Abstract
  • CrossRef
  • Google Scholar