Contrasting genes conferring short- and long-term biofilm adaptation in Listeria

Abstract

Listeria monocytogenes is an opportunistic food-borne bacterium that is capable of infecting humans with high rates of hospitalization and mortality. Natural populations are genotypically and phenotypically variable, with some lineages being responsible for most human infections. The success of L. monocytogenes is linked to its capacity to persist on food and in the environment. Biofilms are an important feature that allow these bacteria to persist and infect humans, so understanding the genetic basis of biofilm formation is key to understanding transmission. We sought to investigate the biofilm-forming ability of L. monocytogenes by identifying genetic variation that underlies biofilm formation in natural populations using genome-wide association studies (GWAS). Changes in gene expression of specific strains during biofilm formation were then investigated using RNA sequencing (RNA-seq). Genetic variation associated with enhanced biofilm formation was identified in 273 genes by GWAS and differential expression in 220 genes by RNA-seq. Statistical analyses show that the number of overlapping genes flagged by either type of experiment is less than expected by random sampling. This novel finding is consistent with an evolutionary scenario where rapid adaptation is driven by variation in gene expression of pioneer genes, and this is followed by slower adaptation driven by nucleotide changes within the core genome.

Keywords: Listeria; RNAseq; biofilm; genome-wide association study; genomics.

Fig. 1.

L. monocytogenes clones isolated in Japan are globally distributed. (a) L. monocytogenes isolates (n=868) are shown on a core-genome phylogenetic tree reconstructed using an ML algorithm implemented in RAxML. The rings, from inner to outer, indicate the phylogenetic lineage of isolates (lineage I is green, lineage II is yellow, lineage III is purple, lineage IV is red); the diversity of CCs present within each lineage; the distribution of Japanese isolates across the phylogeny is indicated in black. Bar 0.1 substitutions per site. (b) The frequency of L. monocytogenes isolates by CC. Only complexes with four or more isolates are included. Bars are coloured according to isolation location (Japan is black, the rest of the world is white).

Fig. 2.

Epidemiological factors associated with biofilm formation in Japan. Box plots showing the association between biofilm-forming ability and (a) temperature, (b) isolation location, (c) isolation source, (d) phylogenetic lineage, (e) serotype, (f) clonal complex and (g) time. Statistical analysis was performed by Student’s t-test comparing biofilm thickness in each condition (*P<0.05; **P<0.01; ***P<0.001; ****P<0.0001).

Fig. 3.

Biofilm formation is distributed throughout the clonal frame of L. monocytogenes . (a) L. monocytogenes isolates (n=108) from Japan are shown on a core-genome phylogenetic tree reconstructed using an ML algorithm implemented in RAxML. The coloured rings indicate the biofilm-forming ability of isolates at 30 °C (inner) and 37 °C (outer). Isolate Lm0132 has been highlighted using an asterisk. Bar, 0.02 substitutions per site. (b) Scatter plot of isolate biofilm formation at 30 and 37 °C. The blue dotted line indicates the line of best fit.

Fig. 4.

Genes with differential expression in transcription analysis are distinct from genes associated with biofilm formation in GWAS. (a) Differential gene expression between planktonic- and biofilm-forming states of Lm0132. The log-fold change in gene expression against its position in the F2365 genome is shown. Red lines denote a significance threshold of 100× change in gene expression. (b) and (c) Manhattan plots of unitigs associated with biofilm formation at 30 and 37 °C, respectively. The significance of each unitig k-mer’s association with biofilm formation against its position in the L. monocytogenes F2365 genome is shown. The red line denotes a genome-wide significance threshold −log₁₀(P)=6.97 (α<0.05 Bonferroni corrected with 436 125 unique variants), and the blue line denotes the suggested significance threshold −log₁₀(P)=3.97. Bar charts represent the number of genes exceeding significance thresholds with a function described in the KEGG database [76]. Dark grey bars represent gene functions (09142: Cell motility, 09145: Cellular community - prokaryotes, 09131: Membrane transport, 09132: Signal transduction, 09123: Folding, sorting and degradation, 99975: Protein processing, 09124: Replication and repair, 09121: Transcription, 09122: Translation, 09175: Drug resistance - antimicrobial, 09171: Infectious disease - bacterial, 09105: Amino acid metabolism, 09110: Biosynthesis of other secondary metabolites, 09101: Carbohydrate metabolism, 09102: Energy metabolism, 09107: Glycan biosynthesis and metabolism, 09103: Lipid metabolism, 09108: Metabolism of cofactors and vitamins, 09106: Metabolism of other amino acids, 09109: Metabolism of terpenoids and polyketides, 09104: Nucleotide metabolism), while light grey bars represent functional classes (09140: Cellular processes, 09130: Environmental information processing, 09120: Genetic information processing, 09160: Human diseases, 09100: Metabolism).

Fig. 5.

Biofilm adaptive changes in gene expression do not correlate with mutational changes identified by GWAS. Scatter plots of overlapping hits exceeding the suggestied significance threshold in GWAS and with a log(change) in gene expression >2. GWAS scores were determined as the proportion of the maximal −log(P-value) observed per experiment. RNA-seq scores were determined as the proportion of the maximal log(change) observed. Coloured dots indicate biofilm formation at 30 °C (blue) and 37 °C (red). The blue line indicates the line of best fit and dotted blue lines indicate 95 % confidence intervals.