High-throughput sequencing analysis reveals genomic similarity in phenotypic heterogeneous Photorhabdus luminescens cell populations

Abstract

Purpose

Phenotypic heterogeneity occurs in many bacterial populations: single cells of the same species display different phenotypes, despite being genetically identical. The Gram-negative entomopathogenic bacterium Photorhabdus luminescens is an excellent example to investigate bacterial phenotypic heterogeneity. Its dualistic life cycle includes a symbiotic stage interacting with entomopathogenic nematodes (EPNs) and a pathogenic stage killing insect larvae. P. luminescens appears in two phenotypically different cell forms: the primary (1°) and the secondary (2°) cell variants. While 1° cells are bioluminescent, pigmented, and produce a huge set of secondary metabolites, 2° cells lack all these phenotypes. The main difference between both phenotypic variants is that only 1° cells can undergo symbiosis with EPNs, a phenotype that is absent from 2° cells. Recent comparative transcriptome analysis revealed that genes mediating 1° cell-specific traits are modulated differently in 2° cells. Although it was previously suggested that heterogeneity in P. luminescens cells cultures is not genetically mediated by, e.g., larger rearrangements in the genome, the genetic similarity of both cell variants has not clearly been demonstrated yet.

Methods

Here, we analyzed the genomes of both 1° and 2° cells by genome sequencing of each six single 1° and 2° clones that emerged from a single 1° clone after prolonged growth. Using different bioinformatics tools, the sequence data were analyzed for clustered point mutations or genetic rearrangements with respect to the respective phenotypic variant.

Result

We demonstrate that isolated clones of 2° cells that switched from the 1° cell state do not display any noticeable mutation and do not genetically differ from 1° cells.

Conclusion

In summary, we show that the phenotypic differences in P. luminescens cell cultures are obviously not caused by mutations or genetic rearrangements in the genome but truly emerge from phenotypic heterogeneity.

Findings

Bacteria constantly encounter different environmental stress conditions, whereupon they have evolved different survival strategies to cope with these challenges. Besides altering the expression of single genes, one of these adaptation strategies is genetic modifications such as the occurrence of DNA methylation or genomic rearrangements to evolve a different phenotype for adaptation (Smits et al. 2006). Phenotypic heterogeneity instead is another strategy, e.g., for bet-hedging to ensure the survival of a bacterial population describing the appearance of different phenotypic cells within a genetically identical population (Avery 2006; Davidson and Surette 2008; Grote et al. 2015), resulting in phase variation mostly correlated with altered gene expressions (Elowitz et al. 2002; van der Woude 2011; Davis and Isberg 2016). Examples of this adapting phenotypic heterogeneity are persister cells (Balaban et al. 2004) as well as the occurrence of competence or sporulation of the Gram-positive bacterium Bacillus subtilis (Veening et al. 2005; Smits et al. 2006). Phenotypic heterogeneity also occurs in the entomopathogenic bacterium Photorhabdus luminescens, which exists in two phenotypically different cell forms, the endosymbiotic primary (1°) cells and the free-living secondary (2°) cells. In its dualistic life cycle, the 1° cell variant colonizes entomopathogenic nematodes (EPNs), which invade insect larvae in the soil. Once inside, the EPNs release the bacteria into the haemocoel, where P. luminescens produces a huge set of toxins to effectively kill the larvae (Forst et al. 1997). During the infective cycle (also after prolonged cultivation in the laboratory), up to 50% of 1° cells switch to the 2° cell variant. The 1° cells exhibit different phenotypes such as biofilm formation, pigmentation, bioluminescence, and the production of secondary metabolites, characteristics that are absent from 2° cells (Akhurst 1980; Forst et al. 1997; Joyce and Clarke 2003; Eckstein and Heermann 2019). Moreover, 2° cells can neither reassociate with EPNs nor support their growth and development anymore and therefore remain in the soil when the EPNs have left the depleted insect cadaver. Recent studies indicated a new fate of these 2° cells in soil, as this cell variant reacts to and interacts with plant roots (Regaiolo et al. 2020). Furthermore, comparative transcriptome analysis revealed that genes responsible for 1° cell-specific phenotypes are downregulated in 2° cells (Eckstein et al. 2019). Although the exact regulation mechanism of this phenotypic switching in P. luminescens is yet unknown, some studies showed that transcriptional regulators play an important role during this event. One of these regulators is HexA, a member of the LysR-type transcriptional regulator family. HexA is involved in the phenotypic switching process of P. luminescens, by directly and indirectly repressing the expression of 1° cell-specific genes (Joyce and Clarke 2003; Langer et al. 2017). Moreover, the RNA chaperone Hfq regulates the expression of hexA mediating higher copy numbers of HexA in 2° cells, suggesting that Hfq is also involved as a global regulator in the regulation of phenotypic heterogeneity in P. luminescens cell populations (Neubacher et al. 2020). Furthermore, XRE-like transcriptional regulators were shown to also control phenotypic heterogeneity in P. luminescens. Indeed, deletion of xreR1 in 1° and xreR2 in 2° cells and insertion of extra copies of xreR1 in 2° cells and xreR2 in 1° cells led to the opposite phenotype in both cell forms (Eckstein et al. 2021). The two-component system AstS/AstR was found to control the timing of phenotypic switching in P. luminescens, since the deletion of astR led to faster switching of 1° cells compared to the wild type (Derzelle et al. 2004).

