Metal Resistance Ecology — Phylogenetic Conservation vs. Environmental Selection

Raw report: microbeatlas_metal_ecology Feeds into: adversarial-research-quality-assurance, callability-limited-comparative-inference, chromosomal-and-integrative-gene-transfer, classifier-database-compatibility-in-taxonomic-quantification, composite-resistance-score-limitations, confirmatory-exploratory-ecological-association-discordance, cultivation-collection-bias-in-ecological-genomics, ecotype-environment-gene-content, environment-embedding-geography, environmental-resistome, functional-marker-validation, gene-function-acquisition-depth, genomic-under-representation, metal-cross-resistance, ontology-and-category-schema-sensitivity, pangenome-integration, pangenome-openness-determinants, phenotype-database-coverage-bias, phylogenetic-confounding-of-pangenome-associations, sampling-depth-and-downsampling-effects, spatial-sampling-effort-confounding, study-batch-confounding-of-environmental-associations, subsurface-hydrogeological-zonation, taxonomic-nomenclature-reconciliation, taxonomic-resolution-dependent-functional-inference, within-species-conservation-between-species-functional-divergence

Overview

This study links genus-level metal-resistance type diversity, inferred from AMRFinderPlus pangenome annotations across 6,789 GTDB species, with global ecological niche breadth derived from a 464,000-sample 16S amplicon atlas (MicrobeAtlas). The metal-resistance repertoire considered spans mercury, arsenic, copper, zinc, cadmium, chromium, and nickel. The report describes the joined resources as a pangenome database of 6,789 sequenced species across the GTDB bacterial tree and a global OTU (operational taxonomic unit) atlas of 98,919 OTUs × 463,972 samples, a scope it says site-specific studies cannot reach. Across 1,264 bacterial genera with ≥ 3 OTUs, niche breadth was phylogenetically conserved, while phylogenetic generalized least squares (PGLS), which models trait associations while accounting for shared evolutionary history, tested whether metal-resistance traits predicted niche breadth beyond phylogeny. The analysis used 606 bacterial genera with AMR data and found that metal type diversity, rather than total AMR gene burden or core AMR fraction, was associated with broader inferred ecological ranges. ₁

The primary evidence is correlational: the cross-sectional comparison cannot distinguish whether metal tolerance enables habitat expansion, broad ecological range increases opportunities for metal-resistance gene acquisition, or a shared factor such as genome size, metabolic versatility, or biofilm capacity promotes both traits. The report therefore concludes that genera with diverse metal-resistance repertoires are, on average, broader ecological generalists in the global microbiome, not that metal-resistance diversity causes ecological generalism. It describes the PGLS coefficient as a partial correlation after phylogenetic correction, not a causal estimate. Distinguishing the scenarios would require time-series genomic data tracking metal-resistance evolution alongside range shifts, or natural-experiment designs such as metal enrichment of defined microbial communities with longitudinal tracking. The report argues that the pattern of results slightly favours habitat expansion or a shared driver over a simple model in which genome size drives everything, because total AMR cluster count is non-significant while metal type diversity is significant, and cluster count scales more directly with total genome size. It also states that metal type diversity remaining independently significant after formal genome-size control (β = +0.022, p = 3.6×10⁻⁴) rules out genome size as the primary confounder. The report nonetheless treats the association as correlational, with its direction undetermined. It cites Qi et al. (2022) as independent evidence that multiple heavy metal contamination favours microbial generalists as network keystones in soil, and reads that soil result as consistent with the direction of its own association. ₁

The report, which calls this the first analysis of its kind, links genus-level metal type diversity, inferred from AMRFinderPlus pangenome annotations across 6,789 GTDB species, to global niche breadth from the 464,000-sample MicrobeAtlas 16S atlas. It describes a positive association that holds beyond what phylogenetic relatedness predicts and calls it robust to genome size, species richness, prevalence filtering, OTU abundance filtering, and all 13 environment categories. Two sensitivity results reported below qualify that broad robustness claim. The strict 5% within-environment prevalence analysis was not significant (p = 0.092), and none of the 7 leave-one-metal-out exclusions stayed significant. ₁

Figure 1 (fig1_lambda_heatmap.png) is a Pagel’s λ heatmap by trait and domain, covering 6 traits × Bacteria/Archaea on a YlGnBu colour scale. The report’s other figures are a PGLS forest plot of all simple models plus the multi-predictor model (Figure 2, fig2_pgls_forest.png, with the 6 simple models as circles and the multi-predictor model as triangles), a primary 2-panel synthesis figure combining a Pagel’s λ heatmap (Panel A) with a PGLS forest plot (Panel B) (Figure 3, fig3_synthesis.png), and a scatter of metal type diversity against Levins’ B_std with the top-5 genera annotated (Figure 4, fig4_metal_types_scatter.png). A robustness summary (Figure 5) shows β across 5 analysis scenarios, the rarefied β distribution over 200 iterations, and an archaeal power curve (file fig5_robustness.png). Other figure files are niche_breadth_distribution.png, showing the distribution of Levins’ B_std across MicrobeAtlas OTUs; pgls_forest_plot.png, a preliminary NB05 forest plot superseded by fig2; and pgls_scatter.png, a partial regression scatter by phylum from NB05. fig_enigma_validation.png has two panels: a Spearman scatter of metal-type diversity against groundwater prevalence (767 genera, ρ=+0.112, p=0.0019) and a Mann–Whitney box plot comparing top-Q4 and bottom-Q1 metal-diversity genera (p=0.007). fig_enigma_trackB.png has three PRJNA1084851 panels: community-weighted mean (CWM) metal-type diversity by well (Kruskal–Wallis p=0.0001), CWM against time point after carbon amendment (ρ=+0.38, p<0.0001), and Sulfurimonas relative abundance declining post-amendment (ρ=−0.32, p=0.0002); these panel captions give rounded versions of the ρ = +0.383 and ρ = −0.323 reported below. An interactive standalone Plotly dashboard (1 MB) maps 3,942 sample grid cells at 0.5° resolution, built from 386,020 samples with GPS coordinates in arkinlab_microbeatlas.sample_metadata and coloured by mean metal gene-type diversity per environment. It also plots all 3,160 genera by niche breadth and metal diversity and shows a Pagel’s λ bar chart with significance colouring for all traits. ₁

Key Findings

Phylogenetic conservation and the primary PGLS association

Bacterial Levins’ B_std niche breadth showed strong phylogenetic signal, with Pagel’s λ = 0.787 and p = 7.9×10⁻¹⁰² across 1,264 bacterial genera with ≥ 3 OTUs in MicrobeAtlas; habitat range, measured as the number of detected environment categories, was more conserved with λ = 0.909 and p = 1.4×10⁻¹⁵⁷. In the primary PGLS analysis of 606 genera, metal type diversity was the only metal AMR predictor surviving Bonferroni correction across 6 simple models (β = +0.021, SE = 0.0056, p = 1.5×10⁻⁴; threshold p < 0.0083). ₁

The PGLS analytical subset comprised 606 bacterial genera with ≥ 3 OTUs. Each genus was present in the GTDB r214 bacterial genus tree and had metal AMR data from the kbase_ke_pangenome dataset, and predictors were z-scored. PGLS was estimated jointly with λ by maximum likelihood (ape::gls with nlme::corPagel). Data were sorted into the tree’s tip-label order before fitting. ₁

A unit test (scripts/test_pgls_ordering.py) checked this ordering step. It generated synthetic data with known phylogenetic signal (λ = 0.8, β = 0.05) and scrambled the row order. It then verified that sorted PGLS recovered the correct λ and β within tolerance, while unsorted PGLS failed. The report records the test as passing. The test checks the software step on synthetic data. It does not test the empirical result. ₁

The multi-predictor model retained a positive metal type effect (β = +0.023, SE = 0.0067, p = 5.5×10⁻⁴), whereas total AMR cluster count (β = −0.002, SE = 0.0065, p = 0.789) and core AMR fraction (β = +0.006, SE = 0.0054, p = 0.288) were not significant. Each simple PGLS model used n = 606 genera. The simple estimates were β = +0.021 for metal types (λ = 0.708, SE = 0.0056, p = 1.5×10⁻⁴, ΔAIC = −12.4; Bonferroni significant), β = +0.010 for AMR clusters (λ = 0.700, SE = 0.0056, p = 0.090, ΔAIC = −0.9; not Bonferroni significant), and β = +0.003 for core fraction (λ = 0.701, SE = 0.0054, p = 0.523, ΔAIC = +1.6). Total AMR cluster count and core AMR fraction therefore did not predict niche breadth in either the simple or the multi-predictor models. The report reads this as showing that the breadth of the metal-resistance repertoire, rather than its depth (many genes for few metals), is what is associated with ecological versatility. The report states that genera with many AMR genes for a single metal type are not ecological generalists, but genera spanning multiple metal resistance families are. It presents this diversity-over-depth result as paralleling findings in other stress response domains. ₁

