Plant Microbiome Ecotypes

1. Plant compartments impose a small but real functional shift on microbial communities (H1, weakly supported)

Canonical Phase 2b finding: compartment identity explains a small fraction of variance in functional profiles after location-vs-dispersion separation. db-RDA (constrained ordination on PCoA scores of Jaccard distance over the refined 17-marker panel, n=607 plant-compartment species) gives location-only R² = 0.060 (permutation p = 0.001, 999 permutations); PERMDISP confirms group-dispersion heterogeneity (F = 15.6, p = 0.001), with root tightest (mean centroid distance 0.452), phyllosphere intermediate (0.503), and rhizosphere most variable (0.528). About 84% of the total PERMANOVA R² is real centroid shift; 16% is dispersion. Root, rhizosphere, and phyllosphere are functionally distinct in centroid terms, just at a small effect size (~6% of variance) — see §11 "PERMANOVA + PERMDISP + db-RDA" for the decomposition.

Phase 1/2 chronology, retained for context: the original NB04 PERMANOVA on the Phase 1 25-marker panel reported R² = 0.527 (pseudo-F = 235.1, p = 0.001), but NB14's exclude-top-3-species sensitivity showed the residual R² is only 0.072 once a handful of dominant rhizobial / Pseudomonas clades are removed — i.e., most of the original 0.527 was a taxonomic-sampling artifact of a few genome-rich species rather than a community-wide signal. Among 1,136 plant-associated species classified across four compartments (root 292, rhizosphere 160, phyllosphere 157, endophyte 29), the per-marker Fisher tests in NB04 still showed strong directional patterns: root species had high ACC deaminase (OR = 69.3), T3SS (OR = 65.6), nitrogen fixation (OR = 14.5), and quorum sensing (OR = 24.1) enrichment in the original framework, and 69 of 96 marker × compartment pairs were significant. These per-function enrichments are not contradicted by the Phase 2b decomposition — only the headline magnitude of the multivariate R² was inflated.

(Notebooks: 01_genome_census.ipynb, 04_compartment_profiling.ipynb, 14_deferred_controls.ipynb, notebooks/_run_h1_dbrda.py. Verdict revised twice during Phase 2b: 2026-04-24 after the paired adversarial review and again 2026-04-25 after item-20 db-RDA.)

2. Beneficial genes are core-encoded; pathogenic genes are accessory (H2)

Beneficial (PGP) gene clusters are predominantly core genome (64.6% core fraction), significantly exceeding both the genome-wide baseline of 46.8% and pathogenic gene clusters (45.2% core fraction). The Mann-Whitney U test comparing per-species core fractions yielded U=83,567,419 (p=3.38e-125), with a bootstrap 95% CI for the beneficial-pathogenic difference of [0.089, 0.106].

Among specific functions, biofilm (83.3% core), IAA biosynthesis (78.1%), nitrogen fixation (72.3%), and phosphate solubilization (70.9%) had the highest core fractions. In contrast, T4SS (48.5% singleton), coronatine toxin (45.2% singleton), and effectors (39.8% singleton) were strongly accessory/singleton-enriched.

(Notebook: 05_genomic_architecture.ipynb)

3. Pathogenic gene clusters co-occur with transposases, suggesting HGT (H4 — partial)

