Pan-Genomic Architecture and Evolutionary Profiling of Buchnera aphidicola

Abstract

This repository contains the computational framework and execution protocols for a high-resolution, genus-level pan-genomic analysis of Buchnera aphidicola, an obligate intracellular endosymbiont of aphids (Hemiptera: Aphididae). The genus Buchnera represents one of the most ancient and deeply co-evolved mutualistic symbioses in the animal kingdom, with an estimated origin dating to approximately 160--280 million years ago. As a consequence of this sustained obligate intracellular lifecycle, B. aphidicola genomes have undergone extreme reductive evolution, exhibiting severe gene loss, pronounced AT-bias, and the near-complete erosion of mobile genetic elements.

Critically, the individual Buchnera lineages harbored by different aphid host species have been reproductively isolated for tens to hundreds of millions of years, resulting in inter-lineage ANI values that frequently fall below the canonical 95% species boundary. Rather than restricting this analysis to a single host-associated lineage, this study deliberately encompasses all available lineages to capture the full scope of pan-genomic decay across the genus. This cross-host, multi-lineage design is essential for understanding the divergent trajectories of reductive evolution that characterize obligate endosymbiosis.

Due to this deep evolutionary divergence and the highly degraded nature of endosymbiont coding sequences, standard pan-genomic clustering algorithms that rely on rigid sequence identity thresholds (e.g., CD-HIT alone at 95% identity) are fundamentally inadequate for resolving orthologous relationships across Buchnera lineages. To overcome this limitation, this study employs the Pan-Genome Analysis Pipeline 2 (PGAP2), which leverages Fine-Grained Feature Networks and synteny-aware graph-based clustering. This approach enables the accurate identification of highly divergent orthologs that would otherwise be mis-classified as strain-specific singletons by conventional methods.

The pipeline encompasses automated data acquisition from NCBI RefSeq, ANI-based quality diagnostics, MMseqs2-based feature clustering, Markov Cluster Algorithm (MCL) community detection, and statistical rarefaction modeling to determine pan-genomic trajectory.

---

Table of Contents

1. Repository Architecture
2. Prerequisites and Computational Environment
3. Genomic Data Acquisition
4. Methodology and Automated Execution
5. Output Artifacts and Biological Inference
6. Results
7. Interpretation and Biological Significance
8. References

---

Repository Architecture

To ensure absolute data provenance and computational reproducibility, the repository is structured chronologically, isolating raw inputs from executable code and all algorithmic derivations. Each directory contains its own README.md that documents the purpose and contents of that particular stage.

Buchnera_PGAP2_PanGenome/
|
|-- 00_raw_genomes/              # Curated GenBank Flat Files (.gbff) retrieved from NCBI RefSeq.
|
|-- 01_scripts/                  # Automation scripts for data retrieval and pipeline execution.
|   |-- 01_raw_genomes.sh        # Downloads and decompresses the input genomic assemblies.
|   +-- 02_run_pgap2_pipeline.sh # Orchestrates all three PGAP2 analytical phases.
|
|-- 02_pgap2_prep/               # Preprocessing outputs: ANI matrices, quality control, and outlier flagging.
|   |-- genome_index/            # Indexed genome references for pairwise alignment.
|   |-- refined_input/           # Symlinks or filtered inputs post-QC.
|   |-- clust_0.98/              # CD-HIT clustering at 98% identity threshold.
|   |-- clust_0.99/              # CD-HIT clustering at 99% identity threshold.
|   +-- tmpvhzn2fuj/             # Temporary working directory (auto-generated by PGAP2).
|
|-- 03_pgap2_main/               # Primary pan-genome map: PAV matrix, gene network, and summary statistics.
|   |-- clust_0.98/              # MMseqs2 clustering at 98% identity threshold.
|   |-- clust_0.99/              # MMseqs2 clustering at 99% identity threshold.
|   +-- tmp5vne8jnj/             # Temporary working directory (auto-generated by PGAP2).
|
+-- 04_pgap2_post/               # Post-processing: rarefaction curves, statistical profiling, and vector graphics.

---

Prerequisites and Computational Environment

The execution of this pipeline requires a dedicated, isolated computational environment to manage the complex dependency chain, specifically the integration of the MMseqs2 clustering engine, the MCL graph clustering library, the DIAMOND aligner, and the R-based visualization suite.

Environment Instantiation

The environment is provisioned via the mamba package manager (installed within a base conda environment) to leverage rapid dependency resolution:

conda create -n pgap2_env -c conda-forge mamba -y
conda activate pgap2_env
mamba install -c bioconda pgap2 -y

