The enumeration of genetic variation within microbial species is crucial for understanding their evolution, adaptation, and pathogenicity. Pangenomes, which encompass the entire gene repertoire of a species, provide a powerful framework for studying this variation. This study delves into Two Pangenomes Compare construction methodologies, specifically highlighting the benefits of a subtype-aware approach over traditional genome-based methods. We demonstrate how incorporating multilocus sequence typing (MLST) into pangenome construction refines the estimation of pangenome openness and improves the accuracy of extrapolating pangenome size, offering a more nuanced perspective on species-wide genetic diversity.
Figure 1: Subtype-Balanced Heaps’ Law Estimates of Pangenome Openness Across 12 Microbial Pathogen Species. (a) Illustrates the total number of genomes analyzed for each species. (b) Shows MLST subtype diversity within each species’ genome collection, measured by unique MLSTs per genome and Shannon diversity. (c) Depicts an example Heaps’ Law fit for genome and MLST count versus running pangenome size. (d) Compares openness values estimated with and without MLST control. (e) Presents openness estimates after MLST control against phylogenetic class.
Enhancing Pangenome Openness Estimation with Subtype-Based Analysis
Pangenome openness, a measure of the propensity of new genomes to introduce novel genes, serves as a key indicator of gene-level diversity within a species. Heaps’ Law, originally applied in linguistics, has been adapted to pangenomics to estimate openness by analyzing the relationship between pangenome size and the number of genomes sequenced. The power law exponent derived from this relationship quantifies openness. However, a critical challenge arises from subtype bias within genome datasets. As revealed by MLST classification, some species exhibit a skewed distribution of subtypes, potentially leading to pangenome openness estimates that are subtype-specific rather than representative of the entire species. This bias can result in underestimations of overall species-wide openness and inaccurate predictions of pangenome size.
To address this limitation, we introduce a refined approach to pangenome analysis that leverages MLST subtypes. We compared two pangenomes compare construction methods for estimating openness using Heaps’ Law:
- Genome-based approach: The traditional method, where genomes are randomly shuffled to generate orderings for Heaps’ Law calculation.
- MLST-based approach: A novel method where, for each ordering, one genome is randomly selected per MLST subtype, ensuring subtype diversity in the analysis.
Our findings demonstrate that MLST-based openness estimates are generally higher than genome-based estimates across the majority of species studied (10 out of 12). Notably, the E. faecium case, known for strong subtype bias, showed a near doubling of the openness estimate with the MLST-based approach. This highlights the improved sensitivity of the subtype-aware method in capturing broader species diversity.
For validation, we assessed the accuracy of both methods in extrapolating pangenome size. Heaps’ Law fits were generated using the first half of genomes (or MLST-representative genomes) and evaluated against the second half. The MLST-based approach exhibited a lower median mean absolute error (MAE) in both the fitted and extrapolated regions for most species (11/12 and 9/12, respectively). This improved predictive power, despite using fewer data points for fitting in the MLST-based method, underscores its robustness and accuracy in pangenome size estimation. While P. aeruginosa and S. enterica showed slightly higher MAE with the MLST-based approach, these species also exhibited overall poorer fits to Heaps’ Law, suggesting a general limitation of the model in these specific cases rather than a weakness of the MLST-based method itself.
Overall, the MLST-based Heaps’ Law approach provides a more consistent and potentially more accurate estimation of pangenome size, especially for datasets with subtype bias. It offers a refined representation of species’ genetic diversity, revealing clustering of openness values according to phylogenetic classification. Gammaproteobacteria species displayed the highest openness, followed by Bacilli and Campylobacter species with intermediate openness, and E. faecium and N. gonorrhoeae exhibiting the lowest openness.
Figure 2: Division of the Campylobacter coli Pangenome into Unique, Accessory, and Core Genomes Based on Gene Frequency. (a) Illustrates the gene frequency distribution P(x) and cumulative distribution F(x), highlighting peaks for “unique” and “core” genes. (b) Shows log-log plots of frequency distributions at low and high frequencies with power function models. (c) Depicts the division of the pangenome into unique, accessory, and core genomes based on the cumulative distribution fit.
Power Function-Based Pangenome Partitioning by Gene Frequency
Beyond overall pangenome characteristics, understanding the distribution of individual gene frequencies provides deeper insights into genetic diversity. Analyzing gene frequencies (the number of genomes in which a gene is present) across pangenomes reveals a consistent pattern: a peak for rare genes and another for highly conserved core genes. The cumulative gene frequency distribution exhibits an inverse sigmoidal shape, suggesting a natural tripartite division of the pangenome:
- Unique genome: Comprising rare, poorly characterized genes.
- Core genome: Encompassing highly conserved, essential genes.
- Accessory genome: Representing uncommon genes and capturing the bulk of gene-level diversity.
While previous approaches have used static thresholds or exponential function fitting for pangenome division, we developed a novel method based on fitting gene frequency distributions to the sum of two power functions. This approach accurately models the observed power-law-like behavior at both low and high frequency extremes. The accessory genome is defined relative to the inflection point of the fitted cumulative distribution. This method achieved high accuracy, with MAE values less than 0.5% of total pangenome size for most species and R² values exceeding 0.99 for the majority. The resulting frequency thresholds for core and unique genes highlight the inherent asymmetry in gene frequency distributions.
This power function-based division yielded core genomes ranging from 4.3% to 34.2% of pangenome size and accessory genomes ranging from 5.0% to 38.1%. Core and accessory genome sizes were comparable, and larger genomes tended to have more open pangenomes. This partitioning approach enables focused analysis on specific gene sets – conserved core genes or diverse accessory genes – rather than the often less informative unique genes.
Figure 3: Functional Enrichment in Core and Accessory Genomes Across 12 Species. (a) Shows the number of genes in core and accessory genomes. (b) Compares core genome size, accessory genome size, and pangenome openness. (c) Depicts functional enrichments by COG category in core and accessory genomes, showing log2 odds ratios and species with significant enrichment.
Functional Specialization of Core and Accessory Genomes
To explore the functional roles associated with gene frequency, we conducted functional enrichment analysis using Clusters of Orthologous Groups (COG) categories and Gene Ontology (GO) terms. Core genomes consistently showed enrichment in metabolic COGs (energy production, amino acid, nucleotide, and coenzyme metabolism) and essential cellular processes (translation, ribosomal structure, post-translational modification). Accessory genomes exhibited weaker but notable enrichment in COGs related to intracellular trafficking and defense mechanisms. Unique genomes, due to their limited characterization, showed no significant functional enrichment beyond unknown function categories.
GO term analysis corroborated these findings, with core genomes highly enriched in ribosome and RNA processing related terms. In contrast, accessory and unique genomes lacked significant enrichment in specific GO terms, except for broad categories like “cellular process.” This functional differentiation suggests that core genomes are enriched in fundamental metabolic and housekeeping functions, while accessory genomes contribute to more niche-specific or adaptive roles.
Further investigation at the ortholog group (OG) level revealed 168 OGs conserved across all 12 core genomes, many of which are essential for growth. These core OGs were predominantly involved in translation/ribosomal functions and core metabolic pathways, including purine metabolism, chorismate biosynthesis, coenzyme A biosynthesis, and fatty acid biosynthesis.
Figure 4: Distribution and Function of Shared Genes in Core and Accessory Genomes. (a) Shows the number of shared genes versus observation frequency in core genomes. (b) Shows the same for accessory genomes. (c) Breaks down the functions of 168 genes shared across all 12 core genomes by COG category, detailing metabolic pathways.
Sequence Diversity in Core Genes: Functional and Positional Insights
Moving to sequence-level diversity, we analyzed the frequency of protein sequence variants within core genes, calculating allelic entropy as a measure of diversity for coding sequences and flanking intergenic regions (5’IG and 3’IG). Limited correlation was observed between coding sequence variability and intergenic sequence variability.
To identify genes with contrasting sequence diversity, we classified the top and bottom 5% of core genes by allelic entropy as “diverse” and “conserved,” respectively. Quantile regression was used to account for gene length effects on coding sequence diversity. Functional enrichment analysis revealed a consistent enrichment of COG J (translation, ribosomal structure and biogenesis) among the least sequence-diverse core genes. Other COG categories showed weak or inconsistent associations with sequence diversity. Intriguingly, MLST genes exhibited a range of sequence diversity, from highly conserved to highly diverse within core genomes.
Figure 5: Functional Enrichment in Core Genes Based on Sequence Diversity. (a) Illustrates the workflow for identifying genes with high or low sequence diversity using coding, 5’IG, and 3’IG variant frequencies and entropy calculations. (b) Shows functional enrichment in genes classified as diverse or conserved by sequence entropy, highlighting COGs with positive mean log2 odds ratios.
Domain-Specific Variation in Conserved Core Genes
At the highest resolution, we examined positional sequence diversity within 76 highly conserved core genes, focusing on domain-level variation. Multiple sequence alignments (MSAs) were generated for protein sequence variants, and domain annotations were obtained. Domain entropy percentiles were calculated to assess domain variability relative to the entire gene.
Analysis of 443 gene-domain pairs identified 27 mutation-enriched domains and 61 mutation-depleted domains across species. Both enriched and depleted domains were functionally diverse, but a significant proportion of mutation-enriched domains were associated with aminoacyl-tRNA synthetases (AARSs). Further investigation of AARS domains revealed a function-dependent pattern of mutation enrichment. Domains involved in editing, anticodon binding, and tRNA binding tended to be mutation-enriched or neutral, while catalytic domains were predominantly mutation-depleted.
Figure 6: Mutation Enrichment in Protein Domains of Conserved Core Genes. (a) Depicts the workflow for computing mutation enrichment in protein domains using MSA entropy and domain entropy percentiles. (b) Shows species-specific mutation enrichment for 443 gene-domain pairs, highlighting significantly enriched or depleted domains. (c) Focuses on gene-domain pairs with significant multispecies mutation enrichment, with AARS-related domains labeled.
Figure 7: Mutation Enrichment in Aminoacyl-tRNA Synthetase Domains. (a) Shows the enrichment of AARS-related features among mutation-enriched or depleted domains. (b) Compares mutation enrichment in AARS domains to function. (c) Presents species-specific mutation enrichment for AARS-associated gene-domain pairs.
Conclusion
This study demonstrates the utility of a subtype-aware, MLST-based approach for constructing and analyzing pangenomes. By comparing two pangenomes compare methodologies, we highlighted the enhanced accuracy of the MLST-based method in estimating pangenome openness and extrapolating pangenome size, particularly in datasets with subtype biases. Furthermore, our analysis of gene frequency, functional enrichment, and sequence diversity provides a multi-layered understanding of genetic variation within microbial pangenomes. The domain-level analysis reveals positional constraints on sequence diversity, especially within essential genes like aminoacyl-tRNA synthetases. These findings underscore the complex interplay of gene conservation, functional specialization, and sequence evolution in shaping microbial pangenomes and offer valuable insights into microbial adaptation and diversity.