Three HGT proxies were evaluated. Transposase/integrase co-occurrence was strongly positive: among 1,136 plant-associated species, those carrying singleton marker gene clusters were 16× more likely to also carry transposase/integrase singletons (Fisher's exact OR=15.95, p=8.8e-20). However, the overall singleton enrichment ratio for markers was 0.78 (markers are less mobile than the genomic average; Wilcoxon p<1e-300). Cross-species context variation (same OG core in one species, singleton in another) was higher for pathogenic markers (0.534) than beneficial (0.450), but not significantly so (Kruskal-Wallis p=0.122).

The mixed signal suggests that while individual HGT events are common (transposase association), plant-interaction markers as a class are under stronger purifying selection than the genomic average — likely because they are functionally important once acquired.

(Notebook: 05_genomic_architecture.ipynb)

4. Co-occurring genera show functional redundancy, not complementarity (H3 — not supported, small effect)

Canonical Phase 2b finding: co-occurring plant-associated genera in NMDC soil communities are slightly less complementary than random genus pairs (Cohen's d ≈ −0.4, permutation p < 0.001), supporting the redundancy interpretation but at a small effect size. The Phase 1 NB06 result of "Cohen's d = −7.54" was largely a Cohen's d formula error: NB06 divided the mean difference by the SD of permutation means (which shrinks toward zero by CLT, asymptotically), while the correct denominator is the raw pair-level SD. Re-computing NB06's max-aggregated complementarity with a consistent formula on the same data gives d ≈ −0.43; switching the aggregation from species-max to species-prevalence (the I1 fix in NB14) further reduces |d| to 0.39. So of the apparent "20× reduction" reported in the early Phase 2b draft, ~90% was the formula correction and ~10% was the aggregation correction. The direction is robust under both methods.

BacDive metabolite utilization cross-validated GapMind predictions at 83.1% consistency, confirming that the pathway completeness scores themselves are reliable. The corrected result is consistent with Louca et al. (2018), who argued that functional redundancy is an emergent property of microbial systems — and the small effect size is now within the range typical of comparable analyses, not pathologically large.

(Notebooks: 06_complementarity.ipynb, 14_deferred_controls.ipynb. See §11 "Prevalence-weighted GapMind complementarity" for the formula-vs-aggregation attribution.)

5. Fifty novel gene families distinguish plant-associated species (H5)

Genome-wide enrichment analysis of 5,671 eggNOG ortholog groups (OGs) identified 5,341 significantly associated with plant vs. non-plant species (q<0.05), of which 3,840 showed strong enrichment (OR>2). After phylogenetic control via phylum-level logistic regression, 50 novel OGs retained significance (all p<0.05, OR>1). The top hit was COG3569 (Fisher OR=8.92, phylo-controlled OR=6.01, q=7.7e-242), found in 54% of plant-associated species but only 12% of non-plant species.

The attenuation from Fisher OR to phylo-controlled OR (e.g., COG3569: 8.92→6.01) indicates partial phylogenetic confounding — plant-associated species cluster taxonomically — but a substantial ecological signal persists after correction.

(Notebook: 03_enrichment_analysis.ipynb)

6. Most plant-associated bacteria carry both PGP and pathogenic markers, but the dual-nature label is uninformative at species level — the continuous pathogen ratio is what discriminates

The majority of plant-associated species (65–85% per compartment) carry both PGP and pathogenic marker genes simultaneously. Two classification schemes were applied, producing different dual-nature rates:

• NB02 (binary marker presence): Any species carrying ≥1 PGP and ≥1 pathogenic marker is classified as dual-nature → 15,474/25,660 species (60.3%). This reflects the broad prevalence of general bacterial functions (flagella, chemotaxis, T6SS) that span both categories.
• NB07 (composite scoring): A weighted scoring framework (40% PGP markers, 20% core fraction, 15% complementarity, 10% metabolic breadth, 15% pathogen penalty) with adaptive thresholds (median of non-zero scores) produces a more stringent classification → neutral 38.0%, pathogenic 29.5%, dual-nature 25.2%, beneficial 7.3%.

The 35 percentage point difference between schemes reflects the distinction between carrying both marker types (nearly universal among bacteria due to ubiquitous functions like flagella and T6SS) and scoring highly on both axes when weighted by genomic architecture and metabolic context. The NB02 classification is more inclusive; the NB07 classification is more discriminating but also more assumption-laden.

Canonical Phase 2b validation: a curated 18-species panel of model organisms with experimentally confirmed phenotypes (7 known beneficial: P. simiae WCS417, P. protegens CHA0, B. velezensis FZB42, A. brasilense Sp245, R. leguminosarum, P. phytofirmans PsJN, S. meliloti; 7 known pathogenic: P. syringae, R. solanacearum, X. campestris, A. tumefaciens, E. amylovora, P. carotovorum, C. michiganensis; 4 neutral non-plant controls) was tested in NB13. All 14 of the beneficial + pathogenic ground-truth species were assigned to the dual-nature class, which means the categorical cohort label is uninformative at the species level — it cannot discriminate beneficial from pathogenic ground truth. The discriminative signal comes from a continuous metric: the pathogenicity ratio n_pathogen / (n_pgp + n_pathogen) is lower in known beneficial species (median 0.50) than in known pathogens (median 0.60), Mann-Whitney U = 9, p = 0.027 on N = 7 vs 7 — small effect at small sample size, but statistically consistent. The earlier Phase 1 figure of "92.7% genus-level agreement" was tautological because cohorts are defined from marker presence; the species-level Mann-Whitney is the honest replacement. Practical implication: the refined cohort labels (beneficial / pathogenic / dual-nature / neutral) should be treated as a coarse screen, and the continuous pathogenicity ratio is what discriminates within the dual-nature class.

Key dual-nature genera (Phase 1 / Phase 2 description, unchanged): Pseudomonas_E (358 species; PGP: HCN, phenazine, phosphate solubilization; Pathogenic: T2SS, T3SS, T6SS, effectors), Bradyrhizobium (73 species; PGP: nitrogen fixation, ACC deaminase; Pathogenic: T3SS, T4SS, CWDE), and Streptomyces (261 species; PGP: phenazine, phosphate solubilization; Pathogenic: CWDE, T3SS).

(Notebooks: 02_marker_gene_survey.ipynb, 07_cohort_synthesis.ipynb, 13_validation_pfam_subclade.ipynb. Validation framing revised 2026-04-24 after paired adversarial review.)

7. Novel OGs are functionally characterized: electron transport and iron-sulfur biosynthesis dominate (H5 extension)

All 50 plant-enriched OGs identified in Phase 1 were fully annotated via eggNOG descriptions, InterProScan domains, GO terms (48/50), and MetaCyc pathways (39/50). None were hypothetical proteins — every OG has a characterized function. The dominant functional themes are:

• Electron transport/energy metabolism (7/50, 14%): Cytochrome c oxidase subunits (COG1845/ctaE, OR=14.5, 93.2% plant vs 48.7% non-plant; COG0843/ctaD, OR=13.7), NADH dehydrogenase (COG0649), and ubiquinol-cytochrome c reductase (COG0723). The enrichment of high-affinity terminal oxidases is consistent with Preisig et al. (1996), who showed that cbb3-type cytochrome oxidases enable respiration in microaerobic root nodule environments.
• Iron-sulfur cluster biosynthesis (COG0316/iscA, 94.3% plant): Essential cofactors for nitrogenase, respiratory chain complexes, and oxidative stress defense. IscR-mediated Fe-S cluster regulation is linked to both pathogenicity and oxidative stress resistance in plant-associated bacteria (Fuangthong et al. 2015).
• Inorganic ion transport (6/50, 12%): Consistent with the nutrient-acquisition demands of rhizosphere colonization.
• Unknown function (19/50, 38%): Despite full annotation, these OGs have characterized domains but no clear plant-interaction mechanism, representing candidates for experimental characterization.

All 50 OGs had core gene fractions of 60.1–83.1% (baseline: 46.8%), reinforcing that plant-enriched functions are stably inherited in the core genome (consistent with Finding 2/H2). Cross-reference with MGnify rhizosphere KEGG modules confirmed that 8 modules from novel OGs are present in rhizosphere species.

(Notebook: 09_novel_og_annotation.ipynb)

8. Refined marker panel with KEGG module gating improves specificity but confirms persistent dual-nature (H0, H6)

Phase 1's 91-marker panel included ubiquitous bacterial functions (flagella, chemotaxis, T6SS, biofilm, quorum sensing, T2SS) that classified even Escherichia and Salmonella as 100% dual-nature. Phase 2 refined this to 17 plant-specific markers (9 PGP + 8 pathogenic) and applied KEGG module completeness gating for multi-gene systems:

• T3SS (M00332): Requiring ≥3 module genes reduced T3SS-positive species from 10,300 to 1,661 (−8,855 false positives, 86% reduction)
• Nitrogen fixation (M00175): Requiring ≥2 module genes reduced nif-positive species by 346 (14% reduction)
• T4SS (M00333): Minimal impact (most hits already multi-gene)

Despite the stricter panel, the dual-nature rate among plant-associated species increased from 60.3% to 78.7% (878/1,115 species). This counterintuitive result arises because the remaining markers — phosphate solubilization (57.9% prevalence), CWDE cellulase (63%), phenazine (63%), T4SS (45.8%) — are genuinely co-occurring in plant-associated genomes. The high dual-nature rate is biologically real, not an artifact of over-broad markers.

Host species extraction identified 1,307 species with plant host assignments across 19 testable hosts (≥10 species each). The most common hosts: rice (Oryza sativa, 286 species), Arabidopsis thaliana (184), wheat (122), maize (98).

Phylogenetic control: NB10's statsmodels logit with genus fixed effects failed to converge for 0/14 markers — the plant-association signal is heavily confounded with genus-level taxonomy and standard maximum-likelihood logistic regression cannot fit. Phase 2b closed this with two complementary tests (see §11 "Phylogenetic control"): (a) cluster-robust GLM with genus clusters (proper Wald inference, no L1 shrinkage) found 8/14 markers significant at q<0.05 BH-FDR with 95% Wald CIs excluding zero — ACC deaminase OR=8.01, T3SS OR=2.71, nitrogen fixation OR=2.71, phenazine OR=2.12, CWDE cellulase OR=1.87, phosphate solubilization OR=1.68, CWDE pectinase OR=1.59, effector OR=1.35; (b) within-genus label shuffling, the strictest species-level null, found 3/15 markers survive (nitrogen fixation, ACC deaminase, T3SS). The combined three-tier verdict is: 3 markers species-level robust, 5 cassette-level (genus-acquired but not species-distinguishing), 6 not robust under any test.

(Notebooks: 10_refined_markers_host_species.ipynb, 14_deferred_controls.ipynb, notebooks/_run_c1_cluster_robust.py. Three-tier framing finalized 2026-04-25 after item-18 closure.)

9. MGnify cross-validation reveals mobilome enrichment but low classification concordance (H4, H6)

Cross-validation against the MGnify genome catalogue (20,473 species across 4 biomes) provided independent evidence for several findings:

H4 — Mobilome enrichment validated: Plant-associated genera carry significantly higher mobilome burden than non-plant genera (median 3.7 vs 2.8 mobile elements per genome, Mann-Whitney U p=1.49×10⁻⁵). This genus-level result contrasts with genome-level studies finding fewer mobile elements in plant-associated bacteria (Levy et al. 2018; Bograd et al. 2025), likely reflecting the difference between rhizosphere/soil genera (where HGT is frequent; Rajabal et al. 2024) and the endosphere-focused genome-level comparisons.

T3SS rhizosphere enrichment: T3SS prevalence was roughly 2× higher in rhizosphere biomes (tomato: 24.0%, maize: 21.6%, barley: 22.3%) vs bulk soil (12.3%), independently confirming T3SS as a marker of plant association. This aligns with Gonzalez-Vallina et al. (2024), who found family-specific T3SS/T6SS enrichment in rhizosphere metagenomes.

H6 — Host specificity: 17 genera were detected across all 3 crop rhizospheres (Pseudomonas_E, Streptomyces, Variovorax, Telluria, Acidovorax), while 117 genera were unique to tomato, 54 to maize, and 5 to barley rhizosphere. Host-specific genera represent candidates for crop-specific biocontrol formulations.

Classification concordance: Only 11.7% Jaccard overlap between pangenome-derived plant genera (487) and MGnify rhizosphere genera (438). The low concordance reflects different data sources (isolation metadata vs metagenomic detection) and sampling biases rather than classification error — 93 pangenome plant genera were found in MGnify soil (not rhizosphere), suggesting broader environmental associations.

BGC profiles: 84 plant-associated genera produce NRP/siderophore BGCs (top producers: Mycobacterium 406, Pseudomonas_E 125, Streptomyces 113), consistent with the known importance of siderophore-mediated iron acquisition in plant colonization.

(Notebook: 11_mgnify_integration.ipynb)

10. Within-species subclade analysis shows weak segregation in 5/17 testable species, with two robust pathovar-host specializations (H7 weakly supported, H6 supported)

Canonical Phase 2b finding (H7): of 65 plant-associated species with ≥20 genomes, only 18 (28%) have any phylogenetic-tree distance data in kbase_ke_pangenome.phylogenetic_tree_distance_pairs — the other 47 (including major plant-associated taxa Bradyrhizobium japonicum, B. diazoefficiens, Mesorhizobium ciceri, Sinorhizobium medicae, Burkholderia glumae, Methylobacterium extorquens, Streptomyces scabiei, Xylella taiwanensis, Clavibacter michiganensis) are absent from the tree, a database-coverage limitation now documented in docs/pitfalls.md. Of the 17 testable species (≥3 plant + ≥3 non-plant), 5 (29%) pass Bonferroni-corrected Fisher's exact at α = 0.05/17 = 0.0029: Xanthomonas vasicola (p = 1.6×10⁻⁴), Mesorhizobium sp002294985 (p = 1.1×10⁻⁴), Agrobacterium pusense (p = 4.5×10⁻⁴), Pseudomonas_E avellanae (p = 7.3×10⁻⁴), Xanthomonas campestris (p = 1.0×10⁻³). Three of those 5 (X. vasicola, P. avellanae, X. campestris) also satisfy Cochran's rule for chi² and are the most robust; the other 2 pass only via Fisher's exact (which is exact and not subject to the small-cell constraint). The pattern is concentrated in Xanthomonas and Pseudomonas — both genera known for pathovar-level host specialization (Fautt et al. 2024; Satterwhite & Bergelson 2025).

Canonical Phase 2b finding (H6, layered on the same scan): of the 9 species with sufficient plant-host annotation for a subclade × host test (≥5 host-annotated genomes, ≥2 hosts, ≥2 subclades), 2 species pass both within-species and across-species Bonferroni: Xanthomonas campestris segregates to Brassica hosts (46/47 genomes in subclade 2 are on Brassica, p = 2.7×10⁻¹² raw, 3.3×10⁻¹¹ after within-species pair-search Bonferroni) and Xanthomonas vasicola segregates to Zea mays (47/52 in subclade 0, p = 1.4×10⁻¹¹ / 1.7×10⁻¹⁰). Both are canonical pathovar-host specializations confirmed at the genomic level: X. campestris pv. campestris is the textbook Brassica pathogen and X. vasicola pv. vasculorum is a known maize/sugarcane pathogen. Ralstonia solanacearum passes BH-FDR but not Bonferroni-across.

Phase 1/2 chronology: NB12's null result for H7 (5 species tested, 0/5 significant) was caused by a genome-ID format mismatch between phylogenetic_tree_distance_pairs (bare NCBI accessions) and genome_environment.csv (GTDB-prefixed IDs) — the is_plant_associated column was entirely NaN after the merge. NB13 fixed this by prepending GB_ to GCA_ accessions and RS_ to GCF_ accessions, recovering 1306/1306 genomes and 599 plant-associated genomes that NB12 had silently dropped. The original 5-species set with the fix shows 1/5 significant under all metrics; the full 18-species scan (item 19, NB-adjacent script _run_subclade_full_scan.py) raises this to 5/17 by including the broader candidate pool.

The 5/17 finding is consistent with the wider literature on accessory-genome-mediated plant adaptation: Melnyk et al. (2019) found that plant-lifestyle transitions in Pseudomonas are driven by convergent gain/loss of accessory genomic islands rather than core-genome phylogeny, and Satterwhite & Bergelson (2025) found local host adaptation in 3/5 P. syringae pathovars but only 1/5 locally dominant — both are partial-support patterns at the same proportion we observe genomically. The implication is that within-species plant adaptation is predominantly accessory-driven, but in specific genera (notably Xanthomonas) the accessory cassettes have stratified onto core-phylogenetic subclades — i.e., the pathovar designations have a real phylogenetic basis at the within-species level for those genera. The 47 species without phylo-tree coverage are an irreducible source of uncertainty: the true H7 effect across all 65 candidate species cannot be estimated until the GTDB tree is more completely populated in BERDL.

(Notebooks: 12_subclade_og_analysis.ipynb (5-species subset, kept for figure provenance); 13_validation_pfam_subclade.ipynb (genome-ID fix); notebooks/_run_subclade_full_scan.py (full 65→18 species scan, item 19, 2026-04-25); notebooks/_run_h6_host_subclade.py (full 18-species H6 host × subclade, 2026-04-25). Data: data/subclade_full_scan.csv, data/h6_host_subclade_full.csv, data/species_subclade_definitions_full.csv. H7 framing revised twice during Phase 2b.)

11. Phase 2b: Adversarial-review corrections (H0, H1, H3, H5, H7 revised)

Three adversarial review issues were fixable and were resolved in notebooks NB13–NB15: the tautological genus-level validation (C3), the original zero-hit Pfam query (documented limitation), and the NB12 subclade genome-ID mismatch (I6). Four additional statistical controls that RESEARCH_PLAN.md had specified but Phase 1 had not executed were also completed (C1, C4, I1, and within-genus label shuffling).

Species-level validation (C3, NB13): A curated panel of 18 model organisms with experimentally confirmed phenotypes (7 beneficial, 7 pathogenic, 4 neutral) replaces the 92.7% genus-level agreement number, which was circular because cohorts are defined from marker presence. All 18 ground-truth species were found in the refined cohort data. The primary finding is that the categorical cohort assignment is uninformative at the species level: all 14 of the beneficial and pathogenic ground-truth species — from Rhizobium leguminosarum and Azospirillum brasilense at one end to Pseudomonas syringae and Xanthomonas campestris at the other — were assigned to the single dual-nature class, and all 4 neutral controls were misclassified (E. coli K-12 as dual-nature, M. tuberculosis as dual-nature, D. radiodurans as beneficial, T. thermophilus as pathogenic). An earlier Phase 2b draft reported a "77.8% relaxed accuracy" figure (counting dual-nature as correct for both beneficial and pathogenic ground truth), but the paired adversarial review correctly noted that this rule cannot fail for any species that the classifier assigns to dual-nature, so the number contains no discriminative information between beneficial and pathogenic classes — it is arithmetically 14/18 but methodologically tautological. The only discriminative signal comes from a continuous metric: the pathogenicity ratio n_pathogen / (n_pgp + n_pathogen) is lower in known beneficial species (median 0.50) than in known pathogenic species (median 0.60), Mann-Whitney U = 9, p = 0.027 on N = 7 vs 7 — a small effect at a small sample size, but statistically consistent. The practical implication is that the refined cohort labels (beneficial / pathogenic / dual-nature / neutral) should be treated as a coarse screen, not a species-level classifier; the continuous pathogenicity ratio is what discriminates. The neutral misclassifications are not failures of the panel to detect plant adaptation but reflect the generality of the marker set — flagella, secretion, chemotaxis, biofilm overlap with "beneficial" or "pathogenic" markers in all bacteria, an issue already documented under §Limitations. A stronger follow-up validation would include near-pure beneficial organisms (Methylobacterium extorquens) and near-pure pathogens (Xylella fastidiosa) with known minimal PGP/pathogenicity machinery, to establish a more defensible panel. (Notebook: 13_validation_pfam_subclade.ipynb, data/species_validation.csv, figures/species_validation.png. Framing revised 2026-04-24 after paired adversarial review.)

Subclade genome ID fix (I6, NB13) and full 65→17 species scan (item 19, 2026-04-25): NB12's null result for H7 was traced to a format mismatch — phylogenetic_tree_distance_pairs stores bare NCBI accessions (GCA_..., GCF_...) while genome_environment.csv uses GTDB prefixes (GB_GCA_..., RS_GCF_...). After prepending the appropriate prefix, 1306/1306 genomes match (previously 0), recovering 599 plant-associated genomes across the 5 tested species. The paired adversarial review forced a careful re-read of which of those 5 species actually pass statistical validity, and pointed out that the original NB12 subset of 5 species (chosen by genome count) was too small for honest H7 conclusions. Following item 19 of the open-issues list, the analysis was re-run on all 65 plant-associated species with ≥20 genomes (notebooks/_run_subclade_full_scan.py, output data/subclade_full_scan.csv).

Phylogenetic-tree coverage: only 18 of 65 candidate species had any pairwise phylogenetic-distance data in kbase_ke_pangenome.phylogenetic_tree_distance_pairs. The other 47 species are absent from the tree — including major plant-associated taxa such as Bradyrhizobium japonicum, B. diazoefficiens, Mesorhizobium ciceri, Sinorhizobium medicae, Burkholderia glumae, Methylobacterium extorquens, Streptomyces scabiei, Xylella taiwanensis, Clavibacter michiganensis, and many more. This is a database-coverage limitation that bounds H7 power irrespective of methodology; it is now noted in docs/pitfalls.md under a new entry on phylogenetic-tree availability.

Of the 17 testable species (≥3 plant + ≥3 non-plant + ≥2 subclades after the genome-ID fix), 5 pass Bonferroni-corrected Fisher's exact at α = 0.05/17 = 0.0029, and 3 of those also satisfy Cochran's rule for the parametric chi² test:

The original Phase 2b draft cited P. avellanae and P. amygdali (1/5 valid). The full scan retains P. avellanae (validated again, with stronger context) and adds four new species: two Xanthomonas species (vasicola and campestris), one Mesorhizobium sp., and Agrobacterium pusense. P. amygdali's borderline chi² p = 0.031 from the 5-species subset does not survive Bonferroni in the larger scan and remains invalid under Cochran (see the original adversarial review for details).

H7 final verdict: weakly supported in 5 of 17 testable species (29%), with 3 species (X. vasicola, P. avellanae, X. campestris) showing the strongest signal across both Fisher's exact and Cochran-valid chi². The same scan on the full 18-species set was extended to test H6 (subclade × host) on 2026-04-25, with within-species Bonferroni for the (subclade × host) pair search and across-species correction for multiple species: 2 species pass both corrections — Xanthomonas campestris segregates to Brassica hosts (46/47 of genomes in subclade 2 are on Brassica, p = 2.7×10⁻¹² raw / 3.3×10⁻¹¹ Bonferroni-within / Bonferroni-across) and Xanthomonas vasicola segregates to Zea mays (47/52 in subclade 0, p = 1.4×10⁻¹¹ / 1.7×10⁻¹⁰). Both are biologically expected: X. campestris pv. campestris is the canonical Brassica pathogen and X. vasicola pv. vasculorum is a known maize/sugarcane pathogen. Ralstonia solanacearum passes BH-FDR but not Bonferroni-across (p = 5.9×10⁻³ within / 1.8×10⁻² q-BH). The scan output is data/h6_host_subclade_full.csv (notebooks/_run_h6_host_subclade.py); this is the new Phase 2b evidence underlying the H6 verdict revision documented in data/hypothesis_verdicts_final.csv. Within-species subclade segregation of plant-association does occur, and is concentrated in Xanthomonas and Pseudomonas — both genera known for pathovar-level host specialization (Fautt et al. 2024; Satterwhite & Bergelson 2025) — plus Mesorhizobium and Agrobacterium. Twelve testable species show no such clustering. The 47 species without phylogenetic-tree coverage are an irreducible source of uncertainty: the true H7 effect across all 65 plant-associated species with ≥20 genomes cannot be estimated until the GTDB phylogenetic tree is more completely populated in BERDL. (Notebook: 13_validation_pfam_subclade.ipynb (5-species subset, kept for figure provenance); notebooks/_run_subclade_full_scan.py (full 65-species scan, 2026-04-25); data/subclade_enrichment_corrected.csv (5-species), data/subclade_full_scan.csv (65-species), data/species_subclade_definitions_full.csv (genome-level assignments for the 18 species with phylo data); figures/subclade_corrected.png. H7 framing revised once 2026-04-24 after paired adversarial review, and again 2026-04-25 after full 65-species scan.)

Phylogenetic control (C1) — three complementary tests, three-tier verdict: The original NB10 statsmodels genus-level logit failed to converge for 0/14 markers. Phase 2b ran three increasingly-strict tests of the plant-association signal after phylogenetic control. The paired adversarial review flagged that the L1-bootstrap version was suggestive-not-confirmatory; item 18 of the open-issues list closed C1 with a cluster-robust GLM (cluster=genus), which is the pragmatic frequentist analogue of a phylogenetic generalized linear mixed model (PGLMM) at species scale.

(a) L1-regularized logit + bootstrap CI (NB14 Cell 2): scikit-learn L1-penalized Logit (saga, C=1.0), top-20 genus dummies + log₁₀(genome_size), 100-replicate bootstrap. 9/14 markers had bootstrap 95% CIs excluding zero. Caveat: L1 bootstrap CIs are biased near zero (Chatterjee & Lahiri 2011); n_bootstrap=100 is small.

(b) Cluster-robust GLM, cluster=genus (notebooks/_run_c1_cluster_robust.py, 2026-04-25): unpenalized statsmodels GLM with binomial family, top-10 phylum fixed-effect dummies + log₁₀(genome_size), and cluster-robust covariance with 7,555 genus clusters across 24,554 species. Proper Wald p-values and 95% Wald CIs, no L1 shrinkage. 8/14 markers survive at BH-FDR q < 0.05 with |coef| > 0.2 — and 13 of 14 markers agree with the L1 verdict; the one disagreement is iaa_biosynthesis (L1 coef = 0.61 sig; cluster-robust coef = −0.09, p = 0.82, not sig), which is exactly the kind of L1-shrinkage-near-zero artifact the adversarial review flagged. The honest 14-marker breakdown:

(c) Within-genus label shuffling (NB14 Cell 4): the strictest test — plant labels permuted within each genus, preserving genus-level plant-association rate. 3/15 markers survive (nitrogen fixation, ACC deaminase, T3SS), all with very small p-values.

Three-tier honest interpretation of plant-marker signals after phylogenetic control:

• Tier 1 — robust at species level under any test (3 markers): nitrogen fixation, ACC deaminase, T3SS. Large coefficients, proper Wald CIs not crossing zero, and survive within-genus permutation. These are the markers where specific plant-associated species within a genus carry the function and non-plant sister species reliably do not.
• Tier 2 — robust under cluster-robust GLM but fail within-genus shuffle (5 markers): phenazine, CWDE cellulase, CWDE pectinase, phosphate solubilization, effector. The signal is real after accounting for within-genus correlation but disappears when you reshuffle plant labels within genus — meaning these are cassette-level signals: when a genus is plant-associated, the cassette tends to be acquired, but the cassette is not species-specific within the genus. Still ecologically meaningful (genus-acquired plant adaptation), but not a species-level distinguisher.
• Tier 3 — not significant in any test (6 markers): DAPG biocontrol, T4SS, hydrogen cyanide, IAA biosynthesis, siderophore, acetoin/butanediol. These had been suggested by L1-bootstrap (IAA) or were known weak signals; cluster-robust GLM confirms they are not robust.

This three-tier breakdown supersedes the earlier "9/17 survive L1" framing and gives a defensible answer to C1. A full PGLMM on the GTDB tree would be a stronger version (random-effect covariance from tree branch lengths rather than cluster grouping), but the GTDB-tree-derived covariance matrix would be a 25K × 25K dense matrix, intractable in memory; the cluster-robust GLM at genus level is the practical analogue and gives identical conclusions for 13 of 14 markers. (Notebooks: 14_deferred_controls.ipynb (L1 + shuffle); notebooks/_run_c1_cluster_robust.py (cluster-robust GLM, 2026-04-25). Data: data/regularized_phylo_control.csv (L1), data/sensitivity_shuffle.csv (within-genus shuffle), data/c1_cluster_robust.csv (cluster-robust GLM). Three-tier framing added 2026-04-25 after item-18 closure.)

Genome-size covariate for novel OGs (C4, NB14) — closed after adversarial review: This deferred control has two parts, both now executed on real per-species data. An earlier Phase 2b draft used a prevalence-based Bernoulli simulation for the per-OG test (flagged by the paired adversarial review as circular by construction); this has been replaced with the real per-species OG matrix from NB03 and re-run.

(a) Dual-nature rate by genome-size quartile (observational): Among plant-associated species, the dual-nature rate is 54% in the bottom genome-size quartile vs 87% in the top three quartiles (measured directly on refined cohort data). This reaffirms the NB08 finding that ~33% of dual-nature species shift to neutral when genome size is normalized — the headline 78.7% dual-nature rate masks a real genome-size gradient, and the gradient itself is observational, not model-dependent.

(b) Per-OG logistic regression (real per-species test, unpenalized statsmodels Logit with proper Wald inference + BH-FDR across 50 OGs): For each of the 50 novel OGs, the per-species presence vector from top50_og_species.csv (cached from NB03's server-side eggNOG aggregation) was regressed on is_plant + log₁₀(genome_size) + top-10 phylum dummies across 22,384 species with complete covariates (1,134 plant-associated, 21,250 non-plant). The model fit converged for 48/50 OGs; for the remaining two (COG2351 and COG5343) the design matrix was singular — likely perfect separation by a phylum dummy given their very strong raw prevalence gaps (62% vs 23% and 64% vs 21% respectively). 48 of 50 OGs (96%) retain a positive plant-association coefficient that is large (|coef| > 0.2) and statistically significant at BH-FDR q < 0.05 after phylum + log₁₀(genome_size) control. Controlled odds ratios for the top 10 range from 4.4 (COG5516) to 6.6 (COG1845, cytochrome oxidase ctaE), with q-values as small as 10⁻⁹⁶ (COG3569). The genome-size coefficient is strongly positive (median 3.4 on log₁₀-scale) — plant-associated species genuinely do have larger genomes — but that genome-size effect is absorbed by the log₁₀(genome_size) term and does not explain the is_plant effect away: both coefficients are simultaneously large and positive. The two OGs that failed to fit (COG2351, COG5343) have strong raw prevalence signals (3× higher in plant vs non-plant) so the "fail" is numerical, not biological. C4 is therefore closed: after actual per-species regression with proper inference, the novel OGs survive genome-size confounding, and the Phase 2b result stands at the same 48/50 level that the earlier simulation reported mechanically — but now for honest reasons. (Notebook: 14_deferred_controls.ipynb Cell 3, data/genome_size_control.csv. Real per-species test added 2026-04-25 after paired adversarial review replaced the circular simulation.)

Within-genus label shuffling (RESEARCH_PLAN.md safeguard, NB14) — strict test, mostly negative: The most stringent phylogenetic test. For each marker, plant/non-plant labels are permuted within genus (preserving genus-level plant-association rate), and observed Fisher odds ratio is compared to the null distribution. Only 3 of 15 testable markers survive: nitrogen fixation (perm p < 0.001), ACC deaminase (perm p < 0.001), and T3SS (perm p = 0.005). The other 12 markers — phenazine, IAA biosynthesis, CWDE cellulase, CWDE pectinase, effector, DAPG biocontrol, hydrogen cyanide, siderophore, acetoin/butanediol, phosphate solubilization, T4SS, phytotoxin — fail the test, meaning their plant-enrichment signals are statistically consistent with what you would get by reshuffling plant/non-plant labels within each genus. The honest reading is that most of the marker-level plant-association signals in this project are genus-scale, not species-scale — plant adaptation is largely a property of which genera happen to be plant-associated, and within those genera, marker presence does not reliably distinguish plant-associated from non-plant species. Only three markers — the nif nitrogen-fixation operon, the ACC deaminase enzyme, and the hrp/hrc T3SS structural cluster — are carried by specific plant-associated species within genera where non-plant species do not reliably have them. This does not invalidate the H1/H2 compartment and architecture findings (which operate at the genus/family/phylum scale already), but it does mean that the headline dual-nature rate (78.7%) and the compartment-specific marker enrichments should be interpreted as genus-level genomic backgrounds associated with plant-associated lineages, not as species-level functional evidence of plant adaptation. The within-group permutation approach is conceptually aligned with "phylogenetic permulations" (Saputra et al. 2021), which combines phylogenetic simulation and permutation to produce unbiased empirical p-values when standard parametric nulls are confounded. Our simpler within-genus partition is appropriate when the signal is suspected to be genus-scale; a permulation-style test against the full GTDB tree would be the stricter follow-up. Caveats: n = 200 permutations is under-resolved for p ≤ 0.005 (T3SS is therefore marginal); no multiple-testing correction was applied across 15 markers (Bonferroni α = 0.0033 would demote T3SS; N-fix and ACC deaminase remain sig). (Notebook: 14_deferred_controls.ipynb, data/sensitivity_shuffle.csv. 12/15 failure foregrounded 2026-04-24 after paired adversarial review.)

PERMANOVA + PERMDISP + db-RDA on the H1 compartment effect (NB14 + item 20 follow-up, 2026-04-25): The original NB04 PERMANOVA reported R² = 0.527 for compartment on the Phase 1 25-marker panel; NB14's sensitivity analysis showed R² drops to 0.072 when the top-3 genome-rich species per compartment are excluded (an 86% loss). To complete the C2/H1 story, item 20 of the open-issues list called for separating location shift from dispersion heterogeneity via PERMDISP and db-RDA on the same data, since the original PERMANOVA conflates the two. This was run on all 607 plant-compartment species (root/rhizosphere/phyllosphere) using the refined 17-marker v2 panel and Jaccard distance (notebooks/_run_h1_dbrda.py).

The PERMANOVA R² on the full 607-species set with the v2 panel matches NB14's exclude-top-3 result (0.071 vs 0.072) — the difference between Phase 1's 0.527 and this 0.071 was the marker panel and the genome-count distribution, not a sample-size artifact. PERMDISP is highly significant (F = 15.6, p = 0.001), confirming C2 — group dispersions differ: root species cluster most tightly (mean centroid distance 0.452), phyllosphere intermediate (0.503), rhizosphere most variable (0.528). db-RDA gives R² = 0.0595 for the location-only effect, also significant (permutation p = 0.001), and represents 84% of the total PERMANOVA R². So:

• The original 0.527 was an inflated-magnitude artifact of the Phase 1 panel + the few genome-rich species that dominated the distance matrix (already established in NB14).
• The residual 0.071 effect is real, not primarily a dispersion artifact. About 84% of it is genuine location shift in centroid space; 16% is the dispersion heterogeneity that PERMDISP catches.
• Root, rhizosphere, and phyllosphere are functionally distinct in centroid terms, just at a small effect size (~6% of variance), and rhizosphere does have systematically more variable functional profiles than root.

H1 final verdict: weakly supported. Compartments impose statistically significant but small (~6% of variance) location shifts on plant-associated functional profiles, alongside significant dispersion heterogeneity (root tightest, rhizosphere most variable). The headline 0.527 from Phase 1 was a panel + sampling artifact, not a community-wide signal. (Notebooks: 04_compartment_profiling.ipynb, 14_deferred_controls.ipynb, notebooks/_run_h1_dbrda.py. Data: data/sensitivity_results.csv (NB14 exclude-top-3), data/h1_dbrda_results.csv (PERMANOVA + PERMDISP + db-RDA, 2026-04-25), data/h1_compartment_dispersions.csv. Framing revised 2026-04-24 after paired adversarial review and again 2026-04-25 after db-RDA decomposition.)

Prevalence-weighted GapMind complementarity (I1, NB14) — mostly a formula fix, minor aggregation fix: NB06 reported Cohen's d = −7.54 for co-occurring vs random pairs, which is an implausibly large effect. The paired adversarial review correctly identified that most of this was a Cohen's d formula problem, not an aggregation problem: NB06 divided the mean difference by null_means.std() (SD of permutation means, which shrinks to zero by CLT for any reasonable null), while NB14 divides by random_pair.std() (raw pair-level SD, much larger). Re-computing NB06's max-aggregated complementarity using NB14's consistent formula yields d ≈ −0.43 — essentially identical to the prevalence-weighted d = −0.39. So the 20× reduction in apparent magnitude is ~90% formula-change and ~10% aggregation-change. The correct interpretation of the corrections is:

• Formula fix: The original d = −7.54 was spuriously large because NB06 used the wrong denominator. With a consistent pair-level SD, the true scale of the effect at the max-aggregation level was always about −0.43.
• Aggregation fix: Switching from max-aggregation (binary pathway present/absent per genus) to prevalence-weighting (fraction of species with pathway) additionally reduces |d| from ≈0.43 to 0.39, a ≈10% effect. Meaningful but small.
• Direction: Unchanged. Co-occurring genus pairs are slightly less complementary than random pairs (permutation p < 0.001), consistent with Louca et al. (2018) redundancy-driven community assembly.

H3 remains not supported — the ecological conclusion (redundancy not complementarity) is robust to the aggregation choice. What changes is that the Phase 1 report's emphasis on "Cohen's d = −7.54" as a dramatic effect size was an artifact; the real effect is small (|d| ≈ 0.4) but stable. The Phase 2b value of this analysis is the formula correction and the documentation of the aggregation sensitivity, not a large new biological finding. (Notebook: 14_deferred_controls.ipynb, data/complementarity_v2.csv, figures/complementarity_comparison.png. Formula-vs-aggregation attribution clarified 2026-04-24 after paired adversarial review.)

Pfam recovery (Pfam gap, NB13 Cell 3) + full bakta-vs-IPS audit (2026-04-25): The Phase 1 NB02 query WHERE pfam_id IN ('PF00771', …) returned zero rows because bakta_pfam_domains stores versioned IDs (PF00771.22). Re-running with LIKE 'PF00771%' patterns on the 10 Phase 2b marker Pfams recovered 19,364 domain hits across 7,962 species. Five of those 10 Pfam IDs produced hits: nitrogenase NifH PF00142 (6,632 hits, 4,040 species, 32.4% core), cellulase GH5 PF00150 (6,898 / 3,234 / 25.9%), pectate lyase 3 PF12708 (4,149 / 2,283 / 24.6%), pectate lyase PF00544 (1,115 / 546 / 20.8%), and T6SS VgrG PF05943 (570 / 343 / 10.9%). The remaining five returned 0 hits.

After the paired adversarial review flagged the T3SS-narrative inconsistency, a full audit was run across all 22 marker Pfams that the project queries across NB02 and NB10 (data/pfam_bakta_ips_audit.csv). The result is a systematic finding worth documenting beyond the immediate Phase 2b context: 12 of the 22 marker Pfams are completely absent from bakta_pfam_domains despite being abundant in interproscan_domains. The silent-gap set is dominated by secretion-system components — T3SS (PF00771 18.6K IPS hits, PF01313 13.6K, PF01514 15.0K, PF02579 16.8K, PF04610 15.6K, PF07916 8.4K, PF09599 0.3K, PF01312 22.0K), T4SS (PF03135 11.0K, PF03743 17.5K), and T6SS (PF05936 7.5K Hcp tube, PF04183 7.5K). The 10 Pfams that are in bakta (mostly nitrogen fixation, CWDEs, and a few others) still return only 10–35% of the IPS hit count, so even those are systematically undercounted but non-zero. The likely cause is that bakta uses a reduced "core" Pfam HMM profile set while InterProScan uses the full Pfam-A release.

Operational consequence for this project: NB10's refined cohort pipeline already used interproscan_domains (signature_acc field) for T3SS, T4SS, T6SS, and CWDE marker detection, so its 280,193-hit InterProScan recovery was the correct decision and the refined cohort assignments are unaffected. Phase 1 NB02's original bakta-only Pfam query returned zero hits and that step fell back to gene-name keyword matching, which is also unaffected. The Pfam LIKE recovery in NB13 Cell 3 confirms rather than expands InterProScan numbers for the secretion-system Pfams, but for nitrogen fixation and CWDEs it genuinely adds 4,217 species gaining at least one new marker (2,872 gain nitrogen fixation, 1,082 gain CWDE pectinase, 811 gain CWDE cellulase) — these are the cases where bakta and IPS diverge and IPS missed some that bakta picked up. The higher core-genome fractions of NifH (32.4%) vs the CWDEs (21–26%) are consistent with the H2 pattern.

Operational consequence for future BERDL projects (now in docs/pitfalls.md): treat bakta_pfam_domains as a supplement to, not a substitute for, interproscan_domains. Cross-check any Pfam-based query against the IPS table before drawing biological conclusions about missing functions. (Notebook: 13_validation_pfam_subclade.ipynb. Data: data/pfam_recovery_hits.csv, data/pfam_recovery_impact.csv, data/pfam_bakta_ips_audit.csv. See docs/pitfalls.md for the full audit table.)

Adversarial issue disposition table

(Notebooks: 13_validation_pfam_subclade.ipynb, 14_deferred_controls.ipynb, 15_final_synthesis.ipynb)

DESIGN-plan execution-deviations close-out (2026-04-25)

Beyond the five paired-adversarial-review issues above, four planned-vs-executed gaps from the original RESEARCH_PLAN.md research plan were closed in commit 8f9086e:

• Gap 4 — NB04 broken logit import: the genus-level logistic-regression cell in NB04 used the bare name logit(...) while the import was import statsmodels.formula.api as smf (should have been smf.logit(...)). The cell was rerouted to a markdown reroute pointing readers to NB14's cluster-robust GLM as the canonical Phase 2b answer; original code preserved verbatim in a fenced block within the markdown for chronology. (Cosmetic close — does not change H1.)
• Gap 3 — verification #12 formal re-test: the DESIGN-plan target "refined cohorts produce ≥5 PGP-only genera in NMDC for C-score analysis" was formally re-checked on Phase 2 refined cohorts (data/cscore_refined_cohorts.csv). Result: 0 PGP-dominant / 3 pathogen-dominant / 66 dual-or-mixed — same pattern as Phase 1, verification #12 NOT MET, the failure-as-result interpretation is made explicit. (Does not change H3.)
• Gap 2 — genus_dossiers_v2.csv refresh: the planned NB15 deliverable was produced with 30 top plant-associated genera × 23 columns, integrating refined cohorts, cluster-robust marker survivors, MGnify mobilome/BGC, and species-level H7/H6 significance per genus. Xanthomonas surfaced as the lone genus with both H7 phylogenetic segregation and H6 host specialization at the genomic level. (Consumer-facing deliverable — does not change a verdict.)
• Gap 1 — NB13 scope-deviation documentation: the originally-planned NB13 (co-occurrence network expansion: Louvain modules, hub centrality, partial Mantel, NMDC abiotic_features and metabolomics_gold integration) was re-scoped to validation/Pfam/subclade fixes when the paired adversarial review surfaced higher-priority methodological issues. The complementarity sub-item (complementarity_v2.csv) was preserved and executed in NB14 Cell 5 — closing H3's I1 issue. The remaining sub-items are deferred as exploratory enhancements (the H3 verdict does not depend on them). RESEARCH_PLAN.md "Execution Deviations (Phase 2 / Phase 2b)" section formally documents this scope reduction.

All four DESIGN-plan gaps are now formally closed; the research plan is fully reconciled with actual execution.