Phenotypic switching in P. luminescens has previously been referred to as phase variation (Akhurst and Boemare 1988). However, this phenomenon has been suggested to be different from the classical bacterial phase variations as both cell forms were suspected to be genetically homogeneous (Forst et al. 1997). Classical phase variation involves reversible genetic events, occurs at a significant frequency, and is almost reversible. Larger DNA rearrangements or modifications, genetic instability, or the loss of plasmids were excluded in P. luminescens 2° cell formation suggesting that the differences between 1° and 2° cells are caused by phenotypic and not genetic heterogeneity in P. luminescens (Akhurst et al. 1992; Forst et al. 1997; Hu and Webster 1998; Forst and Clarke 2002). However, none of the previous studies could provide evidence that heterogeneity in P. luminescens cell populations is due to true phenotypic heterogeneity. For that reason, we analyzed and compared genomes of both P. luminescens subs. laumondii strain DJC 1° and 2° (Zamora-Lagos et al. 2018) [later reclassified as P. laumondii (Machado et al. 2018)] to prove whether the different characteristics of P. luminescens 1° and 2° derive from phenotypic and not from genotypic heterogeneity.

The experimental workflow is schematically presented in Fig. 1. First, 1° cells were aerobically cultivated at 30 °C by shaking at 200 rpm over 11 days in LB broth [1% NaCl (w/v); 1% tryptone (w/v); 0.5% yeast (w/v)] and streaked on LB agar plates. Upon the phenotypic appearance of red pigmentation and bioluminescence, six single colonies of each cell variant were picked, bioluminescent and pigmented colonies as representatives for the 1° and dark non-pigmented colonies as representatives for the 2° variant. Then, the genomic DNA was extracted from overnight cultures using a genomic DNA extraction kit (Süd-Laborbedarf, Gauting, Germany) according to the manufacturer’s protocol, and high-throughput sequencing (HTS) analysis was performed including the laboratory strains of P. luminescens 1° and 2°. For that, library preparation of 50 ng gDNA was performed using Nextra Library Prep Kit (Illumina) according to the manufacturer’s protocol. Libraries were quality controlled with DNA High Sensitivity DNA Kit on Bioanalyzer (Agilent) and quantified on Qubit 2.0 Fluorometer (Thermo Fisher Scientific with the ds HS Assay Kit). Genome sequencing was performed in the Genomics Service Unit (LMU Biocenter, Munich) on Illumina MiSeq with v2 chemistry (2× 250-bp paired-end sequencing). The resulting HTS data were processed using different bioinformatics tools. First, reads were aligned and mapped against the Photorhabdus luminescens subs. laumondii DJC reference genome (GenBank: CP024900.1) using Bowtie 2 (v 2.3.5) and then quality filtered with SAMtools (v. 1.13), allowing alignments with mapping quality > 30 (Li et al. 2009; Langmead and Salzberg 2012). After that, read duplicates were marked via the Picard tool (v. 2.21.4) (http://broadinstitute.github.io/picard/) to avoid distort genome coverage. Qualimap 2 was used for multi-samples quality control of HTS data (Okonechnikov et al. 2015). Pairwise comparison of the genomes was performed using VarScan 2 (v. 2.4.2), where variants with a base and mapping quality > 30 were called and filtered for single nucleotide polymorphisms (SNPs) with a coverage of > 20× (Koboldt et al. 2012). The resulting data were manually inspected for informative SNPs. The sequence data were uploaded to the NCBI sequence read archive under BioProjects PRJNA812858 for 1° clones and PRJNA812795 for 2° clones (https://www.ncbi.nlm.nih.gov). For the laboratory strains of P. luminescens 1° (DJC 1°) and 2° (DJC 2°), both with a coverage of > 90%, no SNPs were detected when compared to the reference genome, whereas the mean GC content of the mapped reads ranged between 42 and 47% (Additional file 1: Table S1). During the library preparation, PCR is one of the principal sources leading to GC content bias in HTS. Indeed, diverse base composition bias in the G and C bases emerges during library preparation upon PCR (Dohm et al. 2008; Aird et al. 2011; Benjamini and Speed 2012). With a coverage of 100%, both variants are genetically identical confirming that the different appearances are due to phenotypic heterogeneity upon different gene expression. However, the sequencing data displayed various SNPs occurring in 1° cells as well as in 2° cells only after prolonged cultivation. On average, with a mean rate of 95%, all sequenced samples displayed high coverage of the reference genome (Additional file 1: Table S1). Generally, we observed only a few mutations among all tested samples, which is in accordance with the observation that P. luminescens, has the lowest mutation rate among bacteria (Pan et al. 2021). Nevertheless, all these mutations are not consistent, as all the P. luminescens 1° as well as 2° samples displayed mutations in different loci (Table 1). In all tested samples, 1–4 SNPs for 1° cells and 1–6 SNPs for 2° cells were detected. Spontaneous gene mutation during replication, presumably also after prolonged cultivation, is not a rare phenomenon; their occurrence was explained to keep a balance between the effects of deleterious mutation rate and metabolic costs (Drake et al. 1998; Denamur and Matic 2006). Generally, the mutation rate in bacteria was described to range between 1 × 10⁻⁶ and 1 × 10⁻⁸ base substitutions per nucleotide per generation (Westra et al. 2017). For example, in Escherichia coli, a deleterious mutation rate of 2–8 × 10⁻⁴ and a beneficial mutation rate of 2 × 10⁻⁹ per genome per replication have been calculated (Kibota and Lyncht 1996; Boe et al. 2000; Imhof and Schlotterer 2001), whereas recent studies reported a low mutation rate with a low base substitution rate of 5.94 × 10⁻¹¹ per nucleotide site per cell division in P. luminescens (Pan et al. 2021). Most of the genes affected from this spontaneous mutation in our analysis code for phage tail fibers in 1° as well as 2° cells. These genes are PluDJC_00175 (phage tail collar domain), PluDJC_15370 (phage tail domain), and PluDJC_15455 (phage tail fiber repeat and collar domain). The latter one occurred to have the most mutations in both cell variants, and in some cases, point mutations led to a base pair exchange not affecting the amino acid sequence. Moreover, in three from the six switched clones (2° cells), we observed single point mutations in different regions of rpoD (PluDJC_19710). However, this mutation was not observed in the 2° cell control genome and the other switched variants. Further genes or promoters displaying mutations are listed in Table 1. Additionally, some of the genes that displayed a mutation in the genome analysis have been further checked on mutations. For that purpose, single 1° and 2° colonies were again picked and aerobically cultivated for 48 h at 30 °C, and genomic DNA was extracted using a genomic DNA extraction kit (Süd-Laborbedarf, Gauting) according to the manufacturer’s protocol. Primers (Table 2) were designed to amplify the respective genes, sequenced, and aligned to the reference genome. Although no mutations were found for hexA and hfq in HTS analysis, primers for both including the promoters of both respective genes were designed to exclude potential mutations in the promoter regions affecting the gene expression. Sequencing data revealed no mutations in the rpoD, PluDJC_00175, hex and hfq promoter regions neither in 1° nor in 2° cells. For PluDJC_15370 and PluDJC_15455, (only few) mutations for both cell variants were detected. These mutations were located in different parts of the respective genes as pointed out in the HTS data, but they were not equally distributed throughout the cell variations and were not found in all switched samples. Even though these mutations are inconsistent and do not lead to genotypic heterogeneity, considering that no phenotypic cell variations have been observed, further work should be investigated to understand the higher mutation rate in P. luminescens phage-related genes (e.g., PluDJC_15370 and PluDJC_15455). Mutations in loci coding for phage subunits as well as loci involved in immunity against phages are known to have a faster mutation rate compared to point mutations in random genomic regions (Bikard and Marraffini 2012).

Fig. 1

Schematic presentation of the experimental workflow for genomic comparison between P. luminescens 1° and 2° variants. In brief, P. luminescens 1° cells were cultivated for 11 days so that a large proportion of single cells were switched to the 2° phenotype and then plated on LB agar. Then, six single colonies of each variant, 1° and 2°, were picked and cultivated; genomic DNA was isolated, analyzed by HTS, and resulting DNA sequence data was mapped against the P. luminescens ssp. laumondii DJC reference genome

Table 1 Point mutations in different gene loci of P. luminescens 1° and 2° clones identified by HTS. Sample ID 1 indicates 1° cells and 2 indicates 2° cells of P. luminescens. A–F indicate the 6 different replicates of each tested phenotypic variant. The info indicates the sum of the coverage with the reference genome (n1), the number of matches with the covered region (n2), and the number of occurring SNPs within n1 (n3) represented as n1:n2:n3

Table 2 Primers used for amplification of the respective genes carrying point mutations identified by HTS analysis. Primers were used for the gene amplification and further sequencing

Taken together, our data prove evidence that variations in P. luminescens subs. laumondii DJC cell populations is truly based on phenotypic heterogeneity. The identified mutations after long-term cultivation are due to spontaneous mutations that are randomly distributed on different genes and not always located in the same genetic area, so that genetic modifications or genomic rearrangements are not involved in phenotypic heterogeneity, i.e., phenotypic switching from 1° to 2° cells in P. luminescens cell populations.