Note on Dependency Pinning

The pgap2 bioconda distribution strictly pins the pandas library to version 1.5.0, which utilizes the DataFrame.map() method. In pandas v2.x, this method was renamed to DataFrame.applymap(). Rather than forcing a potentially destabilizing package upgrade via the solver, a dynamic syntax patch is applied at runtime within the orchestration script (see Methodology, Phase 2.5). This approach preserves the integrity of the entire dependency graph while ensuring forward compatibility with the statistical profiling module.

Core Software Dependencies

Software	Purpose
PGAP2	Pan-genome analysis pipeline orchestration
MMseqs2	Sensitive protein sequence clustering
DIAMOND	Fast protein alignment for ortholog detection
CD-HIT	Initial redundancy reduction
MCL	Markov Cluster Algorithm for graph partitioning
R (+ ggplot2)	Statistical visualization and PDF generation

---

Genomic Data Acquisition

To mitigate the annotation fragmentation frequently associated with de novo gene prediction algorithms (e.g., Prokka, Prodigal) on highly degraded endosymbiont genomes, this pipeline ingests natively annotated GenBank Flat Files (.gbff) directly from the NCBI RefSeq database. These files contain the authoritative, curated structural and functional annotations deposited by the original submitting groups, ensuring maximal annotation fidelity.

Five representative B. aphidicola assemblies were retrieved, spanning multiple aphid host species to capture the breadth of inter-lineage genomic variation across the genus:

Assembly Accession	Host Aphid Species
GCF_000009605.1 (ASM960v1)	Acyrthosiphon pisum
GCF_000010065.1 (ASM1006v1)	Cinara cedri
GCF_000007725.1 (ASM772v1)	Baizongia pistaciae
GCF_000007365.1 (ASM736v1)	Schizaphis graminum
GCF_000225445.1 (ASM22544v1)	Acyrthosiphon pisum (5A)

Data Retrieval Protocol

The acquisition is automated by the 01_scripts/01_raw_genomes.sh script:

cd 00_raw_genomes

wget -q --show-progress https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/009/605/GCF_000009605.1_ASM960v1/GCF_000009605.1_ASM960v1_genomic.gbff.gz
wget -q --show-progress https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/010/065/GCF_000010065.1_ASM1006v1/GCF_000010065.1_ASM1006v1_genomic.gbff.gz
wget -q --show-progress https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/007/725/GCF_000007725.1_ASM772v1/GCF_000007725.1_ASM772v1_genomic.gbff.gz
wget -q --show-progress https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/007/365/GCF_000007365.1_ASM736v1/GCF_000007365.1_ASM736v1_genomic.gbff.gz
wget -q --show-progress https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/000/225/445/GCF_000225445.1_ASM22544v1/GCF_000225445.1_ASM22544v1_genomic.gbff.gz

gunzip -f *.gbff.gz

---

Methodology and Automated Execution

The entire computational methodology is consolidated into a single orchestration script (01_scripts/02_run_pgap2_pipeline.sh). This design standardizes all algorithmic parameters and ensures deterministic data flow across the three analytical phases. The script terminates immediately upon any non-zero exit code (set -e) to preserve data integrity.

Phase 1: Preprocessing and Quality Control (pgap2 prep)

pgap2 prep -i ../00_raw_genomes/ -o ../02_pgap2_prep/

This module performs the following operations:

1. Genome Parsing: Loads all .gbff files from the input directory, extracts CDS features, and filters annotations by feature type.
2. Hierarchical Clustering: Constructs initial sequence clusters at two identity thresholds (98% and 99%) using CD-HIT for rapid redundancy reduction.
3. ANI Computation: Calculates the pairwise Average Nucleotide Identity matrix across all input strains. Strains exhibiting ANI values below the canonical 95% species boundary are algorithmically flagged as outliers. These flags serve as taxonomic diagnostics; in this study, all strains were intentionally retained for downstream analysis (see Methodological Rationale).
4. Quality Reporting: Generates interactive HTML and publication-quality PDF diagnostics (pgap2.preprocess.pdf) summarizing genome statistics, codon usage, and the ANI distance matrix.

Phase 2: Pan-Genome Map Construction (pgap2 main)

pgap2 main -i ../00_raw_genomes/ -o ../03_pgap2_main/

By design, the main clustering module receives the complete, unfiltered genome set from 00_raw_genomes/, bypassing the ANI-based outlier exclusions generated during Phase 1. This deliberate methodological decision elevates the analysis from a strict intra-species profile to a genus-level, cross-host pan-genome (see Methodological Rationale).

This is the core analytical module:

