Genome sequencing and analysis of the first complete genome of Lactobacillus kunkeei strain MP2, an Apis mellifera gut isolate

DNA isolation and genome sequencing

The L. kunkeei strain MP2 reported in this study was isolated in a previous study from the gut of a honey bee (Apis mellifera), collected from commune hives in the Maria Pinto area, Melipilla Province in the Central zone of Chile (Olmos et al., 2014). Collection and use of honey bees for these studies was reviewed and approved by the Bioethics Committee of Universidad Mayor, which is governed by the regulations of the Animal Health Services of Chile.

For DNA extraction, L. kunkeei MP2 colonies were cultured in MRS broth (37 °C, 5% CO₂) and genomic DNA was obtained using a silica-based protocol (Boom et al., 1999). Briefly, bacterial pellet was lysed using a solution composed of SDS 10%, proteinase K (10 mg/mL, Thermo Scientific) and lysozyme (5 mg/mL, Pierce) at 37 °C for 60 min. The lysate was mixed with guanidine chloride 6M and a silica suspension (50% w/v) and incubated for 10 min. The silica was centrifuged, and DNA was recovered, after washes with 70% alcohol, into sterile, free nuclease water. Approximately 13.5 µg of DNA were used to construct sequencing libraries with an average insert size of 20 kb, and sequenced using one SMRT cell (P6-C4 Chemistry) on a PacBio RSII sequencer (Pacific Biosciences) at the UCSD IGM Genomics Center.

Genome assembly and annotation

