Genome entropy and network centrality contrast exploration and exploitation in evolution of foodborne pathogens

A total of 3939 STM isolates collected between 2016 and 2021 and sequenced at the Institute of Clinical Pathology and Medical Research (ICPMR), NSW Health Pathology were included in the analysis. Sequencing was performed as described in a previous study [18]. Demultiplexed sequencing reads with ${\gt}1\times 10^7$ reads per isolate were trimmed [19], based on a minimum quality read score of 20, and then assembled de novo using SPAdes (version 3.13.0 [20]). The quality of assemblies was assessed with Quast (version 5.0.2 [21]); only assemblies with ${\lt}200$ contigs and N50 values of ${\gt}$50 000 bp were included in further analysis. Contigs were annotated with predicted genes using Prokka (version 1.13.3 [22]). The bioinformatic tool, Roary (the pangenome pipeline; version 3.12.0 [23], with default parameters) was used to partition the gene content of all isolates. The software defines the core genes as genes present in 99%–100% of isolates, soft core (in 95%–99% of isolates), shell (in 15%–95% of isolates) and cloud gene content (in 0%–15% of isolates). Single nucleotide polymorphisms (SNPs) were identified from the 3792 383 base pairs defined as the core genome, using SNP-sites [24] without recombination masking [25]. The SNP distance G_ij between isolates i and j is defined as the number of core SNPs by which they differ. The median value of G_ij across all 3939 STM isolates was 383 (range 0–2051).

We first derived a complete undirected network in which nodes represent individual isolates and edges represent genetic distances. Specifically, following an earlier study [18], the edge weight between adjacent nodes i and j in the complete undirected network was set to the pairwise SNP distance, G_ij . We then used the derived undirected network to construct a 2-dimensional space in terms of two quantities measured for each node i: the closeness centrality l_i , and the prevalence defined as the number of close neighbours within genetic distance $G_{\textrm{max}} = 20$ from node i. The constructed centrality-prevalence space relates a structural property of the pathogen (its network centrality) to the functional property (its neighbourhood prevalence). The latter can be interpreted as the pathogen's fitness capturing the average spread and virulence across its genetic neighbourhood [13]. Following prior studies [12, 13], we use the normalised node closeness centrality, computed as the inverse sum of the length of the shortest paths between the node of interest and other nodes, as described below.

2.2.1. Closeness centrality

Node closeness centrality captures the relative importance of node i to other nodes in the network by distance proximity. It is calculated as the inverse average shortest distance from node i to other nodes:

$$\begin{equation} l_i=\frac{N}{\sum_{j=1, i \neq j}^N d_{i,j}} \end{equation} \tag{ 1 }$$

where N is the number of nodes in the network, $d_{i,j}$ is the length of the shortest path between any pair of nodes i and j, defined as the sum of the edge weights $G_{k_1 k_m}$, where $k_1 = i$ and $k_m = j$, on the path with $(m-1)$ edges. If m = 2, then $d_{ij} = G_{ij}$.

If we take inverse of l_i , $\lambda_i = 1 / l_i$, and average λ_i across the entire network, the result yields the average path length L of the network:

$$\begin{align} L =\frac{1}{N(N-1)}\sum_{i=1}^{N}\sum_{j=1, i \neq j}^{N} d_{i,j} \end{align} \tag{ 2 }$$

$$\begin{align} \kern-4.4pt = \frac{1}{N-1} \sum_{i=1}^{N} \left(\frac{1}{N}\sum_{j=1, i \neq j}^{N}d_{i,j}\right) \end{align} \tag{ 3 }$$

$$\begin{align} & =\frac{1}{N-1} \sum_{i=1}^{N} \frac{1}{l_i} = \frac{1}{N-1} \sum_{i=1}^{N} \lambda_i. \end{align} \tag{ 4 }$$

In this study, we employ the average length of the shortest paths to measure how individual nodes contribute to the average genetic variation in the population. The relationship between the closeness centrality and the average path length suggests that the closeness centrality is a more suitable characteristic of genetic diversity of the population on average than alternative centrality measures, such as betweenness centrality [26] or eigenvector centrality [27].

2.2.2. Thresholded undirected network

Using the undirected network, we then constructed a thresholded undirected network (figure 1) where the edges are only inferred between genetically similar isolates, within the genetic distance threshold $G_{\textrm{max}} = 20$. When comparing Salmonella core genomes, distances of 0–5 SNPs are highly indicative of STM isolates collected from the same point source outbreak—if they are collected within a three month period—whereas distances of more than 20 SNPs suggest that isolates are epidemiologically unrelated [28–30]. The choice of $G_{\textrm{max}} = 20$ therefore produces edges between pairs of plausibly-related isolates while excluding edges where the number of core SNPs would imply the divergence of strains occurred well outside of the study period. In our previous study of Salmonella Enteritidis genotype networks [18], we showed that the network analysis was robust to choices of G_max between 10 and 30.

