Core genes can have higher recombination rates than accessory genes within global microbial populations

  1. Asher Preska Steinberg
  2. Mingzhi Lin
  3. Edo Kussell  Is a corresponding author
  1. Department of Biology, New York University, United States
  2. Department of Physics, New York University, United States

Abstract

Recombination is essential to microbial evolution, and is involved in the spread of antibiotic resistance, antigenic variation, and adaptation to the host niche. However, assessing the impact of homologous recombination on accessory genes which are only present in a subset of strains of a given species remains challenging due to their complex phylogenetic relationships. Quantifying homologous recombination for accessory genes (which are important for niche-specific adaptations) in comparison to core genes (which are present in all strains and have essential functions) is critical to understanding how selection acts on variation to shape species diversity and genome structures of bacteria. Here, we apply a computationally efficient, non-phylogenetic approach to measure homologous recombination rates in the core and accessory genome using >100,000 whole genome sequences from Streptococcus pneumoniae and several additional species. By analyzing diverse sets of sequence clusters, we show that core genes often have higher recombination rates than accessory genes, and for some bacterial species the associated effect sizes for these differences are pronounced. In a subset of species, we find that gene frequency and homologous recombination rate are positively correlated. For S. pneumoniae and several additional species, we find that while the recombination rate is higher for the core genome, the mutational divergence is lower, indicating that divergence-based homologous recombination barriers could contribute to differences in recombination rates between the core and accessory genome. Homologous recombination may therefore play a key role in increasing the efficiency of selection in the most conserved parts of the genome.

Editor's evaluation

Homologous recombination is an important process driving the evolution of bacterial genomes. This paper addresses how the rate of homologous recombination varies between core and accessory genes. The authors quantify recombination using an approach based on the decline of linkage between polymorphic sites over distance. The analysis indicates that core genes can be under higher rates of homologous recombination than accessory genes, which is an interesting observation that contrasts with the prevailing wisdom. Together, it suggests that homologous recombination may play a key role in increasing the efficiency of selection in conserved components of the genome.

https://doi.org/10.7554/eLife.78533.sa0

Introduction

Intrinsic to bacterial genome evolution is the process of recombination (Fraser et al., 2007; Hanage, 2016; Smith, 1991; Thomas and Nielsen, 2005), which can enhance the ability of microbes to adapt to their environment (Arnold et al., 2022; Maiden, 1998; van der Woude and Bäumler, 2004). Recombination in bacteria occurs in two forms. One is ‘homologous recombination’, in which fragments of DNA taken up from the environment are incorporated into homologous sites in the genome resulting in allele replacement and no gain or loss of chromosomal DNA (recently this has also been referred to as ‘allele transfer’ [Arnold et al., 2022]). The second is ‘nonhomologous recombination’ or ‘horizontal gene transfer’ (HGT), in which DNA fragments are inserted or removed from the genome, changing the amount of chromosomal DNA and potentially leading to gene gain or loss events (Arnold et al., 2022; Hanage, 2016). Both of these processes obfuscate clonal signals and confound phylogenetic analyses (Feil et al., 2001; Guttman and Dykhuizen, 1994; Spratt et al., 2001), and recombination is so ubiquitous in certain bacteria that whether or not bacterial genome evolution can truly be characterized as tree-like is the subject of much debate (Creevey et al., 2004; Doolittle, 1999; Sakoparnig et al., 2021; Spratt et al., 2001). Recombination rates vary greatly both between and within different bacterial species (Castillo-Ramírez et al., 2012; Chaguza et al., 2016; Chewapreecha et al., 2014; Croucher et al., 2013; Hanage, 2016; Hanage et al., 2009), and this variance is attributed to a complex interplay of selective pressure, molecular mechanisms, and ecological factors which have yet to be disentangled. Along the genome, recombination rates vary among specific genes (i.e., recombination ‘hotspots’) (Didelot et al., 2012; Everitt et al., 2014) and between gene classes such as ‘informational genes’ (those involved with transcription/translation) and ‘operational genes’ (those involved with biosynthesis, metabolism, and regulatory functions) (Jain et al., 1999; Novick and Doolittle, 2019).

This inter- and intragenomic variation in recombination rates has presented difficulties with using individual genes (e.g., 16S rRNA) to define genotypic clusters and species, and approaches such as multilocus sequence typing (MLST) have been developed to address this (Hanage et al., 2005; Hanage et al., 2006a; Maiden et al., 1998). Multilocus approaches attempt to resolve this issue by concatenating sequences from several housekeeping genes, thereby minimizing the chance that local inter- and intraspecific recombination at a single locus distorts our interpretation of clonal relationships. These housekeeping genes are part of the ‘core’ genome, which is found in nearly all strains of a given species and is thought to be less prone to recombination events (at least in the form of HGT) compared to the collections of genes found in only subsets of strains of a given species known as the ‘accessory’ genome (Daubin et al., 2002; Lan and Reeves, 2001; McInerney et al., 2017).

Quantifying the variation of homologous recombination rates across core and accessory genes is important for understanding of how selection and recombination interact to shape sequence diversity, as core genes are thought to be under strong, purifying selection (Bohlin et al., 2014; Bohlin et al., 2017; den Bakker et al., 2013; McInerney et al., 2017; Moulana et al., 2020) whereas many accessory genes are thought to be important for niche specialization and under positive or diversifying selection (Azarian et al., 2020; McInerney et al., 2017; Vernikos et al., 2015). Moreover, characterizing homologous recombination rates in niche-adaptive accessory genes will have far-reaching implications for understanding how bacteria adapt to stresses such as antibiotics, host immune responses, and nutrient conditions (Didelot et al., 2016; Povolo and Ackermann, 2019; Seifert, 1996; van der Woude and Bäumler, 2004). However, measuring homologous recombination rates across accessory genes has been stymied by the frequent gene loss and gain events that they experience, which cause difficulties in applying phylogeny-based approaches to infer homologous recombination rates (Arnold et al., 2018; Ansari and Didelot, 2014; Croucher et al., 2015; Didelot et al., 2010; Didelot et al., 2012; Didelot and Falush, 2007; Didelot and Wilson, 2015; Iranzo et al., 2019; Lefébure and Stanhope, 2007; Marttinen et al., 2012; Mostowy et al., 2017). For this reason, variation in recombination rates has primarily been studied across core genes (Didelot et al., 2012; Everitt et al., 2014; Jain et al., 1999). Whereas previous studies indicate that the core genome experiences fewer HGT events (Daubin et al., 2002; Haegeman and Weitz, 2012; Lan and Reeves, 2001; Lobkovsky et al., 2013; Wolf et al., 2016), direct comparisons of homologous recombination rates between the core and accessory genomes are currently lacking.

