Contingency, repeatability, and predictability in the evolution of a prokaryotic pangenome

E. coli genomes were downloaded from the National Center for Biotechnology Information (NCBI) genome database (28) using the NCBI command line utility “datasets” 12.17.2 accessed 01/05/2022. All annotated E. coli genomes were downloaded provided that they were of the highest completeness level (complete). If a genome had been assembled by both the GenBank and RefSeq methods, the GenBank (GCA) annotation was maintained; however, if only a RefSeq (GCF) annotation was in the database, this was retained. The final dataset consisted of 2,341 genomes. The full list of accession numbers is included in SI Appendix, Supplementary List 1.

Genomes were re-annotated with PROKKA version 1.14.6 (29), employing the “bacteria” annotation mode. We applied an E-value threshold of 10⁻⁹ and mandated a minimum query coverage of 80% for assigning functional annotations. The E. coli pangenome was inferred using Panaroo version 1.2.9 (30) with the mode set to sensitive. It was important to include rare genes, owing to the possibility that they would have important effects on the prediction of other genes.

The gene presence–absence matrix produced by Panaroo was processed so that genes present in more than 99%, or less than 1% of genomes were removed. The matrix was further modified by converting gene names to “1” and empty fields to “0” so that 1 represented presence and 0 represented absence. In addition, genomes with identical patterns of gene presence and absence were collapsed into the same vector. Genes with identical presence-absence patterns (PAPs) across our genome sample were also collapsed into gene family groups (SI Appendix, Table S1). Essentially, this meant that both the features and predicted variable in these analyses are PAPs rather than genes per se, most of which were represented by only one gene but some of which were represented by multiple genes. The list of genomes with nonunique gene repertoires is found in SI Appendix, Table S2.

To minimise the impact of phylogeny on understanding predictability, repeatability, and contingency, we impose the requirement that we only study genes whose distribution is not “clumped” on one or a few clades. A backbone phylogeny was inferred in order to evaluate the distribution of gene content across the tree. Alignments of universal single-copy genes were constructed using MAFFT version 7.490 (31). The resulting 337 gene alignments were concatenated to form a superalignment and a maximum likelihood phylogeny was reconstructed using IQTree 2.2.0 (32) with extended model selection (33) and 1,000 ultrafast bootstrap replicates (34). The tree was rooted at the midpoint. To assess the distribution of each gene on this tree, we calculated Fitz and Purvis’ D statistic (35), retaining only those predicted genes with a D statistic greater than zero, though source nodes could still have D < 0. To contextualise the meaning of the D statistic, we also calculated the parsimony score for each gene using PAUP v4.0a (36). This metric represents the minimum number of times a gene has shifted from present to absent or absent to present along the phylogeny. Filtration of the edges was carried out using Sqlite3 (37). In practice, we only study genes that have changed character state from present to absent or vice versa on the tree at least 8 times and where D is >0.

A model of gene presence or absence was inferred using Random Forests. Prediction models for each PAP were calculated separately, meaning that the number of times that our Random Forests were trained was equal to the number of PAPs in the processed matrix. For each PAP being predicted, the genomes were randomly split into a training set, equal to 75% of the genomes and a test set of the remaining 25%. These sets were stratified by the gene presence–absence state being predicted, meaning that the proportions of genomes with any given gene being present or absent remained approximately the same in the training and the test datasets. These test and training datasets were assigned independently during the prediction of every PAP in every analysis. Decision trees were then generated by taking a random sample of genes with size equal to the square root of the number of genes in the dataset and the gene that most evenly split the training set out of this subset formed the first node in the tree (according to ref. 38). This process was repeated on each side of each decision node until either the maximum depth was reached or the remaining sample of genomes all had the same state for the gene in question. The impact of varying the number of trees produced and their maximum depth was empirically assessed based on the average effectiveness of the analysis in predicting the presence or absence of each gene (SI Appendix, Fig. S1). We chose the maximum depth and number of trees above which the performance of the models on the test set did not increase substantially. For the analyses in this manuscript, including those where the dataset was downsampled, we used 1,000 trees with a maximum depth of 16 nodes per tree.