The three-predictor Levins’ B_std model had λ = 0.706 and AIC = −743.67. With number of environments as the outcome, no simple model passed the Bonferroni threshold. Core fraction had λ = 0.839, β = +0.142, SE = 0.069, p = 0.038, ΔAIC = −2.3, which is a nominal association only. AMR clusters had λ = 0.847, β = −0.071, SE = 0.067, p = 0.290, ΔAIC = +0.9, and metal types had λ = 0.848, β = −0.007, SE = 0.071, p = 0.923, ΔAIC = +2.0. Metal type diversity therefore predicted Levins’ B_std but not the number of environments detected. ₁

The raw OLS correlation between mean metal type diversity and Levins’ B_std was modest, r ≈ 0.21, which reflects the baseline phylogenetic structure in both variables. The PGLS coefficient of +0.021 per SD increase in metal type diversity is the additional covariation that remains after removing the shared-ancestry component estimated by the λ-transformed variance–covariance matrix. The result indicates that breadth across multiple metal types, rather than depth represented by many genes for a smaller number of metals, distinguishes broad-niche genera. The discussion describes the association as β = +0.021 SD per SD, Bonferroni-corrected. It calls the result consistent with a metabolic versatility hypothesis, in which genera with diverse metal tolerance repertoires, for example spanning Hg detoxification, As efflux and Cu oxidation systems, have broader ecological ranges. Metabolic versatility remains a hypothesis, not a tested mechanism. The ‘SD per SD’ wording conflicts with the per-SD-of-metal-diversity description above and with the model specification, in which predictors but not the B_std response were z-scored, so the outcome units of β are unresolved. ₁

Metal AMR traits have intermediate phylogenetic signal

Among 606 bacterial genera in GTDB r214, total AMR cluster count had Pagel’s λ = 0.260 and p = 6.1×10⁻²⁸, core AMR fraction had λ = 0.441 and p = 1.8×10⁻⁸, and metal type diversity had λ = 0.335 and p = 1.1×10⁻²³. The ordering, core fraction λ > metal type diversity λ > cluster count λ, is consistent with, but does not establish, a model in which constitutive (core) metal resistance is mainly vertically inherited, while horizontal transfer and local metal exposure shape accessory gene accumulation and metal-type expansion. In summary form, bacterial niche breadth and habitat range are more phylogenetically conserved than metal AMR traits (λ_niche ≈ 0.79–0.91 vs λ_AMR ≈ 0.26–0.44). The report reads this as ecological range being deeply canalised by lineage identity while metal AMR is more labile. It attributes that lability to AMR genes likely residing on mobile genetic elements and being frequently exchanged among distantly related taxa. This mobile-element explanation is an interpretation, not a direct measurement of horizontal gene transfer in this analysis. In support, the report cites Hemme et al. (2016), who demonstrated recombinational hotspots for metal resistance gene acquisition in contaminated groundwater. It also cites Gillieatt & Coleman (2024), who reviewed the mechanistic basis for AMR gene transfer by mobile genetic elements under metal co-selection. ₁

Nitrification provided a phylogenetic positive control: its signal was λ = 0.939 in bacteria with p = 2.5×10⁻¹⁰² and λ = 1.000 in archaea with p = 2.5×10⁻⁵¹. This contrasts with the intermediate signal for metal AMR and supports the interpretation that deeply inherited metabolic traits and mobile resistance traits have different phylogenetic architectures. The discussion summarises nitrification as λ = 0.94–1.00 in both domains. It describes ammonia- and nitrite-oxidation genes as deep, vertically inherited innovations that define major lineages (Nitrospirota, Thaumarchaeota). The report says this internal comparison strengthens confidence in the analytical pipeline. ₁

In archaea, Levins’ B_std had detectable phylogenetic signal across 132 genera (λ = 0.197, p = 1.1×10⁻⁵), and the number of environments was strongly conserved (λ = 0.898, p = 4.6×10⁻¹⁴). The report notes that archaeal B_std λ is notably lower than in bacteria (0.787), consistent with the phylogenetically sparse and ecologically distinct archaeal genera sampled by short-read 16S amplicon surveys. Archaeal AMR cluster count showed approximately zero signal across 48 genera (λ ≈ 0, p = 1.0). Pagel’s λ was computed with phytools::phylosig (method = ‘lambda’, test = TRUE) against GTDB r214 genus-representative trees. ₁

Standalone Pagel’s λ for bacterial B_std was 0.787, whereas the PGLS-estimated λ was 0.708, a discrepancy of Δλ = 0.079. The report interprets the lower residual λ as expected because metal type diversity itself has phylogenetic signal (λ_metal_types = 0.335) and removes some phylogenetically structured variance; residual λ = 0.708 remained strongly non-zero. The report further reads the drop from 0.787 to 0.708 as a 10% reduction, meaning that metal type diversity accounts for approximately 10% of B_std’s phylogenetic covariance. On this reading, metal-resistance breadth is one of several phylogenetically structured predictors of niche breadth rather than the sole driver. The report gives residual λ = 0.708 with p ≪ 0.001 and finds no evidence of model misspecification from the difference. This variance-share figure is the report’s interpretation of the λ change, not a formal decomposition. ₁

Robustness and sensitivity analyses

The metal type association remained significant after adding pangenome coverage as a covariate, measured as the number of species with AMR data per genus (n_species): β = +0.0204, SE = 0.0057, p = 3.4×10⁻⁴ in the simple model and β = +0.0224, SE = 0.0067, p = 8.3×10⁻⁴ in the model containing all three AMR predictors. The coverage covariate had β = +0.0093, SE = 0.0051, p = 0.068 in the simple model and β = +0.0092, SE = 0.0051, p = 0.070 in the multi-predictor model. The report calls the covariate borderline (p ≈ 0.07), indicating a modest but non-dominant contribution of sampling depth. The limitations section gives the same R1 result rounded to β = +0.020, p = 3.4×10⁻⁴. It specifies that n_species_with_amr entered the PGLS itself (gls() + corPagel), not an OLS model. ₁

To remove genome-count bias, 200 PGLS rarefaction iterations each represented every genus by 1 randomly sampled species. Of the 606 PGLS genera, 383 had >1 species in the AMR database; the other 223 were already singletons, so rarefaction did not change them. All 200 / 200 iterations completed, and the median metal type coefficient was +0.0147, the IQR was +0.0126–+0.0166, the median p-value was 0.0054, 89.5% of iterations had p < 0.05, 57.5% were Bonferroni-significant at p < 0.0083, and the median λ was 0.704. The direction was consistent across all iterations, although the median coefficient was smaller than the full-data estimate of +0.0215. The report attributes this attenuation to reducing multi-species genera, which have higher inferred metal types, to one representative, and it concludes that the signal is not artefactual. The limitations section gives the R2 median coefficient rounded to +0.015, alongside the 89.5% of iterations with p < 0.05. ₁

The strict 5% within-environment prevalence analysis counted an OTU as present in an environment only if it appeared in ≥5% of samples in that category. This was meant to reduce detections driven by contamination, cross-sample bleed-through, or rare transient events. The filter reduced the number of OTUs from 98,919 to 10,433, the number of genera with B_std estimates from 3,160 to 1,120, and the PGLS sample size from 606 to 379. The coefficient remained positive but was not significant: β = +0.0166, SE = 0.0099, p = 0.092, λ = 0.526, compared with β = +0.0215, SE = 0.0056, p = 1.5×10⁻⁴ and λ = 0.708 in the original analysis; the Pearson correlation between original and strict B_std was 0.680. The report attributes the loss of significance primarily to the reduced sample size (n drops 37%: 606 → 379), not to a change in direction or magnitude, and describes the change in β from +0.0215 to +0.0166 as a 23% reduction. That attribution is interpretive, because the standard error also rose from 0.0056 to 0.0099. The report cautions that the 5% threshold is quite conservative: it excludes 89% of OTUs, many of which the report considers genuinely distributed across environments but uncommon within any one. It proposes a threshold of 1–2% as a reasonable middle ground for future work. ₁