Here, we determine how homologous recombination rates differ between core and accessory genes. We focus on Streptococcus pneumoniae, which frequently engages in homologous recombination, has been extensively sequenced (Chaguza et al., 2016; Chewapreecha et al., 2014; Croucher et al., 2013; Hanage, 2016; Hanage et al., 2009), and is clinically relevant (O’Brien et al., 2009). For these reasons, S. pneumoniae provides a key case study for quantifying variation in homologous recombination rates across the core and accessory genome and, by analyzing additional species, we determine the generality of these trends across several bacterial species.

Results

Inferring recombination parameters for analyzed samples and unobserved gene pools using correlated substitutions

To infer the parameters of homologous recombination, we apply and extend the mcorr method (Lin and Kussell, 2017), an approach that avoids phylogenetic reconstruction by using a coalescent-based population genetics model with recombination to capture the statistics of large-scale sequencing data (see Lin and Kussell, 2019 for details). We define a ‘sample’ to be a set of lineages with mean coalescence time Tsample, which mutate and exchange homologous DNA fragments with a much larger, unobserved ‘pool’ of sequences having mean coalescence time Tpool (Figure 1A). By fitting the model to the data, we infer recombination parameters of both the sample (i.e., the specific set of sequences used for analysis) and its pool; sample construction and collection are discussed in the next section. We obtain the pool’s mutational divergence, θpool2μTpool, the pool’s recombinational divergence, ϕpool2γTpool, the sample’s recombination coverage, csample, and the mean recombined fragment size, f¯ (where μ and γ are the synonymous substitution rate and recombination rate, respectively). Together, θpool and ϕpool estimate the number of synonymous substitutions and recombination events that have occurred per site since coalescence of the pool, and csample estimates the fraction of the genome that has been replaced by homologous recombination with the pool since coalescence of the given sample. Unlike the synonymous substitutions, which we assume to be largely neutral, recombination events may be due to selective pressure or neutral drift.

Inferring parameters of homologous recombination from whole genome sequences (WGS).

(A) Schematic depicting exchange of homologous DNA fragments between the analyzed sequences (‘Sample’) and a larger, unobserved reservoir of bacterial genomes (‘Pool’). The coalescence times of the pool and sample are denoted as Tpool and Tsample , respectively. (B) Example correlation profiles of synonymous substitutions for the core genes in samples consisting of 568 S. pneumoniae and 2215 C. jejuni WGS (further description of these datasets are given in the following sections). (C) Recombination parameters inferred from fitting the profiles shown in panel C to a population genetics model (see Materials and methods). Left to right: the pool’s mutational divergence (θpool), the pool’s recombinational divergence (ϕpool), and the sample’s recombination coverage (csample). Error bars are 95% bootstrap confidence intervals (see Materials and methods). Colors correspond to panel B. θpool and ϕpool have units of bp-1 and csample is given as the percentage of genomic sites that have recombined.

The model predicts the conditional probability of a difference at genome site i+l given a difference at site i, which we refer to as the ‘correlation profile’, P(l), where l is the distance between sites in basepairs (bp). We take a set of whole genome sequences (WGS) to constitute a sample, and use alignments of protein coding (CDS) regions to measure profiles of synonymous substitutions for all possible sequence pairs, yielding an average profile for the sample. For each pair of sequences in the sample at genomic position i, a binary variable σi is assigned 1 for a difference or 0 for identity. The correlation profile is given by P(l)P(σi+l=1|σi=1), where i is restricted to third position sites of codons. While the profiles are computed using gene alignments, recombined fragments may span multiple genes, and the mean fragment size in the coalescent model can take values larger than the length of single genes. Fitting the model to these correlation profiles yields the parameters of homologous recombination described above (see Lin and Kussell, 2019 and Materials and methods for details of model and fitting); flat profiles indicate a lack of recombination, while in the presence of recombination P(l) has a monotonic decay with l, and declines more rapidly with increasing recombination (Lin and Kussell, 2017). Example profiles and inferred parameters are shown in Figure 1B, C using sets of S. pneumoniae and Campylobacter jejuni WGS; these sequences are part of the larger datasets used in the analysis that follows. Thus, the method infers parameters of both the sample and the pool of sequences with which it has recombined without using or inferring any phylogenetic relationships, offering a key advantage in determining homologous recombination rates across accessory genes.

Homologous recombination rates are higher in the core versus accessory genes for S. pneumoniae

To study recombination across the core and accessory genome of a given bacterial species, we analyze a wide range of samples from the global microbial population. In our model, each ‘sample’ is composed of collections of genomes from nature; they are statistical samples from a larger, unobserved bacterial population or gene pool. By examining a diverse set of samples from a given species, we are able to infer the distribution of recombination parameters for their unobserved gene pools. However, as many genomes are collected from specific geographic regions or individual hospitals (Chaguza et al., 2016; Chewapreecha et al., 2014; Pelton et al., 2007), a sample can be artificially overweighted by a single collection site accounting for many genomes. Therefore, to obtain a more meaningful set of samples for a given species, and ameliorate some of the compositional biases inherent in nonrandom sampling, we use sequence clusters obtained by clustering WGS solely based on sequence similarity using the average linkage algorithm (Wheeler and Kececioglu, 2007) with the pairwise synonymous diversity (ds) as the distance metric. Such approaches are widely used to define well-resolved genotypic clusters within bacterial populations without the inference of phylogenetic relationships (Hanage et al., 2006a). We then analyze single clusters or pairs of clusters to infer the distribution of pool parameters. For each cluster or cluster pair, we infer both the divergences of the unobserved gene pool with which they recombine (i.e., θpool and ϕpool), as well as the amount of recombination that has taken place within the set of analyzed sequences (i.e., the cluster or cluster pair) since coalescence (csample). For single clusters, we fit correlation profiles measured by averaging over all of a given cluster’s sequence pairs, which yields the parameters of the unobserved pool that the cluster interacts with. For cluster pairs, we fit correlation profiles measured by averaging exclusively over sequence pairs consisting of one sequence from each of the two clusters. By analyzing all possible pairs of clusters, we use diverse sets of sequences as our samples, yielding a parameter distribution that accounts for the full range of potential interactions between different sequences, regardless of their distance in sequence space and agnostic to population structure or sampling biases.