Raw reads (∼1 Gbps) were processed to remove SMRT bell adapters, short and low-quality reads (<80% accuracy) using SMRT Analysis version 2.3. A total of 154,044 filtered reads (average length, 9 Kb) were used for de novo assembly using Celera Assembler version 8.3 (Myers et al., 2000), with self-correction of the PacBio reads (Berlin et al., 2015). Polishing was done using Quiver, using SMRT Analysis version 2.3. Comparisons between the previously sequenced L. kunkeei MP2 genome (Accession number PRJNA257367) (Olmos et al., 2014), as well with the other two available genome sequences (Porcellato et al., 2015; Djukic et al., 2015) were performed using MUMMER (Kurtz et al., 2004). Genome annotation was performed using Prokka version 1.11 (Seemann, 2014). The predicted CDS were classified into EggNOG categories using HMMER version 3.1 (http://hmmer.org) against the EggNOG 4.1 database (Powell et al., 2014) with an E-value cutoff of 1E-05. Genomic islands were annotated using IslandViewer 3 (Dhillon et al., 2015), and possible phage sequences were searched using PHAST (Zhou et al., 2011). Genome visualization was done using Circos version 0.69 (Krzywinski et al., 2009). The genome sequence an assembly is available at NCBI with the accession number PRJNA298292.

Pan-genome analysis

Comparative genomic analysis was performed from a list of selected genomes from Lactobacillus species, as well as other strains of L. kunkeei (Table 1). To avoid possible biases in the comparisons due to different annotation procedures, all of the genomes were re-annotated using Prokka version 1.11 (Seemann, 2014). Comparisons between all the selected genomes, were done using Roary (Page et al., 2015), with a blast identity cutoff of 97% for the comparison between L. kunkeei strains, and a 40% for the comparison between Lactobacillus species. In addition at the genus level, paralog splitting was disabled. Visualization of the pan-genome data was done using Anvi’o (Eren et al., 2015).

Phylogenetic analysis

16S rRNA gene sequences were obtained from the Silva database (Quast et al., 2013), and aligned using the SINA webserver aligner (Pruesse, Peplies & Glöckner, 2012) with default parameters, and using the Bacteria variability profile. For the L. kunkeei strain tree, we used the core genome (genes shared by all L. kunkeei strains), aligning the genes using Prank (Löytynoja, 2014). For both cases, the phylogenetic trees were generated using FastTree v 2.1.7 (Price, Dehal & Arkin, 2010), with the–slow option.

Assembly description

The PacBio reads obtained for L. kunkeei MP2 were assembled using MHAP (Berlin et al., 2015) implemented in the PBcR pipeline (Celera Assembler 8.3) (Myers et al., 2000). This de novo assembly resulted in one contig, representing the complete genome of L. kunkeei MP2 in a single 1,614,522 nt chromosome. A total of 1,468 CDS were predicted in the genome, 67 tRNA and 5 copies of the ribosomal operon. Functional annotation was done using EggNOG V 4.1 (Powell et al., 2014), and the summary of functional categories is shown on Table 2 (gene annotation on Table S1). The %GC content of the genome was 36.9%, and several features of interests, such as the presence of prophage regions, were found. With this assembly, we were able to differentiate the three ribosomal operons that are present in the chromosome, something that was not possible in the previous sequenced genome of this strain (Olmos et al., 2014).

A comparison of the assembly of L. kunkeei MP2 obtained in this work, with the previously obtained using Illumina sequencing (Olmos et al., 2014), is shown in Fig. 1. All of the previous assembled contigs mapped to the current assembly, and several gaps on the sequence were completed in this new version of the genome.

Central Metabolism of L. kunkeei MP2

Prophage insertions

Previous work in other Lactobacillus species, reported the presence of regions with prophage genes in their genomes, including species such as L. rhamnosus, L. gasseri, L. salivarius, L. casei, L. lactis, and L. johnsonii (Ventura et al., 2004; Ventura et al., 2006; Kankainen et al., 2009; Savabi et al., 2014; Baugher, Durmaz & Klaenhammer, 2014). This shows the widespread abundance of prophages in the genomes of Lactobacillus species, a characteristic shared by L. kunkeei MP2. Two regions were identified by PHAST (Zhou et al., 2011), as putative prophage insertions (Table S2). One of them, located in the region between 594,506 and 613,136, was found to be present in all the 23 Lactobacillus genomes used in this work. The second region, located around 32,973–75,092, was found to be unique to L. kunkeei MP2, compared to other strains of L. kunkeei, as well as other Lactobacillus species. In at least one Lactobacillus species (L. gasseri), the presence of these inserted phages has been associated with the horizontal transfer of genes (Baugher, Durmaz & Klaenhammer, 2014), suggesting a possible role for these elements within the genome of L. kunkeei MP2. However, the detailed mechanisms, as well as the possible adaptive consequences of such events, need to be explored in more detail in the future.

Comparison of L. kunkeei MP2 with other L. kunkeei strains

We performed a comparative genomic analysis of MP2 against sixteen publicly available genome sequences of L. kunkeei strains (Table 1). This analysis can provide us with a snapshot of the unique features that are present in this strain, based on its gene content. Recent work, highlighted the important role that genome reduction played in the evolution of L. kunkeei (Tamarit et al., 2015), which suggests that we would expect to find only a few genes truly unique to L. kunkeei MP2, compared to other strains.

Pangenome analysis resulted in the identification of a set of 853 core genes (present in all the strains). For the shell genes, we found that 813 genes are found between 15% to 95% of the analyzed strains, while 1661 genes are present in less than 15% of the strains. Focusing on the MP2 strain, based on this analysis we identified 113 genes that are not present in any of the other strains, representing close to 7.7% of the total number of genes.

Using the pangenome matrix, we can visualize the relationships between the genomes based only on the clustering of the proteins (presence or absence of a gene in a protein group). Visualizations of the results (Fig. 2), show that the strains LAni, LAce, LAan and EFB6, cluster in the same group as MP2, based on their gene content. This correlates with the phylogenetic placement of MP2 within the same group, based on a concatenated alignment of the core genome (853 genes; 807,585 nucleotides) (Fig. 3). The closest strain on the tree is EFB6, which was also isolated from the gut of A. mellifera (although from a larvae affected by European foulbrood) (Djukic et al., 2015). Currently there is no evidence of association between particular strains and a particular species of Apis (Tamarit et al., 2015). A better understanding of the interactions between the gut microbiota and the host needs to incorporate additional variables, such as the surrounding environment (including flowers and fruit), seasonal variations, among other possible elements. These factors could play a strong role in the association between an Apis species and a L. kunkeei strain.

One of the main differences of MP2 compared to the other L. kunkeei genomes, is the presence of multiple phage genes inserted in several parts of the genome. One of these unique phage regions can be found at coordinates 31,034–75,092 (Fig. 1). It is a large region, which includes several phage-related proteins, including structural and replications components. Sequence analysis using Blast, shows that these proteins are related to phages that infect Gram-positive Bacteria, such as Bacillus (Hastings et al., 2013), Listeria (Dorscht et al., 2009), Enterococcus (Yasmin et al., 2010), and Staphylococcus (Chang et al., 2013) (Table S3).

Comparison of L. kunkeei MP2 with other Lactobacillus strains

Phylogenetic reconstructions using the sequence that encodes for the 16S rRNA gene, shows that the closes species to L. kunkeei MP2 is L. sanfranciscensis (Fig. 4). This placement is in agreement with previous diversity analysis performed on acidophilic bacteria (McFrederick et al., 2012), relating the taxonomy of both species.

Whole genome comparisons between L. kunkeei MP2 and other Lactobacillus species, discovered several unique genes. One example is gtfC, which encodes for a glucosyltransferase, which has been extensively studied in Streptococcus mutans, where is expressed in the presence of carbohydrates such as sucrose, D-glucose, D-fructose, among others (Shemesh et al., 2006). GtfC (as well as GftB), is also considered a virulence factor in S. mutants, promoting bacterial adhesion to smooth surfaces and cells (Tsumori & Kuramitsu, 1997). Also, GtfC is part of the synthesis route of a mixture of insoluble and soluble glucans, which are important components of cariogenic biofilms (Yousefi et al., 2012). Considering the rich carbohydrate environment where L. kunkeei can survive, the presence of unique glucosyltransferase genes, such as gftC, could facilitate bacterial colonization of flowers, as well as the honey bee gut.

Another unique gene found in L. kunkeei, encodes for the adapter protein MecA, a pleiotropic regulator of bacterial development. This protein has been shown to affect competence, protein degradation and sporulation in Bacteria, such as Bacillus subtilis (Schlothauer et al., 2003). MecA interacts with the chaperone ClpC, and with the transcription factor ComK, promoting the degradation of this protein during the logarithmic growth phase. The degradation of ComK stops when bacteria enters to stationary growth phase, where the quorum-sensing pheromone ComX promotes the synthesis of ComS, which binds to MecA and prevents the interaction of MecA-ComK (Persuh, Mandic-Mulec & Dubnau, 2002; Prepiak et al., 2011; Wahl et al., 2014). This could have an effect on the biofilm generation capabilities of L. kunkeei MP2, but this needs to be explored experimentally.

At least seven different Lactobacillus species have been characterized in the gut microbiota of A. mellifera, where it has been suggested that they play different roles in the stability of the host functions (Engel & Moran, 2013). L. kunkeei MP2 appears to have a unique set of genes when compared to other strains of L. kunkeei, as well as with other species of Lactobacillus (Table S3), which suggest unique adaptation strategies of L. kunkeei MP2 to the gut of A. mellifera.

We also identified a hypothetical protein with similarities to a low-molecular-weight protein-tyrosine phosphatase (LMPTP), unique to the L. kunkeei MP2 genome, compared to other L. kunkeei strains and other Lactobacillus species. This LMPTP is similar to the YfkJ protein from Bacillus subtilis, which has been involved in the response to ethanol stress (Musumeci et al., 2005). Ethanol, and other organic compounds, are commonly present in the environment, and accumulate in the bacterial membrane affecting its physical-chemical properties, and in consequence, their functions (Weber & De Bont, 1996). This could suggest a better tolerance to organic compounds, such as ethanol, for L. kunkeei MP2, which could help this organism to tolerate unfavorable conditions, and have a unique competitive advantage compared to other Lactobacillus species (De Guchte van et al., 2002).

The diversity of Firmicutes species in A. mellifera could imply a metabolic diversity that could be crucial for honey bee fitness (Engel & Moran, 2013). Comparative genomics of Lactobacillus genomes, have shown that close to 45% of its accessory genome encode for proteins involved in carbohydrate metabolism and transport functions (Ellegaard et al., 2015). With this metabolic diversity found in the accessory genome, is no surprising to find unique genes in the accessory genome of L. kunkeei MP2, when compared to other strains of L. kunkeei, as well as other Lactobacillus species (Table S3). These genes encode for proteins that take part of the degradation of carbohydrates, transport of molecules, transcription, as well as membrane proteins. It is very likely that some of these genes were acquired via horizontal gene transfer from a diverse group of organisms, including those that inhabit the gut of A. mellifera. This has been observed in the adaptation of strains of Gilliamela apicola and Snodgrassella alvi to the guts of the honey bee and the bumble bee (Kwong et al., 2014), as well as in other mammalian guts (Shterzer & Mizrahi, 2015).

Integrative and conjugative elements in MP2

Multiple mobile genetic elements (MGEs), were identified in the genome of L. kunkeei MP2, including prophages, transposons, and integrases. Several of these genes were unique to the MP2 genome, compared to the other draft genomes of L. kunkeei and other Lactobacillus strains. To explore a possible association between MGEs and the unique genes found in the genome of L. kunkeei MP2, we performed a prediction of genomic islands using Island Viewer 3 (Dhillon et al., 2015). With this approach, we found that most of the unique genes are found outside genomic islands (Fig. 1, Table S3). This could suggest either events of gene loss or ancestral transfer events in the genome of L. kunkeei MP2 (Tamarit et al., 2015).

Most of the MGEs found in the genome, had similarities to integrative and conjugative elements (ICEs), which are characterized by their prophage-like mode of maintenance (Burrus et al., 2002). To contrast this result, the uniquely identified genes in the genome of L. kunkeei MP2 were compared against the ICEberg database (Bi et al., 2012) (Table S3). ICEs commonly encode for genes that provide an increased fitness to the host, such as antibiotic resistance genes, phage resistance, and heavy metal transport (Burrus et al., 2002). In the case of MP2 we found genes that have similarities to transmembrane proteins, phage-related proteins, and antibiotic resistance mechanisms, suggesting that the incorporation and stability of these unique genes in the genome of L. kunkeei MP2, is providing an increase in the fitness of this bacterial strain in the gut of A. mellifera. Among the predicted phage-like sequences, we found one coding for a mef(A)/msr(D) resistance protein, with similarity to a sequence from Streptococcus pyogenes, involved in the resistance to macrolides (Iannelli et al., 2014). In the European Union, the usage of antibiotics, and antibiotic-containing compounds, is not permitted. However, macrolides (such as tylosin and streptomycin), are still used as a preventive treatment against Paenibacillus larvae, the causal agent of American foulbrood, in many countries (Reynaldi et al., 2010; Gaudin, Hedou & Verdon, 2012). Thus, if bees were exposed to antibiotics in their diet, it is possible that the gut microbiota may have acquired the necessary molecular mechanisms to adapt and survive in an exposed environment (Tarapoulouzi et al., 2013). This acquisition can be explained by horizontal gene transfer events from the surrounding natural environment (such as soil). Here we can find multiple bacterial and fungal species that commonly produce antimicrobial compounds, and could act as a source of these resistance genes (Alippi, León & López, 2014).

Most of ICEs coding genes are usually present within genomic islands in the host genome (Hacker & Carniel, 2001; Boyd, Almagro-Moreno & Parent, 2009), but in the case of L. kunkeei MP2, none of the predicted ICEs genes were found in the context of genomic islands according to the predictions performed with IslandViewer 3 (Dhillon et al., 2015). This could suggest the presence of previously uncharacterized ICEs, or also our current limitation in the detection of ICEs from Lactobacillus species.

Prediction of horizontal gene transfer events

To predict horizontally transferred genes we used Darkhorse (Podell & Gaasterland, 2007) to analyze the complete genome of L. kunkeei MP2. We did not consider hits to organisms within the same Phylum, to avoid false predictions, although this could lead to ignore real transfer events between more closely related organisms. A total of 19 genes were predicted to have been acquired via horizontal gene transfer (Table S3), with a normalized LPI score cutoff of 0.546. Seven of these genes had matches with the genome of A. mellifera, which a detailed look suggested as a contamination of the genome of A. mellifera with sequences from Lactobacillus species. This has been previously reported for other genome projects (Merchant, Wood & Salzberg, 2014). Only one of the genes predicted to be acquired via HGT was unique to L. kunkeei MP2 when compared to other L. kunkeei strains and other Lactobacillus genomes, which codifies for a hypothetical protein, with a best hit as a phage protein from Halomonas sp. HAL1. None of the predicted genes was found associated with an ICE or a genomic island. Although the apparent lack of genes of acquired via HGT could be explained by the genome reduction that has been observed in L. kunkeei strains (Tamarit et al., 2015), an alternative explanation is the lack of reference genomes of isolates obtained from the gut of A. mellifera. This could limit our ability to predict events of horizontal gene transfer between members of the honey gut microbiota.