Zoom In Zoom Out Reset image size

We highlighted nodes by applying a centrality threshold of $3.3\,\times\,10^{-3}$: nodes with centrality below the threshold are coloured blue or cyan and nodes with a higher centrality are coloured orange or yellow, as shown in figure 1. This centrality threshold is determined by an observed partition of isolates in the centrality/prevalence space (see section 3.1 for more details). The centrality-prevalence space is used for identifying patterns and phases in the evolutionary dynamics of the pathogens.

We also derived a directed SNP sub-network that captures both genetic and temporal proximity, with a directed edge only inferred between isolates if they are (i) genetically close ($G_{\textrm{max}} = 20$), and (ii) collected within a fixed time window ($T_{\textrm{max}} = 30$ days), with the edge direction pointing from the predecessor to the successor. In cases where the isolates are collected on the same day, two directional edges are inferred. We chose the collection date as a proxy to the date of infection because foodborne salmonellosis is characterised by acute symptoms with the onset periods varying between 6 hours to 6 days. In addition, the time window, T_max, is an approximation of the median shedding period commonly associated with foodborne salmonellosis [31]. A visual representation of the directed SNP sub-network is shown in figure 2.

Zoom In Zoom Out Reset image size

A directed path is formed by connecting directed edges and may be considered as an approximation of pathogen adaptations over time. Previous studies of directed MLVA networks [13] suggested a distinction between paths leading to 'successful' or 'unsuccessful' adaptation, by associating these types with the increasing or decreasing connectivity along the path. In this work, we applied this distinction to the directed network derived using the WGS dataset. In order to trace genetic changes over time, we identified $38,698$ sufficiently long (m > 5) directed paths formed by connecting a sequence of directed edges. For each path we compared changes of the neighbourhood prevalence between the start and end nodes, distinguishing positive (successful adaptation) and negative (unsuccessful adaptation) changes. Some examples of successful and unsuccessful paths are shown in figure 3.

Zoom In Zoom Out Reset image size

These directed paths were projected onto the centrality-prevalence space for the complete undirected SNP network, and for each point (i.e. each node) the prevalence changes were averaged across all paths originating from this node. For each point in the centrality-prevalence space (i.e. for each node), we quantify the average changes in prevalence along both successful and unsuccessful paths originating from this node. The resultant estimated probability density of directed paths distinguishes between the pathogens which tend to develop towards a higher or lower neighbourhood prevalence [13]. Thus, the probability density captures the expected changes in prevalence during a subsequent pathogen adaptation (measured via the likelihood of either expanding or contracting genetic neighbourhood).

The pangenome comprises the entire set of $13,555$ genes from all $3,939$ isolates. After excluding the core and soft core genome (genes present in greater than 95% of isolates) and the cloud genome (genes present in fewer than 15% of isolates), the remaining shell genome consisted of 713 genes. The high-resolution WGS dataset enabled us to further explore the evolutionary pathways by looking at the genetic composition of different groups of isolates (e.g. isolates in the low-to-mid and mid-to-high centrality groups). We computed normalised entropy of the shell genome for each group.

The shell genome entropy is computed using gene presence/absence data over a fixed number of temporally ordered isolates (n = 10, note that collection date is not uniformly spaced). See figure 4 for a schematic of the shell genome entropy computation. The normalised Shannon entropy, S, is then calculated as follows [32]:

$$\begin{equation} S =\frac{H}{H_{\textrm{max}}}=\frac{-\sum^{Q}_{i=1} p_i \log p_i }{H_{\textrm{max}}} \end{equation} \tag{ 5 }$$

where Q is the number of genes, p_i is the probability of the gene i being present in the considered group of isolates, and $H_{\textrm{max}} = {\ln}(Q)$ the maximum entropy.

Zoom In Zoom Out Reset image size

In order to verify that the two proposed types of the evolutionary dynamics (exploration and exploitation) can be characterised as phases, separated by a properly defined phase transition, we considered the pathogen adaptation as a thermodynamic process developing in response to changes in the control parameter defined as the node closeness centrality.

We employed the Fisher information to capture the sensitivity of the shell genome probability distribution to changes in the closeness centrality. In general, Fisher information [33] measures the amount of information that an observable random variable $\mathcal{X}$ contains about an unknown parameter θ, given the probability function $p(x|\theta)$ over $x \in \mathcal{X}$ with respect to parameter θ:

$$\begin{align} F[\theta] & = \mathrm{E} \left[ \left. \left (\frac{\partial}{\partial \theta} \log p(x|\theta) \right) ^2 \right| \,\theta \right] \end{align} \tag{ 6 }$$

$$\begin{align} & = \int_{x} \left (\frac{\partial}{\partial \theta}\log p(x|\theta) \right)^2 p(x|\theta) \ \mathrm{d} x \end{align} \tag{ 7 }$$

where $\mathrm{E} [\cdot | \theta]$ denotes the conditional expectation with respect to θ.

There is a well-established association between phase transitions and the Fisher information [34–41]: the Fisher information diverges at critical points, being proportional to the rate of change of the corresponding order parameter [37]. In the context of evolutionary dynamics of pathogens, the system possibly undergoing a transition is a pathogen which may change its adaptation behaviour from explorative to exploitative in response to changes in its node centrality (i.e. control parameter). In general, the node centrality may be expected to change in response to the underlying evolutionary dynamics. Hence, in this study, we consider changes in the node centrality as a process corresponding to varying the control parameter. On the other hand, the prevalence may correspond to the order parameter characterising the outcome of the evolutionary process. The advantage of using Fisher information in identifying possible phase transitions is that the computation of the corresponding order parameter is not required—instead, this model-independent method seeks specific critical points in the centrality space at which the Fisher information peaks.

We computed normalised Fisher information using gene presence/absence data over the isolates ordered by centrality with a fixed centrality increment of $1 \times 10^{-4}$. Following [42–44], the normalised discrete Fisher information is calculated as follows:

$$\begin{equation} F[\theta]=F_0 \sum^{R-1}_{i=1} [\sqrt{p_{i+1}}-\sqrt{p_i}]^2 \end{equation} \tag{ 8 }$$

where R is the number of centrality increments (steps), and F₀ is a normalisation constant:

$$\begin{align} F_0= \left\{\!\!\! \begin{array}{l@{}l} 1 & \textrm{if} \ p_i^*\!=\!1 \textrm{ for } i^*\!=\!1 \textrm{ or } i^*\!=\!R, \textrm{ and } p_i \neq 0 \ \forall i \neq i^* \\ \frac{1}{2} & \textrm{otherwise.}\end{array}\right. \end{align} \tag{ 9 }$$

The constructed centrality-prevalence space (figure 5) shows a clear structure with at least two groups separated by their closeness centrality values: the isolates with low-to-mid centrality ($\leqslant 3.3 \times 10^{-3}$, shown in blue/cyan), and the isolates with mid-to-high centrality (${\gt}3.3 \times 10^-3$, shown in yellow/orange). We draw an additional distinction by separating the isolates by neighbourhood prevalence (set at 100 isolates marked by the dashed line), above which isolates are considered high prevalence.

Figure 5 shows that (i) isolates in the same centrality group (low-to-mid, shown in blue; or mid-to-high, shown in orange) typically belong to the same network component, and (ii) isolates with high neighbourhood prevalence (shown in yellow and cyan) usually cluster together within the network components with some peripheral nodes from the same centrality group.

Zoom In Zoom Out Reset image size

We also related the neighbourhood prevalence and the local clustering coefficient. This allowed us to investigate the role of clusters and hubs which are known to contribute to epidemic propagation [45]. However, the partition of isolates observed in the clustering-prevalence space (see figure A1(b) in appendix A) is not as structured as the one in the centrality-prevalence space.

Previous studies distinguished between the explorative and exploitative pathways in the STM genetic space and argued that a mid-to-high centrality region contains higher-risk STM pathogens [12, 13]. These studies relied on coarse-grained MLVA datasets, and these results have not been verified with high-resolution WGS data. In this work, in order to associate the two groups separated in the centrality-prevalence space (figure 5) with the exploration-exploitation distinction, we compared the shell genome entropy between the groups.

Figure 6(a) shows that the isolates with low-to-mid centrality tend to have a higher normalised shell genome entropy, thus demonstrating a greater genetic diversity, while the normalised shell genome entropy for the isolates with mid-to-high centrality tends to be smaller, showing a reduced genetic diversity (figure 6(b)). This difference confirms the exploration-exploitation distinction: the isolates with low-to-mid centrality tend to diversify, while the isolates with mid-to-high centrality are more stable in their shell genome. Figure B1 in appendix B illustrates the corresponding shell genome data.

Zoom In Zoom Out Reset image size

We also examined a possible transition between two phases (exploration and exploitation) in the centrality space by tracing the normalised Fisher information computed for all isolates ordered by their centrality values. Figure 7 shows that a small change in the centrality value from $2.6 \times 10^{-3}$ to $2.7 \times 10^{-3}$ results in an abrupt spike in the Fisher information. This indicates that the probability distribution of genes in the isolates around this centrality value is most sensitive to the centrality change. In thermodynamic terms, when the closeness centrality is considered as a control parameter, the observed spike in the Fisher information indicates that there is a phase transition between the lower and higher centrality regions. We interpret the resulting evolutionary dynamics as exploration and exploitation respectively. The critical regime is observed towards the end of exploration phase before the population self-organises into more exploitative dynamics.

Zoom In Zoom Out Reset image size

It is also worth noting that many isolates in the critical centrality range ($2.6 \times 10^{-3}$ to $2.7 \times 10^{-3}$) are associated with an outbreak cluster with MLVA profile 5-17-9-13-490. This is an MLVA profile with a substantially different genetic composition relative to the rest of the population. However, the Fisher information peak indicating a critical regime in the centrality space remains robust even when the outbreak isolates are excluded. This is shown in appendix D, by deriving the network for the truncated population (figure D3). We also found that the Fisher information peak is not sensitive to the removal of all isolates in the mid-centrality—high prevalence cluster (see appendix D). Robustness of the Fisher information peak indicates that the critical regime does not appear due to specific outbreaks emerging during the exploration phase. Instead it points to a phase transition in the overall evolutionary dynamics.

Interestingly, the pathogens within the medium centrality range exhibit features of both exploration and exploitation. This is revealed by comparing the centrality-partitions shown in figures 5, D1 and D2, with the partition becoming less sharp as the network is disturbed by the population truncation. Such mixing behaviours are often observed at critical regimes, with phase separation splitting a homogeneous mix into distinct phases (e.g. one homogeneous liquid separates into two different liquids [46]).

We also verified that the observed peak in Fisher information computed for the shell genome (using the gene presence/absence data) cannot be solely attributed to the population structure. In other words, there are no distinct sub-populations (defined in terms of SNP distance across the entire genome) which may have been partitioned by the critical centrality regime, as shown in appendix D (figure D4).

At the high-centrality region, there is a secondary peak observed when the centrality exceeds $4 \times 10^{-3}$ (figure 7). The second peak may indicate that the system is sensitive because it encountered another change in evolutionary dynamics, possibly an adaptation bottleneck (see the following section). We tested the robustness of the two observed peaks by using different centrality increments (steps) varying from a smaller ($5 \times 10^{-5}$) to a larger step ($2 \times 10^{-4}$), as shown in figure C1. While the peak near the critical value remains evident across different increments, the secondary peak in the high centrality range appears more sensitive to the increment value, possibly indicating its transient nature.

We then studied the conjecture that the isolates with mid-to-high centrality can be associated with the higher risk STM genotypes. Figure 8 shows two regions in the centrality-prevalence space: a dense 'expansion' region with a majority of nodes starting successful paths (shown in red), and a smaller 'bottleneck' region with nodes starting unsuccessful paths (shown in blue). Both regions are characterised by high centrality and high prevalence. It is also evident that isolates in the low-to-mid centrality region do not initiate sufficiently long paths, showing minimal neighbourhood prevalence changes.

Zoom In Zoom Out Reset image size

Earlier studies based on MLVA data covering a nine-year period reported that while the average prevalence increases, the centrality of adapting pathogens initially increases (from low to high values) and then decreases (from high to medium values) [13]. Our results show a partial agreement with this finding: the increasing neighbourhood prevalence is observed only in the group of isolates with increasing centrality, from medium to high values. Importantly, our study shows that the isolates with the highest centrality encountered a bottleneck (shown in blue in figure 8), a region forming an emergent 'barrier' for successful paths and slowing evolutionary dynamics. In other words, the paths originating within the bottleneck region did not yield increased neighbourhood prevalence, in contrast to previous studies [13] which reported that a minority of paths may continue through the bottleneck region and attain higher prevalence while decreasing centrality. In our case, the observed bottleneck is related to the decline in the total numbers of salmonellosis cases reported between 2015 and 2019, which has been attributed to the national Foodborne Illness Reduction Strategy [2, 47] which used to contain the STM outbreaks in NSW. The reduction in salmonellosis notifications may also be attributed to increased attention to hand and surface hygiene during the COVID-19 pandemic [2, 48]. The emergence of the bottleneck region explains the secondary peak in the Fisher information, shown in figure 7.