To sample the global population of S. pneumoniae, we used the PubMLST database (see Supplementary file 5) which includes WGS of strains from across the world. For each WGS, we aligned all sequencing reads to a reference genome to create an alignment of all CDS regions (see Materials and methods for details). We then measured ds across all genes (i.e., core and accessory) for each strain pair, clustered based on these distances to make a dendrogram, and split the sequences into flat clusters, where no two sequences within a cluster were more distant than the 10th percentile of pairwise distances, which corresponded to ds=0.015 (Figure 2A, Materials and methods for details; mean ds values within and between clusters for core and accessory genes are given in Supplementary file 3). This resulted in 44 major clusters (where a major cluster has >100 strains) which we used as our samples (946 total cluster pairs; Figure 2B shows an example correlation profile from a cluster), and we inferred recombination parameters for the core (defined here as present in >95% of strains) and accessory genes of each of the samples and the unobserved pools with which they recombine (Figure 2C–E).

Figure 2 with 4 supplements see all
Inference of recombination parameters for the core and accessory genome of S. pneumoniae.

(A) Dendrogram resulting from hierarchical clustering of 26,599 whole genome sequences (WGS) from the PubMLST genome collection for S. pneumoniae using the average linkage algorithm. The dendrogram was cut at the 10th percentile of all measured pairwise distances (ds ~ 0.015) yielding a discrete set of flat clusters. Here, we analyzed the 44 major clusters (those with >100 sequences) which resulted from this cut, encompassing 24,097 strains. In the dendrogram, alternating colors delineate adjacent clusters. (B) Correlation profiles measured across the core and accessory genes of a S. pneumoniae sequence cluster. Open-faced circle is the profile measured from sequencing data and solid line is the fit to the population genetics model. (C–E) Distributions of the pool’s mutational and recombinational divergence (θpool and ϕpool , shown in panels C and D, respectively) and the relative recombination rate of the pool ((γ/μ)pool , shown in panel E). For C–E, the bottom plot depicts empirical cumulative distribution functions (ECDFs) for each parameter and the top plots depict the medians of the bottom plot, where error bars are 95% bootstrap CIs created by sampling the distributions with replacement (n = 1000). Each step in the ECDF corresponds to a pool recombination parameter inferred from a correlation profile measured over a sequence cluster or pair of clusters like that shown in panel B. Core genes are defined as genes found in >95% of strains. We used model selection with the Akaike information criterion to ensure that each profile was well fit (see Materials and methods for details). θpool and ϕpool have units of bp-1 and (γ/μ)pool is unitless. Two-sided p-values were calculated using the Wilcoxon signed-rank test and were p = 1.9e−159, 3.9e−53, and 3.5e−99 for panels C–E, respectively.

Figure 2—source data 1

List of SRA accession numbers corresponding to the raw reads of genomes used from NCBI.

https://cdn.elifesciences.org/articles/78533/elife-78533-fig2-data1-v2.txt

We found that the core and accessory genes of S. pneumoniae have distinct distributions of θpool and ϕpool (Figure 2C, D), and the median of θpool is lower in core versus accessory genes, indicating that core genes are less mutationally diverged than accessory genes. Moreover, we observed that ϕpool is higher for the core genes, implying that core gene pools have experienced a greater number of recombination events per site compared to accessory genes (Figure 2D). We examined this further using the relation ϕpool/θpool=γ/μ, which removes the dependence on coalescence time to yield the relative recombination rates of the unobserved gene pools (from hereon we will append the subscript pool to indicate this, i.e., (γ/μ)pool). The inferred distribution of (γ/μ)pool shows that for S. pneumoniae, the core genome has a higher relative recombination rate than the accessory genome (Figure 2E), again indicating that core genes recombine more frequently. While the inferred distribution of recombination coverage for the samples initially suggests that the accessory genome may recombine more (Figure 2—figure supplement 1A), this can be reconciled by considering the mean size of the recombined fragments (f¯), which is higher for the accessory genes (Figure 2—figure supplement 1B). Taken together, we find that while accessory genes incorporate larger fragments, the core genes of the sampled genomes experience more recombination events on a per site basis (Figure 2—figure supplement 1C).

We tested the robustness of these results in several ways. As there are various thresholds used to define ‘core’ genes (Livingstone et al., 2018; McInerney et al., 2017; Page et al., 2015; Vernikos et al., 2015), we tested how the parameter distributions shifted with different thresholds and found similar trends (Figure 2—figure supplement 2). To determine whether the difference in the abundance of the core and accessory genes (i.e., how many strains a gene appears in) influences our inference of recombination parameters, for each cluster we artificially ‘thinned’ the core genome by replacing each accessory gene sequence with a randomly chosen core gene sequence. This yielded a thinned core genome where each core gene had the same number of sequence pairs as the accessory genome. We found that this thinned core genome showed similar recombination trends as the actual core genome (Figure 2—figure supplement 3).

Aligning raw reads to a reference genome to build consensus genomes allows for access to a much larger set of WGS, as there are comparably fewer whole genome assemblies for S. pneumoniae, which are typically a requirement for reference-free alignment (Ding et al., 2018; Lin and Kussell, 2019; Page et al., 2015) (as of the time of this analysis NCBI GenBank had 81 complete genome assemblies for S. pneumoniae versus the ~26,000 WGS accessed through PubMLST). However, a potential drawback of aligning reads to a single reference genome is that a subset of accessory genes will be missed. New approaches for aligning to reference graphs appear promising in this regard, yet remain computationally expensive for large numbers of reads (Colquhoun et al., 2021). We therefore sought to address this issue by testing how the parameters changed when more accessory genes were included by building a reference pangenome from all the genome assemblies in NCBI GenBank using Roary (Page et al., 2015) and aligned all the raw reads from PubMLST to this (see Materials and methods). Pangenome analysis suggests that this number of genomes should encompass the majority of genes in the S. pneumoniae pangenome (Donati et al., 2010). While the pangenome alignment included more genes (6200 versus 2018 genes), we found very similar trends for the inferred recombination parameters to those inferred using a single reference alignment (Figure 2—figure supplement 4).