For the genome-size analysis, genus-level mean genome size (bp) was retrieved from kbase_ke_pangenome.gtdb_metadata through the KBase Data Lakehouse REST API, using one representative species per genus. This covered 527/606 PGLS genera. Genome size was log-transformed and z-scored, then added as a PGLS covariate alongside metal type diversity. Because 79/606 genera (13%) lacked genome size data at query time, PGLS n dropped from 606 to 527; the report states that results are consistent in direction and significance with the full-data models. Controlling for genome size, the metal type coefficient was β = +0.0218, SE = 0.0061, p = 3.6×10⁻⁴, λ = 0.629, while log(genome size) had β = +0.0263, SE = 0.0069, p = 1.5×10⁻⁴. The baseline genome-size-only model had β = +0.029, p = 2.7×10⁻⁵, AIC = −630.3; the full model had AIC = −641.1, giving ΔAIC = 10.8 in favor of including metal type diversity. ₁

In the three-covariate PGLS with metal type diversity, log(n_species), and log(genome size), using 527 genera with all covariates, metal types had β = +0.0218, SE = 0.0061, t = 3.60, p = 3.5×10⁻⁴; log(n_species) had β = +0.0052, SE = 0.0057, t = 0.90, p = 0.366; and log(genome size) had β = +0.0255, SE = 0.0069, t = 3.68, p = 2.6×10⁻⁴. The model had λ = 0.628, and the metal type effect had BH-FDR q = 0.003. ₁

In the leave-one-metal-out analysis, the number of metal types was recomputed at the species level excluding each of Hg, As, Cu, Zn, Cd, Cr, and Ni in turn, aggregated to genus, and refit by PGLS. Compared with the original all-7 analysis (β = +0.0215, SE = 0.0056, p = 1.5×10⁻⁴), all 7 exclusions retained a positive coefficient but none remained significant: excluding Hg gave β = +0.0082, p = 0.233; As β = +0.0039, p = 0.518; Cu β = +0.0102, p = 0.115; Zn β = +0.0095, p = 0.140; Cd β = +0.0084, p = 0.186; Cr β = +0.0088, p = 0.173; and Ni β = +0.0095, p = 0.140. The report treats this as consistent with a distributed signal, while noting that reducing the maximum diversity score from 7 to 6 also compresses predictor variance and reduces power. All 7 exclusions drop to p > 0.10. Excluding As produced the weakest residual signal (β = +0.004, p = 0.52) and excluding Cu the strongest (β = +0.010, p = 0.12), so no single metal is singled out as the driver. The report adds an interpretive caveat: the loss of significance is consistent either with a genuinely distributed signal, or with an original finding near the significance threshold that the compressed predictor range (maximum diversity 6 rather than 7) pushes below it. The consistently positive direction supports the first reading, but the report states that the alternatives cannot be definitively distinguished without a larger metal type vocabulary. In its annotation caveats, the report also uses S1 to argue that even if one metal category had elevated false positives, the result requires contributions from multiple metal types. On this argument, a single category’s spurious annotations would not explain it. That argument is stronger than the S1 results support, because all seven exclusions are non-significant (p > 0.10). Consistently positive coefficients cannot distinguish a distributed signal from loss of power as the predictor range contracts from seven to six metals. ₁

Leave-one-environment-out analyses remained nominally significant for all 13 environment categories. The weakest result followed exclusion of aquatic environments, with β = +0.0085, SE = 0.0039, p = 0.031 and λ = 0.880; the strongest was after excluding forest, with β = +0.0142, SE = 0.0042, p = 0.0007 and λ = 0.874. The remaining coefficients ranged from β = +0.0105 to +0.0131, with p-values from 0.0025 to 0.0077, except for desert at β = +0.0120, p = 0.0040 and field at β = +0.0118, p = 0.0030. Each model excluded one category from the B_std computation, row-normalised prevalences over the remaining 12 environments, recomputed B_std with J=12, and then reran PGLS. Excluding aquatic, the most heterogeneous and most over-represented category (40,353 OTUs; 40.8% of the dataset), reduced β from +0.021 to +0.0085 (−59%). Excluding the 3 least-sampled environments, flower (898 OTUs, 0.9%), leaf (1,298 OTUs, 1.3%) and plant (2,647 OTUs, 2.7%), gave β = +0.012, +0.013 and +0.010 respectively, all p < 0.008; these values are reported at lower precision than the coefficients above. λ was consistently higher in these models (0.87–0.91 vs 0.71 in the main analysis). The report attributes this to rescaling B_std over fewer categories, which increases the relative weight of shared environments and amplifies phylogenetic covariation. After FDR correction, 12 of 13 leave-one-environment-out models remained significant at q < 0.05, while the aquatic exclusion had q = 0.061. ₁

Within-genus variation did not explain the association: adding the within-genus standard deviation of metal types produced β = +0.0189, SE = 0.0060, p = 0.0016 for mean metal types and β = +0.0074, SE = 0.0055, p = 0.181 for the standard deviation. The report clarifies that the earlier coverage analysis (R1) had already added n_species_with_amr as a PGLS covariate using gls() with corPagel, not only in OLS. This sensitivity analysis (S4) instead asks whether log-transforming the covariate improves its performance. Adding log(n_species) produced β = +0.0195, SE = 0.0057, p = 6.4×10⁻⁴ for metal types and β = +0.0127, SE = 0.0053, p = 0.016 for log(n_species). Log(n_species) was thus significant (p = 0.016), whereas the linear covariate in R1 was borderline (p = 0.07). The report reads this as a log-linear relationship in which well-sampled genera tend to be slightly broader-niched, possibly because they include more ecologically diverse isolates. The metal type effect was not explained away (β drops from +0.021 to +0.020), and the report states that metal types retain a ΔAIC of approximately −10 beyond the log(n_species)-only model. It notes that a partial R² decomposition would require a likelihood-ratio test. In the limitations section, the report specifies that log(n_species) was z-scored and calls the drop in the metal-type β from +0.021 to +0.020 only a small confound. It states that the genome-count and metal-type effects are independent. ₁

Filtering out OTUs with mean relative abundance below 0.01% removed 66% of OTUs, leaving 33,676 OTUs; 23,759 bacterial, non-organellar OTUs covered 576 of the 606 PGLS genera. Pagel’s λ for B_std changed from 0.763 to 0.750 within the same 576-genus subset, a −1.7% change, and remained highly significant with p = 2.8×10⁻⁴⁰. Per-OTU mean relative abundance was estimated from total counts, the number of detected samples, and a mean library size of 87,595 reads/sample across 463,972 samples. Genus-level B_std was then recomputed on the filtered OTU set. For comparison, unfiltered λ was 0.787 across all 1,264 bacterial genera (p = 7.9×10⁻¹⁰²) and 0.763 for the same 576 genera (p = 5.9×10⁻⁴⁰). The report concludes that the phylogenetic conservation of niche breadth is not an artefact of rare transient OTU detections. This filter tested niche-breadth conservation (Pagel’s λ), not the metal-type PGLS association. ₁

Archaeal analysis and multiple testing

The archaeal analysis included 48 genera and was severely underpowered: metal type diversity had β = +0.0145, SE = 0.0198, t = 0.73, p = 0.467, and λ = 0.567. Formal power was 11% at α = 0.05, and the report estimates that 80% power would require n ≥ 702 at α = 0.05 or n ≥ 1,084 at Bonferroni α = 0.0083; the non-significant result is therefore not interpreted as evidence against the association. Only 48 archaeal genera had both MicrobeAtlas representation (≥ 3 OTUs) and AMR data in kbase_ke_pangenome. The archaeal coefficient had the expected positive direction, but its standard errors were >1× the coefficient. The report reads this as severe underpowering rather than a true null effect. It includes the archaeal PGLS only to quantify the statistical gap between current data and what would be needed, not as evidence for or against the hypothesis. The archaeal fits (scripts/pgls_robustness.R, n = 48 genera with AMR data) found no significant effects for any predictor. AMR clusters had β = +0.0146, SE = 0.0209, t = 0.70, p = 0.488, λ = 0.567, and core fraction had β = +0.0002, SE = 0.0211, t = 0.01, p = 0.991, λ = 0.553. The report notes that the archaeal metal-type β (+0.0145) has the same positive direction as the bacterial estimate (+0.0215). It describes the SE as nearly as large as the estimate (t = 0.73), indicating severe underpowering. It also notes that the 48 archaeal genera are biased toward cultured lineages (methanogens, halophiles, thermoacidophiles) and miss the environmentally dominant MicrobeAtlas archaea Thaumarchaeota and Woesearchaeota. The power analysis used a non-central t-distribution with β_true = +0.0145, λ fixed at 0.567, and SE scaling as SE_obs × √(n_obs/n). Projected power at α = 0.05 and at Bonferroni α = 0.0083 was 11% and 3% at n = 48 (current), 18% and 6% at n = 100, 32% and 12% at n = 200, 45% and 21% at n = 300, 56% and 30% at n = 400, 74% and 48% at n = 600, 85% and 64% at n = 800, and 92% and 76% at n = 1,000. These projections are conditional on the observed archaeal effect size and error. They are not measurements. ₁