1. Protein Alignment: Deploys DIAMOND for sensitive all-vs-all protein alignment to construct the raw ortholog distance graph.
2. Graph Partitioning: Applies the Markov Cluster Algorithm (MCL) to resolve protein families from the distance graph.
3. Synteny-Aware Refinement: Constructs a Fine-Grained Feature Network incorporating gene synteny (conserved gene neighborhood) information to disambiguate paralogs from true orthologs, a critical step for organisms with extensive gene family contraction.
4. PAV Matrix Generation: Produces the comprehensive gene Presence/Absence Variation matrix (pgap2.partition.gene_content.csv), the foundational dataset for all downstream statistical analyses.

Phase 2.5: Dependency Fix (Pandas Syntax Patch)

PGAP2_PATH=$(python -c "import pgap2, os; print(os.path.dirname(pgap2.__file__))")
sed -i "s/pav = pav.map(/pav = pav.applymap(/" "$PGAP2_PATH/postprocess/profile.py"

This step dynamically resolves a known compatibility issue between the pandas v1.5.0 API pinned by PGAP2 and the syntax changes introduced in pandas v2.x. The deprecated .map() call on DataFrame objects is safely replaced with .applymap() in the localized PGAP2 source code, without modifying the upstream bioconda package.

Phase 3: Statistical Profiling (pgap2 post profile)

pgap2 post profile --pav ../03_pgap2_main/pgap2.partition.gene_content.csv --outdir ../04_pgap2_post/

This module derives the evolutionary trajectory of the pan-genome:

1. Rarefaction Analysis: Performs iterative combinatorial sampling to model the accumulation of core and pan gene clusters as a function of the number of genomes sampled.
2. Curve Fitting (Heaps' Law): Fits the pan-genome accumulation curve to Heaps' Law power-law model, n = k * N^alpha, where N is the number of genomes and alpha is the growth exponent. An alpha value less than 1.0 indicates an open pan-genome.
3. Core/Accessory Classification: Categorizes gene clusters into strict core, core, soft core, shell, and cloud compartments based on strain frequency thresholds.
4. Visualization: Generates interactive HTML reports and publication-quality vector PDFs for rarefaction curves and paralog statistics.

Execution

cd 01_scripts
chmod +x 02_run_pgap2_pipeline.sh
./02_run_pgap2_pipeline.sh

---

Output Artifacts and Biological Inference

Upon successful execution, the pipeline yields the following principal deliverables:

Preprocessing Diagnostics (02_pgap2_prep/)

File	Description
preprocess_results.html	Interactive HTML report with ANI matrix and QC
pgap2.preprocess.pdf	Vector-quality PDF of preprocessing diagnostics
preprocess.stat.tsv	Per-strain statistics (contig count, gene count, ANI, outlier status)
preprocess.gene_code.csv	Start/stop codon usage frequencies across strains

Pan-Genome Map (03_pgap2_main/)

File	Description
pgap2.partition.gene_content.csv	The PAV matrix (3,940 clusters x 5 strains)
pgap2.partition.gene_content.detail.tsv	Detailed cluster annotations with locus tags
pgap2.partition.gene_content.pav	Binary PAV matrix for downstream statistical tools
pgap2.partition.summary_statistics.txt	Core/shell/cloud classification summary
pgap2.partition.map.gml	Gene feature network in GML graph format

Statistical Profiles (04_pgap2_post/)

File	Description
postprocess_profile.html	Interactive rarefaction and profiling report
pgap2.postprocess_profile.pdf	Vector PDF of rarefaction curves and gene accumulation
pgap2.postprocess_stat_para.pdf	Vector PDF of paralog distribution statistics
postprocess.curve_fit.txt	Fitted mathematical models and pan-genome status
postprocess.rarefaction.tsv	Raw rarefaction sampling data
postprocess.pan_group_stat.tsv	Core/shell/cloud group counts and proportions

---

Results

Preprocessing: ANI-based Taxonomic Diagnostics

Of the five successfully loaded genomes, the ANI matrix revealed substantial inter-lineage divergence consistent with genus-level diversity. The preprocessing module correctly flagged three strains as falling below the canonical 95% ANI species boundary. However, as described in the Methodological Rationale, all five strains were intentionally retained for the primary clustering and statistical profiling phases.

Strain	ANI (%)	ANI Flag	Analysis Status
GCF_000225445.1 (ASM22544v1)	100.00	Representative	Retained
GCF_000010065.1 (ASM1006v1)	0.00	Below threshold	Retained
GCF_000009605.1 (ASM960v1)	86.08	Below threshold	Retained
GCF_000007365.1 (ASM736v1)	80.56	Below threshold	Retained
GCF_000007725.1 (ASM772v1)	77.35	Below threshold	Retained

Note: The anomalous ANI value of 0.00% for GCF_000010065.1 reflects the extreme genomic divergence of the B. aphidicola BCc lineage (Cinara cedri host), which possesses a dramatically expanded genome (~2.7 Mb, 2,657 CDS) relative to the canonical Buchnera AT-rich, ~640 kb genotype. This lineage is known to represent a fundamentally distinct evolutionary trajectory within the genus, having acquired a large complement of mobile genetic elements absent from all other Buchnera genomes.

The full preprocessing diagnostics, including the ANI distance heatmap and genome statistics, are available in 02_pgap2_prep/preprocess_results.html and 02_pgap2_prep/pgap2.preprocess.pdf.

Methodological Rationale: Deliberate Retention of Divergent Lineages

While the preprocessing module correctly identified three strains as taxonomic outliers based on the strict >95% ANI species boundary, these strains were intentionally retained in the primary clustering phase (pgap2 main). This methodological decision was made to elevate the analysis from a strict intra-species profile to a broader, host-divergent genus-level profile, capturing the severe pan-genomic decay that has accumulated across millions of years of isolated co-evolution with distinct aphid host environments.

This design is justified by the following biological considerations:

1. The 95% ANI threshold was developed for free-living bacteria with recombining populations. In obligate endosymbionts with strictly vertical transmission and no lateral gene transfer, lineage divergence progresses continuously without the homogenizing effects of recombination, rendering the species boundary biologically artificial.

2. The scientific question being addressed -- namely, the extent to which independent reductive evolution has erased the shared ancestral gene complement -- inherently requires the inclusion of maximally divergent lineages.

3. Restricting the analysis to the two strains that pass the 95% threshold would reduce the dataset to only two genomes, which is mathematically insufficient for meaningful rarefaction modeling or core/accessory classification.

Pan-Genome Composition

The feature network analysis resolved a total of 3,940 gene clusters across all five strains. The classification of these clusters into the standard pan-genomic compartments is as follows:

Compartment	Threshold	Cluster Count	Proportion
Strict Core	Present in 100% strains	1	0.03%
Core	99% <= strains < 100%	0	0.00%
Soft Core	95% <= strains < 99%	0	0.00%
Shell	15% <= strains < 95%	3,939	99.97%
Cloud	0% <= strains < 15%	0	0.00%

The overwhelming dominance of shell genes (99.97%) is a direct consequence of the genus-level analytical scope. Each Buchnera lineage, having been isolated within its respective aphid host for tens of millions of years, has undergone independent gene loss and pseudogenization. The result is a set of largely non-overlapping gene repertoires, with only a single universally conserved cluster persisting across all five lineages. This near-complete erosion of the core genome is among the most extreme pan-genomic signatures reported for any bacterial genus.

Rarefaction and Pan-Genome Trajectory

The statistical model fitted to Heaps' Law yielded the following parameters:

• Pan-genome fit: n = 998.791 * N^0.854
• Core genome fit: n = 6893.617 * exp(1.967 * N) + 15.185
• Pan-genome status: Open (alpha = 0.854 < 1.0)
• Coefficient k: 998.79 (initial rate of gene cluster discovery per genome added)

The open status indicates that each additional Buchnera lineage sampled from a new aphid host is expected to contribute a substantial number of previously unobserved gene clusters. This is consistent with the hypothesis that host-specific selective pressures drive lineage-specific gene retention, and that the full genetic repertoire of the genus has not been exhaustively sampled.

The full rarefaction curves and paralog distribution statistics are available in 04_pgap2_post/postprocess_profile.html, 04_pgap2_post/pgap2.postprocess_profile.pdf, and 04_pgap2_post/pgap2.postprocess_stat_para.pdf.

---

Interpretation and Biological Significance

The mathematical derivation of an open pan-genome (alpha approximately 0.85) across the Buchnera genus carries significant biological implications that merit careful interpretation within the context of obligate endosymbiosis and the genus-level scope of this analysis.

It is essential to note that the open pan-genome status does not imply active lateral gene transfer in the manner observed in free-living bacteria with genuinely open pan-genomes (e.g., Escherichia coli, Streptococcus pneumoniae). Buchnera has no known mechanism of lateral gene acquisition, no competence machinery, and no extant mobile genetic elements in most lineages. Rather, the mathematical openness of the genus-level pan-genome reflects three intertwined evolutionary processes:

1. Lineage-specific gene retention and loss. Different Buchnera lineages, isolated within their respective aphid host species for tens of millions of years, have undergone independent and largely irreversible trajectories of reductive evolution. The specific complement of ancestral genes retained in each lineage is therefore partially unique, generating the observed pattern of high shell gene content and an expanding pan-genome upon sequential genome addition. The near-total absence of a core genome (only 1 of 3,940 clusters is universally conserved) represents one of the most extreme cases of core-genome erosion documented in any bacterial genus.

2. Host-driven selective pressures. The nutritional requirements of different aphid host species exert distinct selective pressures on the Buchnera genome, favoring the retention of specific biosynthetic pathways (particularly essential amino acid synthesis) while permitting the pseudogenization and loss of others. For instance, the BCc lineage (Cinara cedri) has undergone a qualitatively different evolutionary trajectory, retaining a much larger genome with acquired mobile elements, reflecting the unique metabolic demands of its conifer-feeding aphid host.

3. Extreme AT-bias and sequence divergence. The severe mutational pressure toward AT content (approximately 75% AT in most Buchnera lineages) accelerates sequence divergence, making ortholog detection challenging even when functional conservation is maintained. The use of PGAP2's synteny-aware clustering was essential for resolving these highly divergent orthologs, which would have been missed by identity-threshold-based approaches.

4. Implications for endosymbiont genomics. These findings demonstrate that the concept of a "pan-genome" can be productively extended to obligate endosymbionts, provided that the analytical tools are sufficiently sensitive to handle extreme sequence divergence. The genus-level scope adopted here reveals that Buchnera lineages are not merely "reduced versions of the same genome" but rather represent distinct evolutionary experiments in minimal cellular life, each shaped by the specific ecological niche defined by their aphid host.

5. The sole strict core gene and the minimal genome boundary. Despite the extreme genomic erosion observed across the genus, the PGAP2 feature network successfully isolated a single gene cluster perfectly conserved across all five highly divergent lineages: rpsL, encoding the 30S ribosomal protein S12 (cluster clust_1781). The retention of rpsL as the sole strict core gene elegantly illustrates the boundary conditions of obligate endosymbiosis. While Buchnera has discarded vast swaths of metabolic, regulatory, and structural genes -- relying instead on the host cytosolic environment -- it must strictly conserve the core translational machinery required to synthesize the essential amino acids demanded by the aphid host. RpsL is a critical component of the small ribosomal subunit responsible for ensuring translational accuracy during mRNA decoding, a function that cannot be complemented by any host-derived factor. Its survival against millions of years of relentless reductive mutational pressure, across lineages separated by ANI values as low as 77%, underscores its non-negotiable role in maintaining translational fidelity within the minimal cell. Notably, rpsL is one of the most universally conserved genes across all domains of life and is frequently employed as a deep-time phylogenetic marker, making its identification as the sole Buchnera core gene both biologically expected and methodologically validating.

These findings are consistent with the broader literature on endosymbiont pan-genomics and underscore the importance of employing synteny-aware, network-based approaches such as PGAP2 for organisms at the far end of the reductive evolution spectrum.

---

References

1. Bu, C., Zhang, H., Zhang, F. et al. (2025). PGAP2: A comprehensive toolkit for prokaryotic pan-genome analysis based on fine-grained feature networks. Nature Communications, 16, 9865. https://doi.org/10.1038/s41467-025-64846-5

2. Shigenobu, S., Watanabe, H., Hattori, M., Sakaki, Y., & Ishikawa, H. (2000). Genome sequence of the endocellular bacterial symbiont of aphids Buchnera sp. APS. Nature, 407(6800), 81--86. https://doi.org/10.1038/35024074

3. Moran, N. A., McCutcheon, J. P., & Nakabachi, A. (2008). Genomics and evolution of heritable bacterial symbionts. Annual Review of Genetics, 42, 165--190. https://doi.org/10.1146/annurev.genet.41.110306.130119

4. Perez-Brocal, V., Gil, R., Ramos, S., et al. (2006). A small microbial genome: the end of a long symbiotic relationship? Science, 314(5797), 312--313. https://doi.org/10.1126/science.1130441

5. Tettelin, H., Riley, D., Cattuto, C., & Medini, D. (2008). Comparative genomics: the bacterial pan-genome. Current Opinion in Microbiology, 11(5), 472--477. https://doi.org/10.1016/j.mib.2008.09.006

6. Steinegger, M., & Soding, J. (2017). MMseqs2 enables sensitive protein sequence searching for the analysis of massive data sets. Nature Biotechnology, 35(11), 1026--1028. https://doi.org/10.1038/nbt.3988