To assess the generality of the trends in recombination parameters between the core and accessory genome, we analyzed 11 additional microbial species (encompassing >100,000 genomes including those already analyzed for S. pneumoniae). We used the same procedure described above to cluster sequences and measure correlation profiles for clusters and cluster pairs, and inferred distributions of recombination rates for the unobserved pools for these species. Pairing the inferred recombination rates of core and accessory genomes for each pool, we examined the effect sizes for each species using Cohen’s d statistic (Cohen, 1988; Nakagawa and Cuthill, 2007) to determine the magnitude and direction of the effect (Figure 3A; Supplementary file 1). We found that, in line with our observations in S. pneumoniae, certain species (e.g., Escherichia coli, Shigella flexneri, Neisseria meningitidis, and C. jejuni) showed markedly higher recombination rates for core genes, whereas others showed slightly higher recombination rates in the core genome. Helicobacter pylori had higher recombination rates in the accessory genome which were statistically significant, yet the magnitude of the effect was small (Cohen’s d = 0.18). Among the nine species for which Cohen’s d was nonzero (within 95% confidence intervals) all but H. pylori showed higher recombination rates (d < 0) in the core versus accessory genome. Similarly, of the 10 species with significantly different median recombination rates in core and accessory genomes (within 95% confidence intervals) (Supplementary file 1), all but H. pylori and N. gonorrhoeae showed higher median recombination rates in the core versus accessory genome. Different metrics thus yield largely consistent results regarding differences in recombination rates between core and accessory genomes across species. Lastly, to understand if core gene pools are less diverged than accessory pools as we had found for S. pneumoniae, we also assessed differences in θpool between the core and accessory genome for these same 12 species (Figure 3B; Supplementary file 2). We found that a majority of species had accessory gene pools which were markedly more diverged than core genes, and the few exceptions (e.g., N. gonorrhoeae and N. meningitidis) had small effect sizes.

Figure 3 with 1 supplement see all
Effect sizes (Cohen’s d) for recombination rates and mutational divergence for S. pneumoniae and 11 additional microbial species.

Effect sizes (Cohen’s d) for (A) the relative recombination rate of the pool ((γ/μ)p) and (B) the mutational divergence of the pool (θp) for core/accessory genome pairs for individual clusters and pairs of clusters for S. pneumoniae and 11 additional microbial species. Cohen’s d for paired samples was calculated as the mean paired difference (Xpacc-Xpcore, where X=(γ/μ)p or θp) divided by the standard deviation of paired difference (σ). Error bars are 95% bootstrap CIs, calculated by sampling the distributions with replacement 10,000 times. All effect sizes are listed in Supplementary files 1 and 2, as well as the medians of the distributions and the results of the Wilcoxon signed-rank test for each species. For A, the effect size and 95% bootstrap CI for S. flexneri are −25.5 and [−34.4,−21.1], respectively. When model selection was performed with Akaike information criterion (AIC) (described in Materials and methods) if part of a core/accessory pair was poorly fit, the paired sample was excluded. Full species names, number of strains, number of clusters, and mean synonymous diversity within and between clusters are reported for each species in Supplementary file 3. Legends of supplementary material.

Figure 3—source data 1

Lists of SRA accession numbers corresponding to the raw reads of genomes used from NCBI.

https://cdn.elifesciences.org/articles/78533/elife-78533-fig3-data1-v2.zip

Gene frequency as a recombination barrier in the core and accessory genome

One testable hypothesis consistent with our observation that core genes generally undergo more homologous recombination when compared with accessory genes is that the least frequent genes experience less recombination because they have fewer recombination partners. In other words, gene frequency could act as a recombination barrier. We tested this by binning the genes for the microbial species in Figure 3 into four gene frequency classes consistent with prior delineations of the pangenome as follows (q is frequency across strains) (Koonin and Wolf, 2012; Page et al., 2015): ‘cloud’ genes (q15%), ‘shell’ genes (divided into two segments, 15%<q55%,  55%<q95%), and core genes (q>95%). We then inferred the distributions of recombination parameters within each gene class for each of the 12 microbial species (Figure 3—figure supplement 1A, B). We found that certain species displayed trends consistent with a gene frequency barrier, particularly for the recombination coverage (Figure 3—figure supplement 1A). When considering the recombination rates (displayed as empirical cumulative distribution functions here because of disparate rates across and within species), C. jejuni, S. agalactiae, and P. aeruginosa showed trends most consistent with this hypothesis. For some of the exceptions in which cloud genes showed high recombination rates (e.g., N. gonorrhoeae, N. meningitidis, E. coli, and S. flexneri), this could be related to the high gene replacement rate experienced by rare genes such as ‘ORFans’ (Wolf et al., 2016), or diversifying selection experienced by these genes. Overall, we see a modest positive correlation with both recombination coverage and rates when considering the 50th and 25th percentiles of the matched parameter distributions for all species (Figure 3—figure supplement 1C, D).

Discussion

There is appreciable interest in understanding why microbes have pangenomes and how different parts of the genome have evolved (Lobkovsky et al., 2013; McInerney et al., 2017; Vernikos et al., 2015; Wolf et al., 2016). While it is known that homologous recombination plays a key role in shaping the genome (Fraser et al., 2007; González-Torres et al., 2019; Hanage, 2016; Iranzo et al., 2019; Thomas and Nielsen, 2005), large-scale analysis of this aspect of genome evolution has been limited. This was due to both computational bottlenecks and a reliance on phylogenetic methods (Arnold et al., 2018; Ansari and Didelot, 2014; Croucher et al., 2015; Didelot et al., 2010; Didelot and Falush, 2007; Didelot and Wilson, 2015; Marttinen et al., 2012; Mostowy et al., 2017), the latter of which hindered determination of recombination parameters for accessory genes, whose phylogenies are challenging to ascertain. Here, we expanded our non-phylogenetic, computationally efficient mcorr framework to overcome these obstacles, and present a detailed analysis of the variation in recombination rates between core and accessory genes for S. pneumoniae, along with several other microbial species.