The confirmatory analysis comprised 6 simple PGLS models with Bonferroni threshold p < 0.0083. Across approximately 47 tests, the primary result p = 1.5×10⁻⁴ also survived the more conservative full-family threshold 0.05/47 ≈ 0.0011. Benjamini-Hochberg false discovery rate correction identified 23 of 47 tests at q < 0.05; the primary PGLS had q = 0.0028, the multi-predictor model q = 0.0037, the genome-size and three-covariate results q = 0.0028, and the archaeal result q = 0.55, so the archaeal PGLS remained non-significant. Further bacterial metal type effects surviving FDR correction were the R1 covariate models (q < 0.004), the within-genus SD model (q = 0.007) and the log(n_species) model (q = 0.0037). The leave-one-environment-out models were significant at q < 0.05 for 12 of 13 environments; the aquatic exclusion had q = 0.061 despite nominal significance at p = 0.031. The leave-one-metal-out models did not survive FDR correction. The report’s test inventory is listed below. ₁

• 6 simple confirmatory PGLS models, with Bonferroni correction (α/6 = 0.0083). ₁
• 6 confirmatory standalone Pagel’s λ tests, using an LRT with no additional correction. ₁
• 3 predictors in 1 multi-predictor model, not corrected beyond the model. ₁
• Exploratory-only sets of 9 robustness coefficients, 3 archaeal coefficients, 24 sensitivity coefficients and 1 OTU-abundance-filter λ test. ₁

The 47-test Bonferroni threshold was not the stated correction method and is offered only as an upper bound. The report names the six simple PGLS models (NB05) as the primary confirmatory analysis, but it also lists the standalone Pagel’s λ analyses (NB04) as confirmatory. Its statement that all other analyses are exploratory is therefore overbroad; the analyses it explicitly labels exploratory are the multi-predictor, robustness, sensitivity and archaeal analyses. ₁

The report’s summary table gives the metal-type β and the Benjamini–Hochberg false discovery rate (BH-FDR) q-value for every sensitivity and robustness analysis. The full FDR table is data/all_tests_fdr.csv. The table entries are:

• Simple PGLS (B_std ~ metal_types_z): n = 606, β = +0.0215, SE = 0.0056, p = 1.5×10⁻⁴, q = 0.003.
• Multi-predictor PGLS: n = 606, β = +0.0231, SE = 0.0067, p = 5.5×10⁻⁴, q = 0.004.
• R1 with the n_species covariate: simple model n = 606, β = +0.0204, SE = 0.0057, p = 3.4×10⁻⁴, q = 0.003; multi-predictor model β = +0.0224, SE = 0.0067, p = 8.3×10⁻⁴, q = 0.004.
• R2, rarefied median over 200 iterations: n = 606, β = +0.0147, SE ~0.011, p = 0.005, with no q-value or significance decision reported.
• R3, archaeal PGLS: n = 48, β = +0.0145, SE = 0.0198, p = 0.467, q = 0.548, not significant.
• R4, strict prevalence (≥5%): n = 379, β = +0.0166, SE = 0.0099, p = 0.092, with no q-value reported.
• R5 with the genome size covariate: n = 527, β = +0.0218, SE = 0.0061, p = 3.6×10⁻⁴, q = 0.003.
• R6 with log_n_species and log_genome_size: n = 527, β = +0.0218, SE = 0.0061, p = 3.5×10⁻⁴, q = 0.003. ₁

The table reports some values at a different precision from the text. Its multi-predictor β is +0.0231, against +0.023 in the text. Its primary and multi-predictor q-values are 0.003 and 0.004, against 0.0028 and 0.0037 above. Its archaeal q is 0.548, against 0.55 in the text. Its R2 row gives p = 0.005, alongside the median p-value of 0.0054 reported above. ₁

In the same table, every leave-one-metal-out (S1) exclusion at n = 606 kept a positive coefficient but none was BH-FDR significant. Excluding Hg gave β = +0.0082, SE = 0.0069, p = 0.233, q = 0.313. Excluding As gave β = +0.0039, SE = 0.0061, p = 0.518, q = 0.572. Excluding Cu gave β = +0.0102, SE = 0.0065, p = 0.115, q = 0.187. Excluding Zn gave β = +0.0095, SE = 0.0064, p = 0.140, q = 0.212. Excluding Cd gave β = +0.0084, SE = 0.0064, p = 0.186, q = 0.258. Excluding Cr gave β = +0.0088, SE = 0.0064, p = 0.173, q = 0.255. Excluding Ni gave β = +0.0095, SE = 0.0064, p = 0.140, q = 0.212. This limits any claim that the association is robust to removing individual metals. ₁

All leave-one-environment-out (S2) models used n = 606. Excluding aquatic gave β = +0.0085, SE = 0.0039, p = 0.031, q = 0.061. The report footnotes this result as nominally significant at p = 0.031 but not significant after FDR adjustment (q = 0.061). The other 12 exclusions were all FDR-significant:

• agricultural: β = +0.0107, SE = 0.0040, p = 0.0073, q = 0.016
• desert: β = +0.0120, SE = 0.0042, p = 0.0040, q = 0.011
• farm: β = +0.0115, SE = 0.0042, p = 0.0061, q = 0.015
• field: β = +0.0118, SE = 0.0040, p = 0.0030, q = 0.010
• flower: β = +0.0118, SE = 0.0044, p = 0.0077, q = 0.017
• forest: β = +0.0142, SE = 0.0042, p = 0.0007, q = 0.004
• leaf: β = +0.0130, SE = 0.0043, p = 0.0025, q = 0.010
• paddy: β = +0.0123, SE = 0.0042, p = 0.0036, q = 0.010
• peatland: β = +0.0131, SE = 0.0044, p = 0.0027, q = 0.010
• plant: β = +0.0105, SE = 0.0039, p = 0.0069, q = 0.016
• shrub: β = +0.0126, SE = 0.0043, p = 0.0034, q = 0.010
• soil: β = +0.0113, SE = 0.0039, p = 0.0035, q = 0.010 ₁

Two further sensitivity models at n = 606 were FDR-significant. Adding the within-genus SD of metal types (S3) gave β = +0.0189, SE = 0.0060, p = 0.0016, q = 0.007. Adding log(n_species) (S4) gave β = +0.0195, SE = 0.0057, p = 6.4×10⁻⁴, q = 0.004. The OTU ≥0.01% relative-abundance filter (S5) is reported only as Pagel’s λ = 0.750, n = 576, p = 2.8×10⁻⁴⁰, with no SE or q-value. ₁

Independent Environmental Validation

Track A analyzed 1,624 groundwater samples from arkinlab_microbeatlas (OTU × sample counts) and 767 genera. The analysis subset was restricted to genera with (i) ≥1 OTU detected in groundwater and (ii) non-null metal AMR data (mean_n_metal_types) in genus_trait_table.csv. Metal type diversity correlated positively with groundwater prevalence, with Spearman ρ = +0.112 and p = 0.0019; top-Q4 versus bottom-Q1 metal-diversity genera (213 versus 429 genera, respectively) also differed in groundwater prevalence in a one-sided (greater) Mann–Whitney test, p = 0.007. Median prevalence was 0.81% for top-Q4 genera and 0.62% for bottom-Q1 genera, a ratio of 1.31×. The report calls this effect size modest. It reads the modest size as consistent with the global PGLS β being a partial correlation after phylogenetic structure is removed. ₁

