Community Metabolic Ecology via NMDC × Pangenome Integration

This project performs the first integration of NMDC multi-omics data with the BERDL pangenome
collection to test whether genome-predicted metabolic potential explains observed community
chemistry. For each NMDC sample with both taxonomic profiles and metabolomics measurements, we
compute a community-weighted GapMind pathway completeness score — the mean pathway
completeness across resident taxa weighted by their relative abundance — and test whether these
scores correlate with metabolomics profiles.

The central hypothesis is that Black Queen dynamics are detectable at community scale:
communities with more complete amino acid biosynthesis pathways will show lower ambient amino
acid concentrations, because those communities produce and internally consume amino acids
rather than importing them. A secondary test asks whether community metabolic potential
clusters by habitat type (soil, sediment, marine) more strongly than taxonomic composition
alone.

Data Overview

Ecosystem distribution in the full 220-sample cohort: 126 Soil, 33 Freshwater,
61 with unrecorded ecosystem type. All 33 Freshwater samples lacked paired
metabolomics data in NMDC and were therefore excluded from H1 testing; they were
retained for H2 (PCA) using pathway completeness alone.

Taxonomy bridge quality was high: mean bridge coverage 94.6% (all 220 samples passed
the 30% QC threshold), with 92% of samples mapping ≥85% of community abundance to
GTDB pangenome species.

H1: Black Queen Hypothesis

Community pathway completeness was computed using the binary frac_complete metric
(fraction of taxa in the community whose GapMind score ≥ 5, i.e., "complete" with no
missing steps). This threshold was chosen over frac_likely_complete (score ≥ 4)
because a score of 5 reflects GapMind's unambiguous "complete" call — the taxon
possesses all required pathway steps — which maps directly onto the BQH prediction that
taxa without a complete pathway are the ones expected to depend on environmental supply.
The frac_likely_complete metric is computed in NB03 for sensitivity comparisons.

Scatter plots for the four lowest-p pathways:

Sensitivity: Soil-only stratification

To control for the strong ecosystem-type effect (Finding 2), H1 was re-run on the
125 Soil-ecosystem samples only (excluding the 6 Unknown-ecosystem samples in the H1
dataset). Results: leu remains FDR-significant (r = −0.390, q = 0.022, n = 62 —
identical to the full-cohort result, since the leu samples are all Soil); arg loses
FDR significance (r = −0.264, q = 0.117, n = 78) but remains directionally consistent.
The binomial sign test (11/13 negative, p = 0.011) is unchanged, confirming the
majority-negative direction is robust to ecosystem stratification.

Sensitivity: Study-level blocking

Study-level analysis of the H1 dataset reveals that 95% of samples (125/131)
originate from a single NMDC study (nmdc:sty-11-r2h77870); the remaining 6 samples
are from a second study (nmdc:sty-11-547rwq94). This means cross-study LC-MS protocol
heterogeneity is effectively absent as a confounder for leu: all leu samples (n = 62)
are from one study with consistent analytical methods. For arg, excluding the 6-sample
second study gives r = −0.264, p = 0.019 (n = 78), confirming the arginine signal is
not driven by the minor study.

(Sensitivity analyses: 05_statistical_analysis.ipynb Part 6)

H2: Ecosystem Differentiation

Top PC1 loadings (most positive = Soil direction): glucuronate (0.146), fumarate (0.144),
succinate (0.143), lysine (0.140), cellobiose (0.139). These are all carbon-utilization
pathways with uniformly high positive loadings, indicating that Soil communities have
substantially higher predicted carbon metabolic diversity than Freshwater communities.

H1: Black Queen Signal

The results provide weak but consistent support for community-scale Black Queen
dynamics in environmental soil microbiomes. The finding that 11 of 13 tested pathways
trend in the BQH-predicted direction (p = 0.011 binomial test) is notable, and the
signal is robust in two sensitivity checks: it is unchanged when the analysis is
restricted to Soil-ecosystem samples (soil-only: 11/13 negative, binomial p = 0.011;
leu q = 0.022), and the dataset is effectively single-study — 95% of H1 samples come
from one NMDC study (nmdc:sty-11-r2h77870), so cross-study LC-MS protocol heterogeneity
is not a confounder for leu.

The two FDR-significant pathways — leucine and arginine — are both energetically
expensive to synthesise (leucine: 37 ATP equivalents; arginine: 26 ATP equivalents),
consistent with the BQH prediction that costly functions are most likely to be lost
when environmental supply is reliable. Methionine shows the largest effect size
(r = −0.50) but is underpowered (n = 18, q = 0.117); a larger metabolomics dataset
would likely reach significance for this pathway. The leucine signal was strengthened
by the isoleucine/leucine compound-mapping bug fix (r −0.326 → −0.390, q 0.045 →
0.022): three isoleucine compounds that had been misassigned to leucine had near-zero
completeness-vs-metabolite correlation, diluting the true leucine signal.

The tyrosine outlier (r = +0.42) merits caution: tyrosine can be produced from
phenylalanine via non-biosynthetic hydroxylation, so communities with high
phenylalanine biosynthesis completeness — not captured in our tyr completeness
score — may still supply tyrosine, decoupling the biosynthesis-vs-pool relationship.