We found that despite the high gene replacement rates and sequence variability experienced by some accessory genes (Iranzo et al., 2019; Kuenne et al., 2013; Kuhn et al., 2006; Wolf et al., 2016), core genes experience markedly higher recombination rates than the accessory genes in S. pneumoniae. In 8 out of 12 microbial species, we observe higher recombination rates in core genes relative to accessory genes. In a subset of microbial species, the magnitude of the effect is large (e.g., S. pneumoniae, E. coli, S. flexneri, N. meningitidis, and C. jejuni). While some work has previously shown housekeeping (Vos and Didelot, 2009) and core genes (Lin and Kussell, 2017; Park and Andam, 2020) can experience extensive recombination consistent with our results, our work offers a direct, quantitative comparison between the core and accessory genomes spanning multiple species and >100,000 sampled genomes. In certain species, mechanisms that give rise to higher recombination rates in the core genome are known; for example, some members of the Neisseriaceae (e.g., N. meningitidis and N. gonorrhoeae) and Pasteurellaceae families (e.g., H. influenzae) have DNA uptake sequences and uptake signal sequences accumulated in their core genomes which promote homologous recombination (Frye et al., 2013; Treangen et al., 2008).

We found that in S. pneumoniae, core gene pools had lower mutational divergence relative to accessory gene pools. Furthermore, we found that for the majority of the species we analyzed this was the case, and 9 out of 12 species analyzed had large effect sizes indicating core gene pools have lower mutational divergence (Cohen’s d = 0.92–2.8). Previous studies indicate that the core genome is under strong, purifying selection (Bohlin et al., 2014; Bohlin et al., 2017; den Bakker et al., 2013; McInerney et al., 2017; Moulana et al., 2020), and it is well known that purifying selection reduces mean coalescence times (Charlesworth et al., 1993; Nicolaisen and Desai, 2012). An additional consideration is that some accessory genes are used for niche specialization (Azarian et al., 2020; McInerney et al., 2017; Vernikos et al., 2015); if different alleles of the same accessory gene are adaptive to different niches, this could result in allelic diversity due to balancing selection, which would increase coalescence times of accessory genes. The lower value of θpool in core genes is therefore consistent with shorter coalescence times of core relative to accessory genes. While another possible explanation could involve codon usage biases (Plotkin and Kudla, 2011) causing differences in synonymous substitution rates between core and accessory genes, the selective coefficients for this process are tiny (s ~ 1/Ne Hershberg and Petrov, 2008); Ne=3×108 for S. pneumoniae (Bobay and Ochman, 2018) and thus effectively neutral over the timescales that we can observe in these data (i.e., coalescence times of S. pneumoniae samples and pools). Accordingly, evidence from long-term microbial evolution experiments has shown that on average the synonymous substitution rate is fairly homogeneous across the core and accessory genomes (Maddamsetti et al., 2015; Maddamsetti et al., 2017). Therefore, we attribute the difference in θpool between core and accessory genes to the difference in their coalescence times, which is consistent with core genes being under stronger purifying selection, and additionally with the possibility that a subset of accessory genes are under diversifying selection. We note that because higher levels of recombination in core genes are expected to reduce the effects of background selection (Charlesworth et al., 1993), the relative reduction in θpool observed in core versus accessory genes of S. pneumoniae and other species likely underestimates the difference in strength of purifying selection acting on these gene classes.

Experiments in numerous species show that homologous recombination rates decay with increasing sequence divergence (Fraser et al., 2007; Majewski, 2001; Majewski et al., 2000; Vulić et al., 1997; Zawadzki et al., 1995), and the role of this recombination barrier in bacterial speciation has been analyzed (Falush et al., 2006; Hanage et al., 2006b). This effect is primarily ascribed to the ubiquity of RecA and mismatch repair systems, both of which are essential for recombination and whose molecular mechanisms depend on sequence similarity. Consistent with these experimental and analytical results, the higher recombination rates and lower pool divergence we observe here for the core genome of S. pneumoniae and several other microbial species may suggest that the increased homology of the core genes reduces the barrier for homologous recombination in this part of the genome. This is also in line with modeling, which has suggested that homologous recombination slows sequence divergence in the core genome by homogenizing this part of the genome, whereas accessory genes do not have the opportunity to be homogenized as they are subject to continual gene gain and loss (Iranzo et al., 2019). Thus, we propose that an evolutionary feedback loop may operate in certain species of bacteria: purifying selection acting on core genes causes higher levels of homology at these loci, this increases their homologous recombination rates, which in turn enables more efficient purifying selection by breaking linkages between genes, effectively minimizing Hill–Robertson interference and hitchhiking effects (Barton, 1995; Barton, 2010; Felsenstein and Yokoyama, 1976) across the core genome. Further investigations are needed to examine whether such an evolutionary mechanism can act to fine-tune recombination rate variation across microbial genomes. Moreover, our work suggests that another recombination barrier for accessory genes may be related to their lower abundance in the population, as we found that the levels of homologous recombination are correlated with gene frequency in a subset of species.

In summary, we have presented a quantitative analysis of the variation in recombination rates between core and accessory genomes. As these two parts of the genome are under different forms of selective pressures, this work contributes broadly to our understanding of how selection and recombination act to shape diversity. The expansion of this approach to more specific gene classes could lead to a better understanding of how selective pressure and homologous recombination interplay to shape niche-adaptive genes such as those involved with drug resistance and antigen variation, for which allele shuffling due to homologous recombination can play a central role in their evolution (Bowler et al., 1994; Maiden, 1998; Seifert, 1996).

Materials and methods

Data and code availability

Request a detailed protocol
  • Lists of SRA accession numbers corresponding to the raw reads used to build the multi-sequence alignments analyzed here are included as source data. All SRA files, reference genomes, and complete genome assemblies are available through NCBI. All sequence collections used are listed in Supplementary file 5. For the PubMLST sequence collections, PubMLST was used to identify WGS (by filtering for strains in the ‘Genome Collection’ of each species where the sequence length is at least that of the reference genome), then the raw reads were downloaded from NCBI using their SRA numbers. Accession numbers for reference genomes used for each microbial species are also listed in Supplementary file 5.

  • All original code has been deposited at GitHub and is publicly available (links provided throughout Methods).

Description of coalescent-based population genetics model with recombination

Request a detailed protocol

The details of the model derivation and parameters are given in the ‘Supplementary Notes’ of Lin and Kussell, 2019. Here, we provide the model equations necessary for the fitting procedure used in this work.

The measured synonymous diversity of the sample, ds, is a linear combination of the pool diversity, dp d(θp) (due to recombination), and the clonal diversity of the sample, d(θs), as follows:

(1) ds=csd(θp)+(1cs)d(θs),

where cs is the recombination coverage in the sample, and d(θ) is given by the classic population genetics expression for heterozygosity (Wakeley, 2009) (see table below). The predicted form of the correlation profile is

(2) P(l)=cs,0(l)d(2θs)+cs,1(l)dp+cs,2(l)Qp(l)/ds,

where the functional forms of cs,i(l) (i=0,1,2) and Qp(l) are given in the table below. As described in Lin and Kussell, 2019, we can reexpress Equation 2 in terms of the three parameters θs, ϕs, and f¯, determined by fitting the profiles. From these, the pool parameters are obtained using the following relations: θp=ds(1+ ϕswf¯+θsa~)θs(1dsa~)(ϕswf¯+θsa~)dsa~ and ϕp=θpϕs/θs ; see table below for values of constants a~ and w.

TermExpressionDescription
θs or θpMutational divergence of the sample or pool
ϕs or ϕpRecombinational divergence of the sample or pool
f¯Mean fragment size of homologous recombination
lDistance (bp) between two loci
The following constants are used below: a~=4/3 and w=2/3.
d(θ)θ1+θa~Synonymous diversity (heterozygosity) for divergence θ
csϕswf¯1+θsa~+ϕswf¯Recombination coverage of the sample
cs,1(l)2ϕswl1+2θsa~+ϕsw(f¯+l)Probability that the most recent event was a recombination affecting only one locus
cs,2(l)ϕsw(f¯l)1+2θsa~+ϕsw(f¯+l)Probability that the most recent event was a recombination affecting both loci
cs,0(l)1cs,1cs,2Probability that the most recent event at either locus was not a recombination
Qp(l)2(θp1+θpa~)2(1+θpa~+ϕpwl1+2θpa~+2ϕpwl)Probability of a difference at both loci in the pool

Calculation of correlation profiles

Request a detailed protocol

For a given sample of aligned sequences, we measured the substitution profile, σi(k,g), for each pair of sequences k, at each position i along each gene g, where σi(k,g) = 1 for a difference and σi(k,g) = 0 for identity. In all calculations below, we only consider positions i that are third position sites of codons yielding synonymous substitutions. We computed the pairwise synonymous diversity of each gene g as:

(3) ds(g)=σi(k,g)¯

where the bar denotes averaging over all sequence pairs k and the bracket denotes averaging over all positions i. The joint probability of synonymous substitutions for two sites separated by a distance l was calculated for each gene as:

(4) Qs(l,g)=σi(k,g)σi+l(k,g)¯

We averaged these over all genes to obtain the sample diversity, ds=(1/ng)gds(g), and the joint probability of substitutions at two sites in the sample, Qs(l)=(1/ng)gQs(l,g), where ng is the total number of genes. The correlation profile is then given by P(l)=Qs(l)/ds.

For calculations of within-pool parameters, all possible sequence pairs within a cluster were considered. This was done with the original command-line (CLI) program mcorr-xmfa program in the mcorr package which takes as input a single extended multi-fasta (XMFA) file. For calculations of between-pool parameters using pairs of clusters, only sequence pairs where each sequence was from a different cluster were considered. This was done with the CLI program created for this paper called mcorr-xmfa-2clades, which uses two XMFA files (one from each sequence cluster) as inputs. These can be found in the mcorr GitHub repository: https://github.com/kussell-lab/mcorr.

Fitting of correlation profiles and model selection using Akaike information criterion

Request a detailed protocol