Genus-level groundwater prevalence was defined as the mean fraction of the 1,624 groundwater samples in which ≥1 OTU from that genus was detected (count > 0). Fold enrichment was defined as groundwater prevalence % divided by non-groundwater prevalence % across 462,348 samples. Among genera with AMR data, the quartile thresholds were Q1 ≤ 1.0 metal types and Q4 ≥ 1.5 metal types. ₁

The groundwater-specific fold-enrichment analysis was null. Across 767 genera, Spearman ρ = +0.042 and p = 0.242. The Mann–Whitney comparison of fold enrichment between top-Q4 and bottom-Q1 metal-diversity genera (213 versus 429 genera) had p = 0.181. The report therefore treats Track A as supporting a prevalence association but not demonstrating groundwater-specific enrichment or uniquely confirming a metal-contamination mechanism. The report calls this a meaningful null. The environmental-specificity metric asks whether metal-diverse genera are preferentially enriched in groundwater beyond their global baseline, and it shows no detectable signal. Track A thus confirms a modest prevalence difference (0.62% → 0.81% median prevalence for bottom-Q1 versus top-Q4) but not groundwater-specific enrichment. The report offers two explanations for the null and does not distinguish between them. Either (a) the prevalence advantage is a general breadth effect rather than groundwater specificity, or (b) the 1,624-sample groundwater subset is underpowered to detect a specificity signal of the observed effect size. The report nonetheless describes Track A as an independent, field-observation-based cross-validation of the PGLS finding, because it uses a different data source (groundwater surveys versus a multi-environment atlas). It calls the result consistent with the idea that diverse metal resistance is associated with greater ecological range in chemically challenging subsurface environments. Given the null specificity result, this interpretation does not establish the proposed mechanism. ₁

The report’s candidate-genera table gives three values for each genus: groundwater prevalence in percent, fold enrichment and metal-type diversity. The values are Azospira 5.21, 2.3×, 3.0; Enterobacter 4.80, 1.1×, 3.1; Citrobacter 0.89, 2.8×, 3.7; Enterococcus 0.71, 2.8×, 3.3; Noviherbaspirillum 1.01, 1.4×, 3.0; and Franconibacter 0.39, 1.6×, 3.5. Pseudomonas, with a metal-type diversity of 2.1, was not in the top-15. The report notes that Citrobacter and Enterococcus are notably enriched in groundwater (2.8× fold enrichment), which it calls consistent with the ENIGMA experimental predictions. It describes Azospira as a denitrifier detected in 5.2% of groundwater samples, with the highest prevalence among high-metal-diversity genera and an established role in groundwater nitrogen cycling. These per-genus enrichments are descriptive only, because the genus-wide fold-enrichment association was not significant. ₁

Track B processed all 133 PRJNA1084851 samples into a 133-sample × 24,295-OTU table, with 14,021 OTUs assigned genus-level taxonomy and 1,637 OTUs representing 16.8% of reads joinable to genus metal-diversity data. Community-weighted mean metal-type diversity ranged from 1.01 to 1.83, with mean 1.28 ± 0.19. The samples come from a Water Research 2024 study (doi:10.1016/j.watres.2024.121460), “Reproducible responses of geochemical and microbial successional patterns in the subsurface to carbon source amendment”. That study deposited 133 MiSeq 16S amplicon runs under PRJNA1084851 (SRR28246018–SRR28246150). The report’s bibliography describes PRJNA1084851 as 133 MiSeq 16S samples from the ENIGMA-funded ORFRC in Tennessee. The pipeline used cutadapt 5.2 for primer trimming, vsearch 2.30.6 for paired-end merging and 97% OTU clustering, and SILVA 138 for taxonomy. The final table had a median of 56,908 reads/sample. CWM was computed for each sample as the abundance-weighted mean of genus-level mean_n_metal_types across OTUs with AMR data, which covered 16.8% of reads per sample on average. ₁

PRJNA1084851 is the ENIGMA ORFRC 16S validation dataset, and its full-pipeline processing is reported as complete: 133/133 samples processed and 24,295 OTUs. The resulting data/enigma_otu_table.csv (133 samples × 24,295 OTUs) used 97% OTU clustering and SILVA 138 taxonomy. Raw FASTQs were deleted after processing. The end-to-end script scripts/prjna1084851_pipeline.py runs ENA FTP download → cutadapt primer trim → vsearch merge → 97% OTU clustering → SILVA 138 taxonomy → OTU table, and it processed 81 samples via ENA. A separate recovery script, scripts/prjna1084851_recover.py, handled 52 PRJNA1084851 samples that were not mirrored to ENA FTP. It used NCBI fasterq-dump v3.4.1, cutadapt and vsearch merging to produce merged FASTQs for the main pipeline. ₁

Community-weighted mean metal-type diversity differed across 8 wells with Kruskal–Wallis H = 29.10 and p = 0.0001 (133 samples), increased with time point with Spearman ρ = +0.383 and p < 0.0001 across 133 samples, and increased within FW216 with ρ = +0.576 and p = 0.0001 across 40 samples. EVO09 versus EVO17 (52 versus 81 samples) differed with Mann–Whitney p = 0.0012; the U statistic is not reported. FW215 and FW216 had the highest median values, 1.37 and 1.36, respectively, compared with 1.23 for GP01, 1.22 for GP03, 1.14 for MLSB3, and 1.09, 1.08, and 1.08 for MLSA3, FW202, and MLSG4. According to BioSample attributes, the design comprised two experiments, EVO09 (52 samples, Jan 2009) and EVO17 (81 samples, 2017–2018). Together they span 8 groundwater monitoring wells (FW202, FW215, FW216, GP01, GP03, MLSA3, MLSB3, MLSG4) at the ENIGMA ORFRC, Oak Ridge, Tennessee. The report does not state which wells each experiment sampled. All samples were collected at 0.2 µm filtration. Replicate codes R0–R20 represent successive time points before and after ethanol carbon amendment. ₁

Wells FW215 and FW216 are described as lying within the U/NO₃ contamination plume at ORFRC (oak-ridge-field-research-center). The report reads their higher CWM as consistent with the core hypothesis that metal-contaminated subsurface environments select for genera with broader metal resistance repertoires. The pattern is consistent with such selection but does not establish it. The report also suggests that carbon amendment progressively selects for metal-diverse communities, with the effect strongest in FW216 (ρ = +0.576, p = 0.0001). It proposes that this may happen because amendment creates reducing conditions that shift the geochemical speciation of redox-active metals (U, Fe, S). This mechanism is a hypothesis, not a demonstrated cause. ₁

Sulfurimonas, a sulfur-oxidizing chemolithotroph with 80 OTUs and 96.2% prevalence, was the dominant flagged genus in the same Track B dataset. Its relative abundance nonetheless decreased with time (Spearman ρ = −0.323, n = 133, p = 0.0002). The report says this decline is likely because carbon amendment depletes nitrate and oxygen, which are Sulfurimonas’ preferred electron acceptors. The report interprets this as a counterexample showing that ecological dominance can be driven by electron-acceptor availability rather than metal resistance. It shows that not all ecologically dominant groundwater taxa are metal-diverse, although the electron-acceptor explanation remains an interpretation. Nitrospira (25 OTUs, 88.7% prevalence) had its highest abundance in the FW215 plume well. Candidatus Jorgensenbacteria (Patescibacteria/CPR; 122 OTUs, 99.2% prevalence) was ubiquitous. Citrobacter and Thermodesulfovibrio were not detected. The report attributes this either to abundance below detection or to genus names mismatched under SILVA nomenclature, and does not distinguish between the two. Overall, the report interprets Track B as supporting the core PGLS finding at the within-study level, with Sulfurimonas as an informative counterexample. In its reading, contaminated plume wells (FW215, FW216) host communities with measurably higher metal-type diversity than background wells. It also reads the temporal CWM increase after carbon amendment as pointing to active ecological selection for metal-diverse communities under geochemical stress. That selection reading is interpretive. ₁

Experimental Validation Framework

The study identified 435 candidate OTUs: 215 in the top 10% for both Levins’ B_std and mean metal-type diversity, using B_std ≥ 0.718 and n_metal_types ≥ 2.0, plus 220 nitrifier positive controls. Eight priority OTUs were selected for microcosm experiments, with Levins’ B values from 0.770 to 0.999 and metal-type values from 2.1 to 3.7; the recommended stresses included Ag, Te, Cu, Cd, As, Hg, and As + Hg co-stress. The 435 candidates were selected from 98,919 MicrobeAtlas OTUs. They are listed in data/candidate_otu_list.csv, which carries a top10pct_both flag and a nitrifier_role field. The 220 nitrifier controls come from Nitrosomonas, Nitrosopumilus, Nitrososphaera, Nitrospira, Nitrobacter, Nitrococcus, Nitrosospira and Nitrospina. For these genera the metal AMR–niche breadth link is not expected to apply, because nitrification is phylogenetically entrenched (λ = 0.94). A separate file, data/honorable_mentions.md, lists 2,186 genera that have metal AMR gene annotations but zero candidate OTUs, because none were in the top 10% by both metrics. ₁