For each gene (or set of genes with the same PAP), a prediction of either its presence or its absence was obtained for each genome in the test set according to the model generated using the training set. For each gene, four performance metrics were taken. These were Accuracy, Precision, Recall, and F1 score, all defined according to ref. 39. As our classification algorithm has two classes, we recorded a version of recall, precision, and F1 score separately for both the gene presence and for the gene absence classes.

In addition to performance statistics, the Gini importance (40) of each gene in predicting each other gene was added to an n by n matrix where n is the number of PAPs in the dataset. Here, the Gini importance is the contribution of a predictor PAP in separating the dataset into genomes where the test gene is present and those where it is absent averaged over all trees in the Random Forest. This Gini importance value can be used as a measure of the strength of the impact of the predictor gene on the presence/absence state being predicted. This statistic is appropriate when the predictor variables have two classes, as is the case here (41). All machine learning algorithms were implemented using the scikit-learn Python module version 1.0.1 (38). Code is available at https://github.com/alanbeavan/pangenome_rf (SI Appendix, Fig. S2).

In principle, it is possible that a gene’s presence may be predictable purely due to coincidental correlations arising by random transfer and loss across the phylogeny. The number of predicted genes occurring in this way can be thought of as a baseline level of false positives. We simulated the presence or absence of genes on each tip of the tree according to an nonreversible All rates Different model which allows for a different rate of 1-to-0 and 0-to-1 transitions (see ref. 42). For each gene family, the substitution matrix was derived from the gene presence and absence patterns observed across the phylogeny. This was achieved using the “fitDiscrete” function from the “geiger” R package (43). Following this, a presence and absence pattern for each gene was simulated using the “simulate_mk_model” function, which is part of the castor R library (44). For simulation, the probability of the gene being present at the root of the tree was equal to the root probabilities inferred by IQTree (32). We employed ancestral state reconstruction and extended model selection and used the phylogeny inferred in this study as a fixed topology. Simulations were conducted 100 times and then analyzed in a manner similar to the empirical data. However, for computational efficiency, only those Presence-Absence Patterns (PAPs) with a high F1 score (as determined in the post-processing of the network) had their D statistic estimated.

A mechanistic explanation for prokaryotic pangenome origin and evolution is emerging (64, 65). Here, we specifically focussed on detecting pangenome-wide repeated and predictable patterns of evolution using a Random Forest approach. In effect, we have asked whether gene content evolution is predictable and whether within-species evolution is constrained by intragenomic forces. The benefit of the Random Forest approach is twofold. First, it allows us to test our models of gene presence and absence on an independent test set that is not used to generate the model. This means that genes can be classified not only according to the genes that influence their likelihood of occurrence but by whether their prediction is generalisable. This is a study that shows that gene presence or absence is predictable based only on other genes in the genome. Second, Random Forests can consider complex relationships as well as simple pairwise correlations. For example, a hypothetical gene A may predict the presence of gene B only in the absence of gene C. These two aspects of the method differentiate it from previous methods (17–22).

Though it was necessary to filter the dataset to remove rare and almost universal gene families (almost 50% of the accessory genome), we identified strong predictors for approximately 30% of the remaining gene families, roughly half of which could be explained by correlation with the phylogeny, leaving 16.7% of pangenome as both predictable and distributed widely across the species. While this leaves much of the accessory genome in a category of nonpredictable, given the current data and method of analysis, it must be viewed in the context that we have only used a single species pangenome and therefore might ask how the inclusion of additional species and datasets might lend additional predictive power. Given that downsampling our dataset reduced the proportion of gene families that could be well predicted, it is likely that the addition of more genomes would aid prediction, and the effect of broadening the taxon sampling would be of interest in future studies. The pattern of repetition and predictability that we observe across more than 2,044 gene families is compatible with a model of deterministic evolution and more difficult to reconcile with an evolutionary process dominated by contingent events. This compares with those gene families that we were not able to predict accurately for which we cannot rule out contingency on outside factors driving their evolution. Furthermore, given the way in which we have analysed the data, the reasons for observing such widespread predictability stem from intragenomic natural selection causing biases in the cohort of other genes that are acquired, retained, or avoided by a genome. Another possibility is that two (or more) genes are found together because they are both selectively beneficial in the same environment. Hence, when a lineage colonises this environment, both genes are selectively recruited independently. The reasons for co-occurrence and avoidance are largely speculative and would require complex experiments to decipher. Biased presence–absence patterns have been noted previously (19, 21, 22, 66).