The basic fitting procedure used was previously described in Lin and Kussell, 2019. We used the python package LMFIT (Newville et al., 2014) version 0.9.7 (https://lmfit.github.io/lmfit-py/) to fit P(l) to the analytical form given in Equation 2; in this work, we used the least-squares minimization with Trust Region Reflective method instead of the default Levenberg–Marquardt algorithm and increased the maximum number of function evaluations from 104 to 106. As described in the text, as an additional test of goodness of fit, we compared the fit of the data with Equation 2 where all parameters vary freely (which we refer to as the ‘full-recombination model’), to the ‘null-recombination model’ in which we set cs,1=cs,2=0, which yields P(l)=d(2θs)=2θs/(1+2θsa~), that is, a constant correlation profile which is fit by taking the average over l of the measured values, yielding a single parameter, θs. We perform model selection between the null- and full-recombination models by evaluating the Akaike information criterion (AIC) for each model, then computing the Akaike weight for each model, which can be roughly interpreted as the probability that a given model yields the best prediction of the data (Bois, 2020; Wagenmakers and Farrell, 2004). The AIC was computed using LMFIT using AIC=n ln(χ2/n)+2Nv , where n is the number of data points, Nv is the number of parameters being varied in the fitting, and χ2=i=1nResidi2 is the sum of the squared residuals for each data point i. The Akaike weight for model j{,r} can then be calculated as:

wj=e0.5(AICjAICmin)e0.5(AICAICmin)+e0.5(AICrAICmin) 

where AIC and AICr denote the AIC values for the null- and full-recombination models, respectively, and AICmin=min{AIC ,AICr}. The evidence ratio, which corresponds to the likelihood of one model being favored over the other in terms of minimizing Kullback–Leibler discrepancy (Wagenmakers and Farrell, 2004), can be computed as wr/w , where wr and w are the Akaike weights of the full- and null-recombination models, respectively. Fits to correlation profiles where wr / w < 100 were not included in the analyses presented in the main text. We also did not include fits to correlation profiles which yielded unphysical parameter values (θpool<0, ϕpool<0), or extreme values of recombination coverage where parameter fitting becomes less reliable due to too few recombination events (csample<1%) or the over-fitting of a flat profile (csample>99%). The CLI program mcorrFitCompare was built for this study to fit correlation profiles to both models and can be found in the mcorr GitHub repository: https://github.com/kussell-lab/mcorr.

Clustering using pairwise synonymous diversity

Request a detailed protocol

In this study, we extended the mcorr method to measure the synonymous diversity separately for each sequence pair, denoted by dpair. This is the same as the calculation of ds described in ‘Calculation of correlation profiles’ without averaging over sequence pairs k. The CLI program mcorr-pair measures dpair for all possible sequence pairs from an XMFA file. We also built a CLI program called mcorr-pair-sync, which calculates a subset of all pairwise diversities from a given set of sequences. We primarily relied on the latter program, as it allows for parallelization of jobs on a high-performance computing cluster. The CLI programs mcorr-dm and mcorr-dm-chunks were built to collect outputs from mcorr-pair and mcorr-pair-sync, respectively, and write them to a square distance matrix for use with standard clustering algorithms. Lastly, we built a program (clusterSequences) relying on the python scipy package to cluster sequences using the unweighted pair group method with arithmetic mean or average linkage method with dpair as the distance metric. All necessary code to perform these calculations can be found in the github repository: https://github.com/kussell-lab/mcorr-clustering.

Generation of multisequence alignments for core and accessory genes using WGS

Request a detailed protocol

For all collections of WGS (Supplementary file 5), with the exception for H. pylori (where the multisequence alignment was already assembled), we used reference-guided alignment to build consensus genomes from raw reads for each sequence, then extracted the CDS regions of each gene to make extended multi-fasta (XMFA) files. For the collections from the PubMLST database, we identified WGS by filtering for all sequences in the ‘Genome Collection’ for a given organism which had lengths greater than or equal to the length of the reference genome. We then exported tables which included the corresponding SRA numbers for each strain, and downloaded the reads from NCBI using the corresponding SRA. For the other sequence collections used, we used their associated NCBI Bioproject ID to obtain the SRA numbers for each strain, and downloaded reads in the same way. Reads were mapped against a reference genome (listed in Supplementary file 5 for each microbial species) using SMALT (version 0.7.6; https://www.sanger.ac.uk/tool/smalt-0/), and consensus genomes were created using SAMtools mpileup (Li et al., 2009). We then extracted the alignments of CDS regions from these consensus genomes using genomic coordinates provided from the general feature format (GFF) file of the reference genome to create an XMFA file of CDS regions for each gene. We filtered out any gene alignments where the gene sequence was 2% gaps. We then measured the percentage of the entire strain collection which had each gene (based on presence/absence of a gene sequence) to determine if a gene should be considered core or accessory. Code to perform this analysis is provided in the GitHub repository: https://github.com/kussell-lab/ReferenceAlignmentGenerator.

Generation of multisequence alignment to pangenome reference genome for S. pneumoniae

Request a detailed protocol

We downloaded all complete genome assemblies from NCBI for S. pneumoniae as of December 18, 2020 (81 genomes), and used Prokka (Seemann, 2014) to reannotate the assemblies and generate general feature format version 3 (GFF3) files for use with Roary (Page et al., 2015). Roary was then used to generate a multi-FASTA for each gene CDS region in the pangenome (which we refer to as the ‘pangenome reference’). We aligned reads to this pangenome reference in the same manner described in ‘Generation of multisequence alignments for core and accessory genes using WGS’, then collected each gene alignment to create an XMFA file for the pangenome. We split the XMFA into files for the core and accessory genome by measuring the percentage of strains which had each gene. In this case, because alignment quality was not as high as when we aligned reads to a single reference genome (as described in the main text), we counted a gene as present if there was only a partially aligned sequence for the gene, and absent if there was no aligned sequence. Code to perform this analysis can be found in the following GitHub repository: https://github.com/kussell-lab/PangenomeAlignmentGenerator.

Generation of artificially ‘thinned’ core genome alignment for S. pneumoniae

Request a detailed protocol

To create the artificially ‘thinned’ core genome alignment described in the main text, we took the core and accessory genome alignments for S. pneumoniae generated as described in ‘Generation of multisequence alignments for core and accessory genes using WGS’, and we replaced each of the accessory gene alignments for each sequence cluster with randomly selected core gene alignments from the same cluster until all accessory gene alignments were replaced. A CLI program called ReduceCoreGenome is used to generate these thinned core genome alignments (https://github.com/kussell-lab/ReferenceAlignmentGenerator).

Statistical analysis

Request a detailed protocol

For the 95% bootstrap CIs appearing in Figure 1C, the same procedure was used as in Lin and Kussell, 2019. In brief, bootstrap replicates of the set of genes were created by resampling the list of all genes in the genome with replacement, then recalculating ds, Qs, and P(l) for each replicate. We then performed the same fitting procedure which was done on the actual dataset to infer recombination parameters on each of the bootstrap replicates to generate 95% confidence intervals. The resampling was done 1000 times.

For all other 95% bootstrap CIs, the procedure used for bootstrapping can be found in the figure legends of the corresponding figures. The process of model selection using AIC is described in ‘Fitting of correlation profiles and model selection using Akaike information criterion’. The number of major clusters and sequences used for each microbial species is given in Supplementary file 3.

Data availability

Lists of SRA accession numbers corresponding to the raw reads used to build the multisequence alignments analyzed in this manuscript are included as Figure 2—source data 1 and Figure 3—source data 1. All SRA files, reference genomes, and complete genome assemblies are available through NCBI. All sequence collections used are listed in Supplementary file 5. For the PubMLST sequence collections, PubMLST was used to identify whole genome sequences (by filtering for strains in the 'Genome Collection' of each species where the sequence length is at least that of the reference genome), then the raw reads were downloaded from NCBI using their SRA numbers. Accession numbers for reference genomes used for each microbial species are also listed in Supplementary file 5. All original code has been deposited at GitHub and is publicly available at https://github.com/kussell-lab (copies archived at swh:1:rev:4adea90557f1e592123e073b46a859d879143a2a, swh:1:rev:1d34860183a5dd55332e08edf9503c4ce6daebdf, swh:1:rev:49bf399a9d1ceae722c0c8c4aeb8376f2644d40f, and swh:1:rev:2228e8ee2df7339d28ef8b9107381d2e4767ac11).

The following previously published data sets were used
    1. Thorell K
    (2018) Dryad Digital Repository
    Data from: Rapid evolution of distinct Helicobacter pylori subpopulations in the Americas.
    https://doi.org/10.5061/dryad.8qp4n

References

  1. Book
    1. Azarian T
    2. Huang IT
    3. Hanage WP
    (2020) Structure and Dynamics of Bacterial Populations: Pangenome Ecology
    In: Tettelin H, Medini D, editors. The Pangenome: Diversity, Dynamics and Evolution of Genomes. Springer International Publishing. pp. 115–128.
    https://doi.org/10.1007/978-3-030-38281-0
    1. Barton NH
    (2010) Mutation and the evolution of recombination
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 365:1281–1294.
    https://doi.org/10.1098/rstb.2009.0320
  2. Book
    1. Cohen J.
    (1988)
    Statistical Power Analysis for the Behavioral Sciences
    Lawrence Erlbaum Associates.
    1. Falush D
    2. Torpdahl M
    3. Didelot X
    4. Conrad DF
    5. Wilson DJ
    6. Achtman M
    (2006) Mismatch induced speciation in Salmonella: model and data
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 361:2045–2053.
    https://doi.org/10.1098/rstb.2006.1925
    1. Hanage WP
    2. Fraser C
    3. Spratt BG
    (2006a) Sequences, sequence clusters and bacterial species
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 361:1917–1927.
    https://doi.org/10.1098/rstb.2006.1917
    1. Hanage WP
    2. Spratt BG
    3. Turner KME
    4. Fraser C
    (2006b) Modelling bacterial speciation
    Philosophical Transactions of the Royal Society of London. Series B, Biological Sciences 361:2039–2044.
    https://doi.org/10.1098/rstb.2006.1926
    1. Smith JM
    (1991) The population genetics of bacteria
    Proceedings of the Royal Society of London. Series B 245:37–41.
    https://doi.org/10.1098/rspb.1991.0085
  3. Book
    1. Wakeley J
    (2009)
    Coalescent Theory: An Introduction (1st ed)
    Macmillan Learning.

Article and author information

Author details

  1. Asher Preska Steinberg

    Department of Biology, New York University, New York, United States
    Contribution
    Conceptualization, Resources, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8694-7224
  2. Mingzhi Lin

    Department of Biology, New York University, New York, United States
    Contribution
    Conceptualization, Software, Investigation
    Competing interests
    No competing interests declared
  3. Edo Kussell

    1. Department of Biology, New York University, New York, United States
    2. Department of Physics, New York University, New York, United States
    Contribution
    Conceptualization, Resources, Formal analysis, Supervision, Funding acquisition, Validation, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    edo.kussell@nyu.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0590-4036

Funding

National Institutes of Health (R01-GM097356)

  • Edo Kussell

Simons Foundation (Simons Foundation Awardee of the Life Sciences Research Foundation)

  • Asher Preska Steinberg

The funders had no role in study design, data collection, and interpretation, or the decision to submit the work for publication.

Acknowledgements

This work was supported in part by NIH grant R01-GM-097356. Asher Preska Steinberg is a Simons Foundation Awardee of the Life Sciences Research Foundation. We gratefully acknowledge Nobuto Takeuchi for useful discussions and feedback on the manuscript; the New York University (NYU) high-performance computing cluster for resources, and its staff for technical support.

Version history

  1. Preprint posted: September 14, 2021 (view preprint)
  2. Received: March 10, 2022
  3. Accepted: June 30, 2022
  4. Accepted Manuscript published: July 8, 2022 (version 1)
  5. Version of Record published: September 5, 2022 (version 2)

Copyright

© 2022, Preska Steinberg et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 3,606
    views
  • 515
    downloads
  • 14
    citations

Views, downloads and citations are aggregated across all versions of this paper published by eLife.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Asher Preska Steinberg
  2. Mingzhi Lin
  3. Edo Kussell
(2022)
Core genes can have higher recombination rates than accessory genes within global microbial populations
eLife 11:e78533.
https://doi.org/10.7554/eLife.78533

Share this article

https://doi.org/10.7554/eLife.78533

Further reading

    1. Microbiology and Infectious Disease
    Brian G Vassallo, Noemie Scheidel ... Dennis H Kim
    Research Article

    The microbiota is a key determinant of the physiology and immunity of animal hosts. The factors governing the transmissibility of viruses between susceptible hosts are incompletely understood. Bacteria serve as food for Caenorhabditis elegans and represent an integral part of the natural environment of C. elegans. We determined the effects of bacteria isolated with C. elegans from its natural environment on the transmission of Orsay virus in C. elegans using quantitative virus transmission and host susceptibility assays. We observed that Ochrobactrum species promoted Orsay virus transmission, whereas Pseudomonas lurida MYb11 attenuated virus transmission relative to the standard laboratory bacterial food Escherichia coli OP50. We found that pathogenic Pseudomonas aeruginosa strains PA01 and PA14 further attenuated virus transmission. We determined that the amount of Orsay virus required to infect 50% of a C. elegans population on P. lurida MYb11 compared with Ochrobactrum vermis MYb71 was dramatically increased, over three orders of magnitude. Host susceptibility was attenuated even further in the presence of P. aeruginosa PA14. Genetic analysis of the determinants of P. aeruginosa required for attenuation of C. elegans susceptibility to Orsay virus infection revealed a role for regulators of quorum sensing. Our data suggest that distinct constituents of the C. elegans microbiota and potential pathogens can have widely divergent effects on Orsay virus transmission, such that associated bacteria can effectively determine host susceptibility versus resistance to viral infection. Our study provides quantitative evidence for a critical role for tripartite host-virus-bacteria interactions in determining the transmissibility of viruses among susceptible hosts.

    1. Biochemistry and Chemical Biology
    2. Microbiology and Infectious Disease
    Carlo Giannangelo, Matthew P Challis ... Darren J Creek
    Research Article

    New antimalarial drug candidates that act via novel mechanisms are urgently needed to combat malaria drug resistance. Here, we describe the multi-omic chemical validation of Plasmodium M1 alanyl metalloaminopeptidase as an attractive drug target using the selective inhibitor, MIPS2673. MIPS2673 demonstrated potent inhibition of recombinant Plasmodium falciparum (PfA-M1) and Plasmodium vivax (PvA-M1) M1 metalloaminopeptidases, with selectivity over other Plasmodium and human aminopeptidases, and displayed excellent in vitro antimalarial activity with no significant host cytotoxicity. Orthogonal label-free chemoproteomic methods based on thermal stability and limited proteolysis of whole parasite lysates revealed that MIPS2673 solely targets PfA-M1 in parasites, with limited proteolysis also enabling estimation of the binding site on PfA-M1 to within ~5 Å of that determined by X-ray crystallography. Finally, functional investigation by untargeted metabolomics demonstrated that MIPS2673 inhibits the key role of PfA-M1 in haemoglobin digestion. Combined, our unbiased multi-omic target deconvolution methods confirmed the on-target activity of MIPS2673, and validated selective inhibition of M1 alanyl metalloaminopeptidase as a promising antimalarial strategy.