Sample counts confirm that the candidate OTUs span 13 environment types globally; the full list is in data/candidate_otu_list.csv. Eight OTUs were prioritised for experimental follow-up by a composite score (Levins’ B_std × n_metal_types), each representing a distinct metal-type signature. The shortlist is:

• Rank 1: OTU 97_56843 (Klebsiella), Levins’ B 0.947, 3.5 metal types, Ag or Te stress
• Rank 2: OTU 97_51587 (Enterococcus), Levins’ B 0.948, 3.3 metal types, Cu or Cd stress
• Rank 3: OTU 97_14443 (Citrobacter), Levins’ B 0.854, 3.7 metal types, As or Hg stress
• Rank 4: OTU 97_3668 (Franconibacter), Levins’ B 0.785, 3.5 metal types, Ag or As stress
• Rank 5: OTU 97_12431 (Noviherbaspirillum), Levins’ B 0.770, 3.0 metal types, As stress
• Rank 6: OTU 97_66129 (Serratia), Levins’ B 0.923, 2.4 metal types, Hg stress
• Rank 7: OTU 97_80244 (Aeromonas), Levins’ B 0.981, 2.2 metal types, As + Hg co-stress
• Rank 8: OTU 97_39338 (Pseudomonas), Levins’ B 0.999, 2.1 metal types, Cu stress ₁

Citrobacter tops metal-type diversity (3.7 types) and is linked to silP/A, arsC/D, chrA and merA/P. The report describes Enterococcus (Ag, As, Cd, Cu, Hg, Te) and Klebsiella (Ag, As, Hg, Te) as the broadest multi-metal resisters. Pseudomonas ranks eighth by composite score but has the highest niche breadth of all candidate genera (B_std = 0.999). It also carries a complete copA/pco Cu-resistance operon. ₁

Each of the 8 shortlisted OTUs has a specific, falsifiable prediction in data/hypotheses_refined.md. The general form is that under a specified metal stress the OTU should increase in relative abundance by ≥2-fold versus metal-free controls after 7–14 days. The stated rationale is that the OTU carries resistance determinants that neutralise metal toxicity while susceptible competitors are inhibited. These predictions are untested. The report also gives falsification criteria. If an OTU does not increase, the candidate explanations are that (a) its resistance genes are on accessory elements absent in the soil ecotype, or (b) the metal concentration exceeds the allele-specific minimum inhibitory concentration (MIC). ₁

The proposed experiment uses 4–6 metal-stress treatments, 5 replicates, and 7- and 14-day time points, with relative-abundance fold-change versus metal-free controls as the primary endpoint. The recommended pre-treatment sequencing design uses 5 replicates and ≥ 50,000 reads per sample to confirm detectable OTUs at ≥ 0.1% relative abundance; at 1% initial abundance, the report estimates coefficient of variation ≈ 10% and power >99% for detecting a 2-fold increase with 3 replicates at 10,000 reads. The detection limit is ≥ 0.1% relative abundance (≥ 10 reads at 10,000 reads/sample). The report states that the binding constraint is multiple testing, not sensitivity. The pre-treatment sequencing is of the ENIGMA inoculum. The experiment would then be restricted to the 3–5 confirmed OTUs, with false discovery rate (FDR, the expected share of false positives among discoveries) correction applied within that set only. These power figures are design estimates, not results. ₁

If OTU-level tracking fails, the backup is genus-specific 16S qPCR (Pseudomonas, Citrobacter, Enterococcus) or a defined synthetic community (SynCom) of 8–12 priority taxa in place of amplicon sequencing. ₁

The report cites three external studies to support the design. Farkas et al. (2023, Geomicrobiology Journal) tracked prokaryotic community composition by 16S rRNA amplicon sequencing in As(III) microcosm experiments and found acute composition shifts detectable at community level within 7 days. The report treats this as consistent with its 7–14-day design. Kou et al. (2018, Frontiers in Microbiology) showed selective OTU-level enrichment under Pb/Zn/Cu in a field trial, which the report takes as support for the premise that metal-tolerant taxa reliably increase in relative abundance. Yuan et al. (2020, Current Microbiology) found that copper specifically reduces OTU richness in freshwater microcosms. From this, the report proposes, without testing, that reduced competition would amplify the selective advantage of Cu-resistant taxa such as Pseudomonas OTU 97_39338. ₁

Caveats and Limitations

MicrobeAtlas niche breadth is a sequencing-effort proxy rather than confirmed ecological range. Detection across 13 environment categories is affected by sampling intensity, primer bias, geographic and temporal non-uniformity, and missing environments; Levins’ B_std corrects for unequal environment sizes but cannot correct for missing environments or detection-probability differences. Genus aggregation also means that a broad genus score can reflect different species occupying different habitats rather than one organism being a generalist. The report specifies that the primer bias involves V3–V4 or V4 regions and that the non-uniformity concerns SRA/ENA sample contributions. Because of genus aggregation, the PGLS β reflects covariation between inferred niche breadth and inferred metal tolerance capacity, not confirmed physiology. ₁

The strict 5% prevalence analysis preserved the direction of the association but reduced the PGLS sample size by 37%, from 606 to 379, and produced p = 0.092. Primer bias, missing environments, and the ecological meaning of sparse detections remain unresolved without a multi-primer, multi-region survey design. The R4 strict-filter coefficient was β = +0.0166 (p = 0.092, n = 379), compared with the original β = +0.0215, p = 1.5×10⁻⁴ (n = 606). In this limitations section, the report rounds these coefficients to β = +0.017 versus +0.022 original. It says the filter confirms the signal direction but reduces significance because 37% fewer genera pass, and concludes that the effect is not explained by rare transient detections. Direction alone does not establish robustness under the filter, because the filtered result was non-significant (p = 0.092). ₁

The 48-genus archaeal analysis is biased toward cultured methanogens, halophiles, and thermoacidophiles and misses environmentally dominant groups including Thaumarchaeota and Woesearchaeota. The report states that ≥200 archaeal genera with AMR data are needed for adequate power. Its own formal power analysis, however, estimates n ≥ 702 for 80% power at α = 0.05 under the observed effect and error assumptions, so the ≥200 adequacy figure is in tension with the formal requirement. The report’s own power table projects only 32% power at α = 0.05 for n = 200. Reaching ≥200 archaeal genera with AMR data would require expanded archaeal pangenome sequencing and annotation in kbase_ke_pangenome, which the report identifies as new data. ₁

Pangenome coverage is uneven: the Spearman correlation between sequenced genomes per genus and inferred metal type diversity was r = 0.35 with p = 2.6×10⁻¹⁹, and single-genome genera had lower metal type counts than multi-genome genera with Mann–Whitney p = 5.1×10⁻¹⁰. Coverage and rarefaction analyses substantially addressed this concern, but a formal rarefaction curve across multiple species-per-genus thresholds remains outstanding. ₁

Environment categories are heterogeneous, especially aquatic, which combines marine, freshwater, estuarine, hydrothermal-vent, and pond environments. Aquatic represented 40,353 OTUs, or 40.8% of the dataset; excluding it reduced β from +0.021 to +0.0085 but retained nominal significance at p = 0.031. A taxon detected in both a marine coastal sample and a freshwater lake sample counts as broad-niched, although this could reflect closely related but ecologically distinct ecotypes rather than a single generalist population. Excluding any environment other than aquatic gave p < 0.01. The report states that the result is not contingent on the heterogeneous aquatic category. That statement rests on nominal significance, because the aquatic exclusion does not survive false discovery rate (FDR) correction (q = 0.061), while 12 of 13 exclusions survive at q < 0.05. Finer subdivision into marine, freshwater, and saline categories requires re-analysis of MicrobeAtlas sample metadata, which is not available in the data extraction used for this project. Because MicrobeAtlas does not provide a finer environment classification, the analysis cannot test sub-category heterogeneity such as marine versus freshwater within aquatic. Splitting categories would require re-running the niche breadth computation with a revised environment taxonomy. The report’s quantitative bias note calls the over-representation of aquatic environments a potential source of downward bias when the category is included. Yet the weaker β after excluding aquatic (+0.0085 vs +0.021) suggests, without establishing, that aquatic OTUs contribute disproportionately to the association; the two readings are in apparent tension. ₁