H2: Ecosystem Metabolic Niche

The strong ecosystem separation in PC space is consistent with the metabolic niche
hypothesis: habitat physicochemistry (aquatic vs. terrestrial) selects for
fundamentally different community metabolic configurations, not just different taxa.
The dominance of carbon utilization pathways on PC1 (not amino acid pathways)
suggests that carbon substrate availability is the primary axis differentiating
ecosystem metabolic type at the community level, with amino acid biosynthesis
playing a secondary (PC2-level) role.

The finding that 17/18 aa pathways differ by ecosystem type is consistent with
Ramoneda et al. (2023), who found that auxotrophy rates for amino acids are
taxonomically structured and environmentally distributed; freshwater taxa may be
under different selective pressures regarding amino acid biosynthesis retention
than their soil counterparts.

Literature Context

• H1 direction aligns with Morris et al. (2012), who proposed that biosynthetic
  gene loss is selectively favoured when the environment reliably supplies the metabolite
  ("Black Queen"). Our data extend this to community-weighted scores at the landscape
  scale across diverse NMDC soil samples.

• H1 is consistent with Mallick et al. (2019) and Noecker et al. (2016), both
  of which demonstrated that community metabolomics is partially predictable from
  genomic/metagenomic functional profiles. Our approach using GapMind pathway
  completeness as the predictor is methodologically distinct (pathway-level completeness
  rather than flux-based or PICRUSt-style models) but arrives at comparable conclusions.

• H2 aligns with Danczak et al. (2020), who established a meta-metabolome ecology
  framework showing that community assembly processes govern which metabolites
  accumulate. Our PCA result extends this to genome-predicted metabolic potential:
  not just what metabolites are present, but what the resident taxa are capable of
  making, clusters strongly by ecosystem.

• Ramoneda et al. (2023) used GapMind to show that amino acid auxotrophy
  distributions are non-random across bacterial phyla and environments. Our
  community-weighted extension shows that these taxon-level patterns aggregate into
  measurable ecosystem-level signals in the NMDC dataset.

• Gowda et al. (2022) showed that metabolite dynamics in model communities are
  predictable from gene content. Our result for leucine (r = −0.33) and arginine
  (r = −0.30) at environmental scale, though weaker, is directionally consistent with
  this finding from controlled communities.

Novel Contribution

This study is, to our knowledge, the first application of GapMind community-weighted
pathway completeness scores to predict environmental metabolomics across a large
cross-habitat NMDC dataset. The integration of 305M GapMind pathway records (27,690
GTDB species × 80 pathways) with NMDC multi-omics for 220 samples spanning two
habitat types establishes a scalable approach for testing metabolic ecological hypotheses
from public data lakehouses.

The finding that community-level amino acid biosynthesis completeness negatively
predicts ambient amino acid concentrations — even in noisy, cross-study environmental
data — suggests that community metabolic potential is a meaningful predictor of
ecosystem chemistry beyond taxonomic composition alone.

Limitations

• Metabolomics technical heterogeneity: The 175 metabolomics samples derive from
  multiple NMDC studies. However, study-level analysis of the H1 subset reveals that
  95% of H1 samples (125/131) originate from a single NMDC study
  (nmdc:sty-11-r2h77870), so cross-study LC-MS protocol heterogeneity is largely
  absent as a confounder for leu and arg. The second study (n = 6) does not
  materially affect the results. For future analyses with broader multi-study coverage,
  study-level random effects could be modelled to improve power.

• Freshwater samples absent from H1: All 33 Freshwater samples lacked paired
  metabolomics in NMDC, so H1 is effectively a soil-only test. Whether BQH dynamics
  operate in freshwater communities at the same scale is untested.

• Abiotic feature coverage: Abiotic features (pH, temperature, total organic carbon,
  etc.) were absent (all NaN) for the 174-sample analysis matrix, preventing partial
  correlation tests controlling for environmental gradients. These variables could
  confound both H1 and H2 results.

• GapMind represents genomic potential, not expression: Pathway completeness
  indicates whether the genes are present, not whether they are actively expressed.
  Transcriptomic data would be needed to determine whether expressed biosynthesis
  correlates more strongly with metabolite pools.

• Compound identification: Metabolite-to-pathway matching relied on string-based
  compound name matching. A substring collision bug (isoleucine compounds matching the
  "leucine" pattern) was identified by automated review and corrected in NB04 cell-14
  using first-match-wins; results reflect the corrected run. Isoleucine is now included
  as a 13th testable pathway (r = −0.057, ns). KEGG compound IDs in the NMDC data are
  sparsely populated (2% annotation rate), limiting KEGG-based matching; cysteine,
  histidine, and lysine remain untestable due to absent compound detections.

• Chorismate metabolomics proxy: Shikimic acid and 3-dehydroshikimic acid are used
  as metabolomics proxies for chorismate pathway activity. These are upstream
  intermediates, not chorismate itself; their ambient concentrations reflect precursor
  availability rather than chorismate pool size. The chorismate BQH correlation
  (r = −0.038) should be interpreted cautiously.