Throughout this study, we have been using the heuristic that homology is closely related to functional similarity and that the effects that genes have on one another are due in some way to the encoded functions of the proteins. Furthermore, we assume that these functions remain constant over time, despite evolutionary changes in the gene sequences. This is of course unlikely to be completely accurate (67, 68). Our dataset almost certainly contains genes that are placed into the same family but have different functions and confer different fitness effects. This limitation would certainly weaken our ability to make accurate predictions, though only a small portion of the dataset has verified function. Indeed, mutation can change the function of genes without changing the gene family to which they are assigned (69, 70). It is outside the scope of this paper to assess the role of point mutation on the predictability of bacterial genome evolution, but we hope that this approach may be applied to aspects of genome evolution outside of HGT.

We note that a lot of gene families are not well predicted given the current dataset and method of analysis. Our dataset might be too small to supply enough power for statistical inference, or indeed, intragenomic fitness effects do not overcome random genetic drift. Downsampling the dataset produced a significant reduction in prediction power. While the outcome with significantly larger datasets is uncertain, the link between dataset size and predictive power suggests that increasing the dataset might enhance prediction accuracy. Future inclusion of other factors such as external environment, gene expression levels, protein interaction, phenotypic, or modification data may also aid prediction. Recently, the development of gene-specific evolutionary rate models has shown a significant level of metabolic predictability (71). Of course, a gene whose presence or absence is nonpredictable in the current analysis might continue to be nonpredictable in enormous datasets, precisely because it is not impacted by co-occurring genes.

We have taken great care to minimise the confounding effect of genome relatedness. By applying the D-score filter (35), all gene families in our predictive model have a history of being gained or lost at least eight times across the pangenome, and furthermore, the distribution of any gene family cannot be “clumped” or restricted to just one part of the backbone tree. This is not a perfect way to eliminate the effect of phylogeny on the associations, but by coupling this approach with a very high threshold for predictability, we end up with an approach that shows correlations that are not just because a pair of genes happen to be in the same clade. Furthermore, by simulating gene presence–absence on the tree we inferred, we observe that chance relationships can only account for ~1.5% of genes’ predictability. Stricter D score cutoffs would reduce this false discovery rate but at the expense of failing to identify some genes that are truly predictable. Similarly, decreasing the cutoff would identify a larger number of truly predictable genes but is likely to raise the false discovery rate.

We arbitrarily chose a conservative Gini importance cutoff of 0.01 though there is no standard procedure for choosing a cutoff. Gini importance is a measure of the reduction in ambiguity of the test variable (gene presence or absence), with each node that comprises the prediction variable (a predictor gene), averaged over the trees in the forest. A Gini importance of 0.01 means that the ambiguity of the predicted gene is reduced by an average of 1% in each tree (bearing in mind it will only be sampled in a subset of trees). Gini filtering limited the number of associations to the strongest 33,426 connecting 4,067 nodes, after filtering before filtering by model performance and D-statistic. A weaker cutoff value of 0.005, for example, would have included 70,581 edges connecting 6,830 nodes. The distribution of edge strengths suggests that the number of edges will increase exponentially as that threshold is reduced linearly (SI Appendix, Fig. S3). Similarly, we chose an arbitrary, conservative, cutoff of 0.9 for accuracy and F1 score for both classes. Applying an accuracy threshold of 80% and an F1 score threshold of 0.8 would have resulted in the inclusion of more than 70% additional genes and relationships. Under these criteria, the number of predictable genes with a D statistic greater than zero would have nearly doubled, increasing from 2,044 to 3,704.

Linkage clearly plays an important role in coordinating the cooccurrence of sets of genes, but it cannot explain all the results. Also, it is not clear whether linkage is the cause or a consequence of co-occurrence. According to the selfish operon theory (63, 72), we would expect two genes that provide a fitness benefit when found together in a genome to evolve to be physically closer together on the chromosome, so they are less likely to be separated via recombination events. The most likely explanation is that close linkage is both a cause and effect of being found together. However, in addition to gene family co-occurrence where linkage plays a role, there are thousands of examples of gene families cooccurring in a repeated manner where linkage clearly has not played a role throughout the evolution of the pangenome, either because the genes are separated by a long stretch of DNA or are on different genetic elements.