The Gaussian PGLS treats metal type diversity, an integer count from 1–7, as a continuous z-scored predictor; the report considers this standard when the count is an independent variable but notes that the Gaussian assumption applies to the residuals of B_std, not to the distribution of metal types. The report says inspecting PGLS residuals for normality is warranted but does not report that this check was performed. It adds that the main concern would arise only if metal types were modelled as a response (e.g., what predicts metal type count?), where a Poisson or negative binomial PGLS (phylolm with Poisson family, or MCMCglmm) would be appropriate; these are proposed for that reverse regression, not as replacements for the current B_std-response PGLS. Pagel’s λ for the count trait also assumes Brownian-motion evolution (Gaussian increments), an approximation for an integer trait. The report names discrete-state (Mk) or threshold models as alternatives for this λ analysis, while noting they are not standard for this type of continuous-trait PGLS analysis. ₁

The GTDB r214 genus-representative tree, constructed from 120 concatenated marker genes, was used for all phylogenetic analyses, but sensitivity to NCBI, SILVA-based 16S, or other GTDB trees was not tested. Genus-level pruning discards within-genus phylogenetic structure, although the within-genus metal-type standard-deviation test was non-significant at p = 0.18. ₁

AMRFinderPlus annotations use hidden Markov model (HMM)-based detection against the NCBI Bacterial AMR Reference Gene Database. For metal resistance, AMRFinderPlus uses curated HMMs built from experimentally validated reference sequences (e.g., merA/merB for mercury, arsABCDR for arsenate, copABCD for copper). HMM gathering thresholds are set during model calibration, and no additional identity threshold is applied because HMM scores are not identity values. Cross-reactive metalloproteases or ABC transporters with structural similarity to metal resistance proteins could be misclassified, although AMRFinderPlus applies taxon-specific cutoffs to reduce this. The “other” metal category (17% of annotated clusters) captures metal-associated genes not classified into the 7 primary types and may include weaker hits. The report says excluding it would produce a 6-metal analysis yet lists seven metals (Hg, As, Cu, Zn, Cd, Cr, Ni), so that description is internally inconsistent. No subset was manually validated. Spot checks of highly prevalent gene families (arsD, merP) showed biologically expected distributions, high in contaminated-environment genera and lower in strict aerobes, which the report treats as only informal support. The recommended follow-up is manual BLAST verification of 100 randomly sampled annotated gene clusters. ₁

Track B community-weighted means cover only 16.8% of reads because the remaining approximately 83% could not be joined to genus-level GTDB AMR annotations. If uncovered genera have systematically different resistance profiles, the community-weighted means could be biased upward; correlation between per-sample covered-read fraction and community-weighted mean is required as a diagnostic. The 16.8% reflects OTUs whose genus-level SILVA taxonomy could be matched to GTDB AMR annotations. The report argues that the between-sample CWM variance (1.01–1.83) is consistent with a real biological signal. It makes the same argument for the correlation with known contamination status (FW215/FW216 in the plume). It states, however, that robustness to variation in coverage fraction has not been formally tested. ₁

The study’s robustness and sensitivity analyses were post-hoc exploratory analyses rather than confirmatory tests, and the analysis plan was not preregistered. The report says the robustness analyses were conducted after the main analysis in response to anticipated reviewer concerns, and should be read as reassurance rather than confirmatory evidence. It describes the primary metal-types-to-B_std result as “Bonferroni-corrected” with p = 1.5×10⁻⁴. That wording does not establish whether 1.5×10⁻⁴ is an adjusted p-value or an unadjusted p-value assessed against the Bonferroni threshold. The report also identifies unresolved needs for HGT-burden predictors, site-level contamination metadata, co-resistance clustering, finer environment classification, alternative phylogenetic mixed models, AMRFinderPlus validation, and a public Zenodo archive with a citable DOI. ₁

The future-directions list does not fully match the completed analyses. Future Direction #1 proposes pangenome rarefaction: resampling genera to equal genome coverage (e.g., n = 1 or 2 species) and repeating PGLS to formally control for the genome-sampling depth confound. It says this would require access to the species-level AMR table stratified by genome count. The report does not reconcile this proposal with the R2 one-species-per-genus rarefaction reported above, so the proposal should not be read as evidence that no rarefaction was performed. Future Direction #2 proposes expanding archaeal pangenome coverage in kbase_ke_pangenome to environmental metagenome-assembled genomes (MAGs, genomes reconstructed from metagenomes), particularly Thaumarchaeota. It would repeat the archaeal PGLS once ≥ 100 archaeal genera with AMR data are available. This ≥ 100 threshold differs from both the report’s ≥200 adequacy figure and its formal n ≥ 702 power requirement. Future Direction #3, the genome size covariate, is marked as completed as Robustness Analysis R5. Genus-level mean genome size from kbase_ke_pangenome.gtdb_metadata was added as a PGLS covariate, and the metal type diversity effect is described as independent of genome size (ΔAIC = 10.8). ₁

Proposed follow-up analyses and review responses

• Future Direction #4 proposes genus-level pangenome openness, the proportion of genes in the accessory genome, as an HGT-burden predictor of niche breadth alongside metal type diversity. The rationale is that a broader accessory genome may be the underlying driver of both traits. This was not tested. ₁
• Future Direction #5 proposes site-level validation. It would test whether MicrobeAtlas genera with high metal type diversity are enriched in metal-contaminated sample metadata (where available: mining-impacted, biosolid-amended, estuarine) versus clean reference environments. The report says this would directly validate the environmental co-selection hypothesis. It remains proposed. ₁
• Future Direction #6 proposes grouping metal types by genetic co-occurrence (e.g., Cu-Zn co-resistance vs Hg-As). It would then test whether specific co-resistance combinations predict niche breadth more than the raw count. ₁
• Future Direction #7 proposes MCMCglmm with a Poisson or negative binomial family to test the reverse regression direction: whether metal type count, as a response variable, is predicted by niche breadth with phylogenetic structure accounted for. It was not performed. ₁
• Future Direction #8 proposes re-running the niche breadth computation with a finer environment taxonomy built from MicrobeAtlas sample metadata, for example splitting “aquatic” into marine, freshwater and estuarine. This would test whether category heterogeneity inflates B_std for certain genera. ₁
• Future Direction #9 proposes manual BLAST verification of 100 randomly sampled annotated metal AMR gene clusters against reference sequences. This would quantify the hidden Markov model (HMM)-based false-positive rate for each metal type category. No such validation is reported. ₁
• Future Direction #11 proposes depositing notebooks, derived data CSVs and a conda environment lockfile in a Zenodo repository, so that all PGLS and Pagel’s λ results can be verified independently. The derived files (data/*.csv, figures/*.png) would be archived under a citable DOI before manuscript submission. No deposit is reported. ₁

Future Direction #10, ENIGMA field validation, is marked fully addressed, with Tracks A and B complete. The report restates Track A as 1,624 groundwater samples in which metal-diverse genera are significantly more prevalent (Spearman ρ = +0.112, p = 0.0019; Mann–Whitney p = 0.007). It restates Track B as a full amplicon pipeline run on all 133 PRJNA1084851 samples (24,295 OTUs; 133/133 QC pass). In Track B, CWM metal-type diversity is higher in plume wells FW215/FW216 (Kruskal–Wallis p = 0.0001) and increases with time after carbon amendment (Spearman ρ = +0.383, p < 0.0001 across all wells; ρ = +0.576, p = 0.0001 in FW216 alone). The report concludes that both tracks independently confirm the PGLS prediction in field data. That conclusion is stronger than the evidence supports. Track A’s groundwater-specific fold-enrichment test was null, Track B CWM covers 16.8% of reads, and neither track establishes the proposed metal-selection mechanism. ₁