• Genus-proxy-ambiguous taxon assignment: For ~1,352 Centrifuge taxa that match
  multiple GTDB clades (same genus, multiple species), one representative clade was
  selected by alphabetical tiebreaking on gtdb_species_clade_id. This is now
  documented in NB03 cell-17. The sensitivity of community completeness scores to this
  choice is bounded: genus_proxy_ambiguous taxa account for ~6.5% of mapped abundance.

• Sample size imbalance: Glutamine (n = 4) and proline (n = 9) could not be tested
  due to insufficient metabolomics coverage. These pathways would be informative
  given their biological importance.

1. Abiotic gradient controls and multi-study replication: The H1 signal could
   be strengthened (or confounders identified) by including available abiotic variables
   (pH, temperature, total organic carbon) in a partial correlation or mixed-effects
   model once NMDC abiotic data is populated for the metabolomics-overlap samples. The
   current H1 dataset is effectively single-study; replicating this analysis across
   additional NMDC studies with paired taxonomy + metabolomics would test whether the
   leu/arg signals generalise beyond the current cohort.

2. Extend to Freshwater samples: Obtaining or generating metabolomics data for the
   33 Freshwater NMDC samples would allow a direct test of whether BQH dynamics differ
   between aquatic and terrestrial communities — a key open question given their
   strikingly different pathway completeness profiles.

3. Transcript-based validation: Pairing metatranscriptomics with metabolomics for
   a subset of NMDC samples would allow comparison of expressed pathway completeness
   vs. genomic potential vs. metabolite abundance, directly testing whether expression
   (not just gene presence) drives the BQH signal.

4. KEGG compound ID expansion: Improving KEGG compound annotations in
   metabolomics_gold (3 currently missing aa pathways: cys, his, lys) would
   expand the H1 test set from 13 to 16 testable pathways, improving statistical power.

5. Cross-habitat partial correlation (soil subtype analysis): The 61 "Unknown"
   ecosystem samples likely represent a mix of soil subtypes or sediment environments;
   resolving their habitat identity (via NMDC ENVO annotations or study-level metadata)
   could improve H2 power and reveal finer-scale metabolic niche structure.

Sources

Generated Data

• Morris JJ, Lenski RE, Zinser ER (2012). "The Black Queen Hypothesis: Evolution of
  Dependencies through Adaptive Gene Loss." mBio. 3(2):e00036-12.
  DOI: 10.1128/mBio.00036-12

• Noecker C, Eng A, Srinivasan S, Theriot CM, Young VB, Jansson JK, et al. (2016).
  "Metabolic model-based integration of microbiome taxonomic and metabolomic profiles
  elucidates mechanistic links between ecological and metabolic variation." mSystems.
  1(1):e00013-15. DOI: 10.1128/mSystems.00013-15

• Mallick H, Franzosa EA, McIver LJ, Banerjee S, Sirota-Madi A, Kostic AD, et al.
  (2019). "Predictive metabolomic profiling of microbial communities using amplicon or
  metagenomic sequences." Nature Communications. 10(1):3136. PMID: 31316056

• Danczak RE, Chu RK, Fansler SJ, Goldman AE, Graham EB, Tfaily MM, et al. (2020).
  "Using metacommunity ecology to understand environmental metabolomes." Nature
  Communications. 11(1):6369. PMID: 33311510. DOI: 10.1038/s41467-020-19989-y

• Price MN, Deutschbauer AM, Arkin AP (2020). "GapMind: Automated Annotation of Amino
  Acid Biosynthesis." mSystems. 5(3):e00291-20. DOI: 10.1128/mSystems.00291-20

• Price MN, Deutschbauer AM, Arkin AP (2022). "Filling gaps in bacterial catabolic
  pathways with computation and high-throughput genetics." PLOS Genetics.
  18(7):e1010156. DOI: 10.1371/journal.pgen.1010156

• Ramoneda J, Jensen TBN, Price MN, Fierer N, Braendle C (2023). "Taxonomic and
  environmental distribution of bacterial amino acid auxotrophies." Nature
  Communications. 14(1):7608. DOI: 10.1038/s41467-023-43435-4

• Gowda K, Ping D, Mani M, Kuehn S (2022). "Genomic structure predicts metabolite
  dynamics in microbial communities." Cell. 185(3):530-546. PMID: 35085485.
  DOI: 10.1016/j.cell.2021.12.036

• Danczak RE, Goldman AE, Chu RK, Toyoda JG, Borton MA, et al. (2021). "Ecological
  theory applied to environmental metabolomes reveals compositional divergence despite
  conserved molecular properties." Science of the Total Environment. 788:147409.
  DOI: 10.1016/j.scitotenv.2021.147409

• Arkin AP, Cottingham RW, Henry CS, Harris NL, Stevens RL, Maslov S, et al. (2018).
  "KBase: The United States Department of Energy Systems Biology Knowledgebase."
  Nature Biotechnology. 36(7):566-569. DOI: 10.1038/nbt.4163