Genes avoiding being in the same genome is an established phenomenon, with Bruns et al. (73), for example, showing that within the genus Salinospora biosynthetic gene clusters encoding iron siderophores avoid one another. In that case, the clusters encoded iron chelators with very different structures, but near-identical chemical properties, and avoidance most likely stems from functional redundancy. In Fig. 2, we outline two situations where avoidance is apparent. In PAPs J-M, both sets of genes encode proton antiporters. It is not clear whether the two sets of genes can functionally complement one another, in which case avoidance might be caused by redundancy in the same way as the iron chelators in Salinosopra. It is also possible that the functioning of the two antiporters either results in toxic effects when present together or they are competing for the same cellular resources, which would obviously result in a reduction in fitness for the organism that possessed both sets of genes, compared with a close relative that contained only one set of genes. It is also possible that one of the sets of these genes is preferable to the other in some environments, but in others, the reverse is true. For example, multidrug resistance can be conferred by the presence of several sets of genes, none of which are functionally related but each of which is selected by the presence of different antimicrobial drugs used in a treatment. Indeed, ecology has been shown to be a major driver of HGT (74). The situation in Fig. 2 PAPs A-I is more enigmatic, with no obvious functional similarity between pac and symE. However, they are found separately in genomes far more often than expected, and furthermore, the appearance of one of these genes in a genome occurs simultaneously with the removal of the other and vice versa at least 30 times across the phylogeny (Fig. 2 PAPs E-F).

Using nomenclature and analogies from ecology, the diversity of motifs embedded in the Random Forest network analysis allows us to identify putative mutualistic, commensal, and competitive classes of relationship, (see ref. 9). In this sense, the pangenome exists as a broad ecosystem, with individual genomes acting as evolving localities, where genes can either potentiate or alternatively reduce the likelihood of the presence of another gene. Just like a macroecological system, we see groups of genes with strong reciprocal co-occurrence patterns which we call putative mutualisms. Sometimes, these are just pairs of genes, but they are often groups of more genes which all co-occur reciprocally. Focussing on putatively commensal relationships, we identify specific examples where one or more genes appear to make a genome more hospitable to another or several other genes. In these cases, the more abundant gene is highly likely to arrive in the genome first or simultaneously with the less abundant gene. We also see many cases where the arrival of a gene into a lineage is concomitant with the loss of another gene, and this pattern is repeated across the phylogeny. There are several possible underlying reasons causing genes to predict the presence or absence of other genes, including functioning in a common pathway or process, redundancy, physical linkage, and shared evolutionary advantage in a new environment.

Ecosystems are known to be dynamic; they also tend to be resilient and to be somewhat resistant to overall change (75). In our analysis of the E. coli pangenome, we see features that are consistent with this perspective. We see a very dynamic system of gene gains and losses; we see repetitive gains of the same cohorts of genes, and we see the establishment of sets of relationships that are persistent through time and across the phylogeny. Due to the diversity of the E. coli pangenome, each time a gene is recruited to a new genome, it finds itself in a different, and sometimes substantially different, genetic background. Nonetheless, we observe repeated, predictable patterns of evolution following a gene’s transfer.

With respect to Gould’s “replaying evolution’s tape” thought experiment, our results lead us to suggest that it is likely that rewinding the tape back to the start of E. coli evolution would still result in hundreds or thousands of predictable events taking place that are not contingent on those highly unlikely events unique to each replaying of the tape. It is doubtful that the exact same evolutionary trajectories would play out, but several motifs would emerge over time.

Other machine learning algorithms such as neural networks may also be able to improve predictions by finding more abstract or subtle patterns in the data. Employing these types of methodologies, combined with the identification of further predictor variables and their application to diverse, large datasets of bacteria, archaea, or eukaryotes, represents promising future directions for enhancing our understanding of pangenome evolution.

We would like to thank Mike McDonald for comments and suggestions on this article. We gratefully acknowledge the financial support provided by the Biotechnology and Biological Sciences Research Council for this research.