The report’s reviewer-response table records the status of each critique. Environment-category validity (Q1a) is only partially addressed: leave-one-out was done (S2), but sub-category splitting requires new metadata. Primer bias and differential sequencing depth (Q1b) are not correctable with available data and are only acknowledged. On tree sensitivity (Q2b), GTDB r214 is called the best available tree. An alternative-tree comparison was not performed, which the report calls a known gap. Within-genus heterogeneity was tested via S3 and was non-significant. The λ discrepancy (Q2c) is described as directly addressed, with Δλ = 0.079 expected given the predictor’s own phylogenetic signal. On the genome-count covariate (Q3a), the response states that R1 and S4 both use gls() + corPagel (PGLS), not OLS, and that β remains significant after the PGLS covariate. On metal-type classification (Q3b), the single-metal question is called addressed because all 7 S1 exclusions were positive and no one metal is the driver. Formal AMRFinderPlus validation was not performed; only spot checks were done. All seven S1 exclusions were also non-significant, so the “addressed” status for the single-metal question rests on direction alone. ₁

On genome size (Q4a), the response states that the formal test was performed as R5. Metal types had β = +0.022 (p = 3.6×10⁻⁴) independent of genome size, with ΔAIC = 10.8. On causation and directionality (Q4b), the response states that the 3 scenarios are not distinguishable and that the cross-sectional design limitation is explicit. It adds that causal phrasing was removed from Key Findings and Interpretation. On multiple comparisons (Q5a), it reports that 47 tests were tabulated and BH-FDR was applied across all of them. The primary result survives both Bonferroni (α/47=0.0011) and FDR (q=0.003). On rarefaction and within-genus variance (Q6b), the response states that S3 confirms genus means are not inflated by outlier species and that R2 rarefaction shows a consistent positive direction. This sits uneasily with Future Direction #1’s call for pangenome rarefaction, so it is unclear which sampling-depth control the report considers outstanding. ₁

On data archiving (Q7a), the response states that a public Zenodo archive has not yet been prepared. The materials are in the institutional repository only, and Zenodo archiving was added as Future Direction #11. On PGLS ordering (Q7b), the response states that scripts/test_pgls_ordering.py generates synthetic data, scrambles the row order and verifies that the sorted model recovers the correct λ and β. The test passes. On post-hoc analyses (Q8a), the response states that all robustness and sensitivity analyses are labelled exploratory. It says the primary confirmatory result (6 PGLS + Bonferroni) stands independently. On preregistration (Q8b), the report acknowledges that the plan was not preregistered and that the confirmatory/exploratory distinction was documented post hoc. On a concrete ENIGMA prediction (Q8c), it calls the item fully addressed. It cites top-quartile metal-diversity genera as 31% more prevalent in groundwater (Mann–Whitney p=0.007). It also cites CWM metal diversity as significantly higher in contaminated FW wells (KW p=0.0001) and increasing with carbon amendment time (ρ=+0.38, p<0.0001) across 133/133 ENIGMA ORFRC samples. The microcosm predictions themselves remain untested. ₁

Data, Code and Outputs

Derived data files (data/*.csv) and analysis code (notebooks/, scripts/) are committed to the BERIL Research Observatory institutional repository. A public Zenodo archive with a citable DOI has not yet been prepared (Future Direction #11). Raw source data reside in the KBase Data Lakehouse and are accessible to KBase users. The arkinlab_microbeatlas collection (tables otu_metadata, otu_counts_long, sample_metadata) supplies 98,919 OTUs × 463,972 samples for niche breadth and taxonomy. The kbase_ke_pangenome collection (tables bakta_amr, gene_cluster, gtdb_species_clade) supplies per-species metal AMR gene annotations and GTDB taxonomy. The GTDB r214 reference tree is publicly available at https://gtdb.ecogenomic.org/. The generated OTU-to-GTDB-genus mapping, data/otu_pangenome_link.csv (NB03), has 22,357 rows. The full FDR table is data/all_tests_fdr.csv. The bibliography cites Finn et al. (2020, FEMS Microbiology Ecology 96(8): fiaa131) for MicroNiche, an R package for assessing microbial niche breadth and overlap from amplicon sequencing data. MicroNiche is an external work, not a result of this project. ₁

Slots Into

• metal-cross-resistance — The PGLS result distinguishes breadth across metal types from total AMR gene burden and supports a cross-metal resistance ecology interpretation. ₁
• environmental-resistome — The study integrates AMRFinderPlus metal-resistance annotations with global environmental occurrence data and independent groundwater validation. ₁
• environment-embedding-geography — MicrobeAtlas environment categories, groundwater prevalence, and ENIGMA well-level gradients connect resistance repertoires to environmental distributions. ₁
• gene-function-acquisition-depth — Intermediate phylogenetic signal for metal AMR, contrasted with strongly conserved niche and nitrification traits, informs the balance of vertical inheritance and horizontal acquisition. ₁
• ecotype-environment-gene-content — The association between genus-level metal-resistance diversity and inferred niche breadth provides a testable link between environmental range and gene-content variation, while retaining explicit causal and aggregation caveats. ₁
• composite-resistance-score-limitations — Total AMR cluster count and core AMR fraction fail to predict niche breadth while metal type diversity does, so the way a resistance repertoire is summarized changes the inference. ₁
• confirmatory-exploratory-ecological-association-discordance — The confirmatory metal-type association becomes non-significant under strict prevalence filtering and under leave-one-metal-out exclusions, and core fraction is only nominally associated with number of environments. The underpowered archaeal analysis is explicitly treated as evidence neither for nor against the association. ₁
• phylogenetic-confounding-of-pangenome-associations — The modest raw correlation, compared with the PGLS coefficient after shared ancestry is removed, illustrates phylogenetic control of a pangenome-trait association. ₁
• sampling-depth-and-downsampling-effects — Species-count covariates, one-species-per-genus rarefaction and strict prevalence filtering test whether sampling depth drives the association. ₁
• ontology-and-category-schema-sensitivity — Leave-one-environment-out recomputation over 13 MicrobeAtlas categories, the heterogeneous aquatic category and the lack of finer environment classes show how the environment schema shapes the niche-breadth association. ₁
• taxonomic-resolution-dependent-functional-inference — Genus-level pruning discards within-genus phylogenetic structure, tested only partially by the non-significant within-genus metal-type SD covariate. ₁
• functional-marker-validation — Metal-resistance calls rest on curated AMRFinderPlus HMMs with possible cross-reactive false positives, a weak-hit “other” category and no manual validation beyond informal arsD/merP spot checks. ₁
• chromosomal-and-integrative-gene-transfer — The lower phylogenetic signal of metal AMR is attributed to mobile-element exchange by interpretation, without direct measurement of horizontal gene transfer. ₁
• pangenome-integration — The project joins a 6,789-species pangenome database to a 98,919-OTU global atlas to test a genus-level trait association at global scope. ₁
• spatial-sampling-effort-confounding — MicrobeAtlas detections conflate habitat use with sampling intensity, primer bias and uneven geographic and temporal sample contributions. ₁
• callability-limited-comparative-inference — Levins’ B_std cannot correct for missing environments or detection-probability differences, and the underpowered archaeal PGLS is treated as uninformative, not as a null. ₁
• cultivation-collection-bias-in-ecological-genomics — The archaeal subset favours cultured methanogens, halophiles and thermoacidophiles, and primer bias and missing environments remain uncorrectable with current data. ₁
• genomic-under-representation — Only 48 archaeal genera had AMR data, and the report says ≥200 would be needed for adequate power, requiring new pangenome sequencing and annotation. ₁
• taxonomic-nomenclature-reconciliation — The analysis depends on a 22,357-row OTU-to-GTDB-genus mapping that links 16S OTUs to pangenome taxonomy. ₁
• subsurface-hydrogeological-zonation — The Track B well-level CWM gradient ties resistance repertoires to subsurface zones: plume wells FW215 and FW216 are highest, CWM rises over time within FW216, and Sulfurimonas and Nitrospira show well- and time-specific patterns. ₁
• study-batch-confounding-of-environmental-associations — The EVO09 versus EVO17 CWM difference compares sample groups within one study, so group structure must be separated from metal selection. ₁
• classifier-database-compatibility-in-taxonomic-quantification — Non-detection of Citrobacter and Thermodesulfovibrio may reflect genus-name mismatches under SILVA nomenclature rather than true absence. ₁
• adversarial-research-quality-assurance — A synthetic-data unit test verified that PGLS recovers known λ and β only when data are sorted to tree tip order. ₁
• pangenome-openness-determinants — The report proposes, but does not test, genus-level pangenome openness (accessory-genome proportion) as an alternative predictor of niche breadth. It raises the hypothesis that a broader accessory genome may drive both metal-resistance diversity and niche breadth. ₁