Study and batch confounding can reverse environmental generalization

Environmental associations estimated across heterogeneous studies can reflect study, batch, classifier, or unmeasured sample-processing structure rather than portable environmental effects. The Metadata Correlates of Eukaryotic Contamination in NMDC Prokaryote-Targeted Metagenomes provides a direct example: eukaryotic read fractions appeared environmentally predictable under random cross-validation but failed to generalize when entire studies were held out. Missing laboratory metadata further prevents direct separation of ecological effects from extraction, filtration, host-depletion, library-preparation, and related processing effects. ₁

Core claim

When environmental categories are strongly nested within studies, an environmental predictor can act as a proxy for study-specific protocols, populations, sequencing workflows, or unmeasured sample-processing factors. Random train/test splits allow closely related samples from the same study to occur in both partitions, whereas study-held-out validation tests whether an association transfers to a new batch or study. This makes study-batch-confounding-of-environmental-associations a direct extension of environment-embedding-geography and callability-limited-comparative-inference. ₁

The Ecotype Reanalysis — Environmental vs Human-Associated Species supports the distinction between an observed environmental signal and a transferable environmental effect, while refining it beyond study-held-out read-fraction prediction. In a genome-level comparison, environmental species had a median environment–gene-content partial correlation of 0.051, versus 0.084 for human-associated species; the one-sided Environmental > Human-associated test gave U=1536 and p=0.83. ₂ The continuous fraction-of-environmental-genomes analysis likewise found no relationship with correlation strength (rho=-0.085, p=0.25). ₂ Thus, the clinical composition bias was real but did not explain the weak environmental gene-content association, consistent with the possibility that environmental structure can be genuine yet non-portable or poorly aligned with the measured biological outcome. ₂

The original Ecotype Correlation Analysis supports this qualification at broader species scale: among 172 species with sufficient data, the median environment–gene-content partial correlation was 0.0025 versus 0.0143 for phylogeny, with phylogeny dominating in 60.5% of species; significant positive, negative, and null environmental effects occurred in 7.0%, 2.3%, and 90.7% of species, respectively. ₃ This refines the page’s interpretation: weak whole-genome environmental association is not by itself evidence of batch confounding, because environmental adaptation may act on specific gene subsets, while it also supports the hypothesis that AlphaEarth embeddings do not capture all ecologically relevant variation. ₃

Evidence from eukaryotic read fractions

The analysis quantified GOTTCHA2 relative eukaryotic abundance across 2,759 NMDC ReadbasedAnalysis runs from 9 studies. ₁ GOTTCHA2 detected eukaryotic reads in 77% of the runs, with a median eukaryotic fraction of 2.7%, a mean of 13.3%, and 20% of runs exceeding 20% eukaryotic reads. ₁ These fractions were classifier- and database-dependent: Kraken2 and Centrifuge produced approximately 0 domain-level Eukaryota signal in this deployment, whereas GOTTCHA2 was plastid- and eukaryote-aware. ₁ The GOTTCHA2 values should therefore be interpreted as relative or ordinal contamination measures rather than calibrated absolute fractions. ₁

The apparent environmental structure was substantial: aquatic freshwater samples had 99.5% eukaryotic detection and a plastid share of 1.00, terrestrial soil had 55.7% detection and a plastid share of 0.43, and plant-root samples had 100% detection and a plastid share of 0.03. ₁ Matrix differences were significant in a Kruskal–Wallis test with H=77.8 and p=1.3×10⁻¹⁷, and all pairwise matrix contrasts were significant after BH-FDR, where FDR means false discovery rate. ₁ These results establish strong observed environmental differences, but they do not identify which differences arise from sample ecology and which arise from unrecorded laboratory processing. ₁ Ecosystem type was not an independent confirmation of the matrix result: after excluding the Unknown/missing-metadata bucket, ecosystem_type reduced to the same three biomes (Soil / Freshwater / Roots) and the test was identical to the matrix test; an earlier version that retained the Unknown bucket, which had the highest median eukaryotic fraction, reported an inflated p≈10⁻⁵⁰ that was a missingness artifact and has been removed. ₁

Each biome was approximately 80–100% nested within a single study, preventing a clean separation of environment from study or batch in cross-study univariate comparisons. ₁ A gradient-boosted model of the logit-transformed eukaryotic fraction (the euk_logit response) using study_id alone achieved R²=0.24 under random cross-validation, while an environment model achieved R²=0.35 under the same validation scheme; all R² values quoted here for these models are for that euk_logit response rather than the untransformed fraction. ₁ Under GroupKFold, meaning cross-validation in which complete studies are held out as groups, the environment model achieved R²=−0.30 out of study, and the model combining environment with sequencing achieved R²=−0.39. ₁ The out-of-study environment model therefore performed worse than predicting the mean, while its out-of-study detection AUC was 0.56, approximately chance. ₁

This contradicts the interpretation that the stronger random-validation environment score demonstrates broadly transferable environmental prediction. ₁ The contrast supports the hypothesis that the pooled cross-study signal was substantially batch-driven and encoded study composition or processing history rather than a portable environmental effect. ₁ Missing protocol fields prevent direct adjustment for wet-lab variables that may differ systematically among studies. ₁

The report’s own numbers are not fully consistent on this point, and the discrepancy is left unresolved here rather than reconciled. Its prose states that environment explains no more variance than study_id alone, while its random cross-validation table reports an environment R²=0.35 against a study_id-only R²=0.24. ₁ The report’s discoveries text likewise equates the out-of-study R²=−0.30 with the study_id-only R², contradicting the table’s study_id-only random-CV R²=0.24; the two figures come from different cross-validation schemes and the inconsistency remains unresolved. ₁

The cross-study test is also underpowered. NMDC’s read-based taxonomy spans only approximately 9 studies, with one soil study accounting for approximately 43% of runs, and no on-system source provides per-sample eukaryotic read fraction across many studies; the within-study result is therefore demonstrated for one soil study only and may not extend to aquatic or host-associated collections. ₁

The report places this pattern in prior literature, citing Salter et al. (2014, BMC Biology, PMID 25387460) for biomass and batch as the master variable. ₁

Adding sequencing platform and depth did not improve study-held-out prediction. ₁ Measured sampling depth — not sequencing depth — was negatively associated with eukaryotic fraction in a whole-collection (cross-study) statistic, with Spearman ρ=−0.29, p=5.2×10⁻⁷, and n=292 runs with measured depth: shallower samples carried more eukaryotic DNA, consistent with surface plant or algal input. This statistic was not batch-controlled. ₁ Measured-depth runs came from a handful of non-soil studies, while the dominant NEON soil study recorded zero non-null depth values; the depth result is therefore suggestive only, not evidence of a controlled sampling-depth effect. ₁

Why within-study analysis changes, but does not settle, the conclusion

A batch-controlled analysis of 1,186 runs from the dominant NEON soil metagenome study used one sampling program with a constant protocol and batch. ₁ Within that study, local vegetation was associated with eukaryotic fraction at H=119.1 and p=7.6×10⁻²¹ across 11 env_local_scale levels — medians were highest in sedge/forb herbaceous soil at 23%, emergent wetland at 14%, dwarf scrub at 13%, and evergreen forest at 2%, and approximately 0 in deciduous forest, cropland, and pasture — and geography was associated across 47 sites at H=310.4 and p=2.4×10⁻⁴⁶. ₁ A model using local environment and geography achieved within-study five-fold R²=+0.17 ± 0.06, contrasting with cross-study out-of-study R²=−0.30. ₁

This supports a genuine fine-scale environmental signal when batch is held approximately constant and refines environment-embedding-geography: environmental metadata can be predictive, but apparent portability depends on validation and sampling structure. ₁ It does not establish that the result is free of protocol confounding or causal, because local environment and geography may still track sub-batches such as sampling campaigns. ₁ Nor does it establish portability to aquatic or host-associated collections. ₁

The ecotype reanalysis supports this portability qualification but does not reproduce the original absolute correlation scale: its all-species median partial correlation was 0.081, compared with 0.003 in the original analysis, a reported 27x difference. ₂ The reanalysis used all genomes with embeddings, including up to 3,505 genomes per species, rather than diversity-maximizing downsampling with a maximum of 250 genomes; the report states that the absolute values are not comparable, while the within-method environmental-versus-human comparison remains valid. ₂ This refines the page’s methodological warning: validation structure matters, but genome selection and extraction methodology can also change the apparent strength of an environmental association.

The original ecotype analysis supports the same caution but adds a different whole-genome result: it analyzed 13,381 genomes across 224 species and produced correlations for 172 species, finding that phylogeny generally dominated gene-content similarity and that environmental effects were often statistically non-significant. ₃ Its environmental embeddings covered only 28.4% of genomes, and geographic coordinates were often missing or imprecise; partial correlations also assume linear relationships between distance matrices. ₃ These limitations refine the interpretation of a null environmental effect: the result is consistent with weak or poorly measured ecological structure, but does not distinguish that explanation from genuine subset-specific adaptation or unmeasured confounding. ₃

Batch can dominate cluster structure in pooled multi-cohort data

Batch dominance is not confined to environmental metadata, and this supports the page’s core claim in a clinical modality. Pooling hmp2 and franzosa metabolomics on the 122 m/z-bridge metabolite panel, as the digest describes it, and clustering with PCA (principal component analysis) plus K-means at K=4 produced clusters that separated completely by cohort, not by diagnosis; PC1 explained 79 % of total variance and was essentially the cohort batch effect, and the cross-cohort LOSO (leave-one-study-out) ARI (adjusted Rand index, where 0 means chance-level agreement) was 0.000 against the 0.113 taxonomic-ecotype LOSO baseline. ₄ The project report gives the same outcome but describes the panel differently: a pooled 326-subject × 111-unique-feature m/z-bridge panel, with cluster structure cohort-driven rather than diagnosis-driven, PC1 explaining 79 % of variance, and cross-cohort LOSO ARI = 0.000 versus the 0.113 taxonomic baseline. ₅ The two panel descriptions — 122 metabolites in the digest, 111 unique features in the report’s method — are recorded as stated and are not reconciled here. _4,5

The cluster composition shows the nesting directly: Cluster 0 (86 subjects) was Franzosa 100 % with 46 Control + 30 UC + 10 CD, Cluster 1 (58) was HMP2 100 % with 29 CD + 23 UC + 6 nonIBD, Cluster 2 (134) was Franzosa 100 % with 78 CD + 46 UC + 10 Control, and Cluster 3 (48) was HMP2 100 % with 21 CD + 7 UC + 20 nonIBD, so every cluster was cohort-pure despite mixed diagnoses. ₅ PC1 separated HMP2 (PC1 ≈ +12) from Franzosa (PC1 ≈ −5) cleanly while PCs 2–15 collectively explained ~10 % of variance; within each cohort the clusters did carry diagnosis information (cluster × diagnosis χ²(6)=88.3, p=7e-17), but the pooled cluster boundaries were mostly co-aligned with cohort boundaries. ₅ The project’s figure NB09d_metabolite_ecotype_stability.png records the same structure in three panels — a PCA scatter colored by metabolite K=4 cluster with PC1 79 % variance labeled as the cohort batch, a per-cluster diagnosis-composition stacked bar at χ²=88, p=7e-17, and a stability summary giving cross-cohort LOSO ARI 0.000 against the taxonomic 0.113 baseline. ₅

This refines how stability metrics should be read: the within-pooled bootstrap ARI of 0.937 is misleading, because it measures how reproducibly pooled K-means recovers the cohort-batch structure under subsampling rather than whether the clustering is biologically informative. ₅ More generally, within-pooled bootstrap stability metrics measure batch reproducibility, not biological reproducibility, when batch dominates; the within-cohort metabolite-feature framework may still be useful for clinical translation inside a single laboratory’s analytical pipeline, and the documented failure is specifically the cross-cohort-portability question. ₄

The project’s own taxonomic-ecotype stability numbers make the same contrast quantitatively: bootstrap 80 %-subsample ARI was 0.13–0.17 on this data, while LOSO ARI across 8 independent sub-studies was 0.00–0.28 with mean 0.113, revealing per-substudy variation that bootstrap masks. ₅

Alignment varied sharply by sub-study rather than failing uniformly: LifeLinesDeep aligned well at ARI 0.21 with 85 % agreement and HansenLBS at 65 % agreement, while AsnicarF reached 38 % and VilaAV 17 %. ₅

The mechanism is modality-specific. Taxonomic relative-abundance values from metaphlan3 are unitless and compositionally constrained, a constraint that acts as a natural normalization across cohorts, whereas metabolite-feature clustering inherits absolute-intensity scale differences between LC-MS runs (instruments, ionization tuning, solvent batches) and requires explicit batch correction such as ComBat, SVA, RUV, or quantile normalization before clustering. ₅ Metabolite-feature ecotype clustering therefore failed as a primary clinical-translation basis without external batch correction, and taxonomic-feature ecotypes remained the project’s primary basis. ₅ These results feed cross-cohort-microbiome-portability and ecotype-clustering-validity.

When case and control come from different studies, the confound is unadjustable

Batch is not always something a model can absorb. curatedmetagenomicdata pools samples across ≈ 51 source studies, and in its ecotype-assigned slice the healthy and disease sub-studies do not overlap. ₅ The pitfalls digest records the same structure for that collection’s CMD_HEALTHY and CMD_IBD buckets, which bundle samples from ~50+ published studies: every healthy sample comes from a healthy-cohort study (LifeLinesDeep_2016, AsnicarF_2021, YachidaS_2019, …) and every CD (Crohn’s disease) sample from an IBD-cohort study (HallAB_2017, VilaAV_2018, LiJ_2014, IjazUZ_2017, NielsenHB_2014), so there is no sub-study that contains both CD and HC (healthy control) samples. ₆

The counts are explicit: in the integrated mart scoped to ecotype-assigned samples — 8,489 total, drawn from approximately 51 source studies — 45 sub-studies have ≥ 10 HC samples, 5 sub-studies have ≥ 10 CD samples, and 0 sub-studies have both. _5,6 A pooled CD-versus-HC contrast is therefore structurally unidentifiable rather than merely noisy, because the sub-study random effect perfectly predicts diagnosis, and it cannot be adjusted via a standard random-effects model on the pooled-cohort contrast; the project verified this empirically, with statsmodels.mixedlm silently failing to converge on every battery species. _5,6

The failure mode is silent. A pooled log_abundance ~ diagnosis + (1 | substudy) LME (linear mixed-effects model, which treats sub-study as a random intercept) is unidentifiable because the sub-study random effect perfectly predicts diagnosis, so the model either fails to converge or absorbs all of the CD effect into the random intercept; statsmodels.mixedlm with lbfgs returned without error but with no usable fixed-effect estimate on this design. ₆ The project’s own LME output table data/nb04c_lme.tsv is correspondingly empty, retained specifically to document that the pooled CD-versus-HC LME is structurally unidentifiable. ₅

The reporting consequence is that a pooled effect must be flagged rather than adjusted: study, cohort, sequencing center and sample preparation are all collinear with diagnosis in this design, so the effect should be quoted as-is with that explicit caveat and sanity-checked against the within-substudy meta-analysis, not “adjusted” by a sub-study random effect. ₆ Both documents cite Vujkovic-Cvijin et al. (2020), “Host variables confound gut microbiota studies of human disease,” Nature 587(7834):448–454, PMID 33149306, for this class of structural confounding. _6,5

This generalizes beyond the IBD corpus: the same unadjustable structure arises in any pooled public-dataset analysis — SRA study inventories, Qiita studies — where the case and control cohorts were collected by different groups, and the within-cohort CD-versus-nonIBD contrast is the design-consistent alternative. ₆

Pooled analysis can change the interpretation of a directional call

The same entanglement can change what a species-level call means. A v2 preliminary IBD report (2026-03-28) called clostridium-scindens CD-enriched at log₂FC = +2.67 from pooled Mann-Whitney differential abundance (DA, a test of whether a taxon’s abundance differs between groups). ₄ The digest describes C. scindens as a well-documented protective species — a secondary bile-acid producer via 7α-dehydroxylation, a TGR5 activator, present in ~79 % of healthy individuals, and an inhibitor of C. difficile — which is why the pooled enrichment was initially read as a false call. ₄

Three mechanisms were proposed for that presumed false call: a compositional or relative-abundance artifact, in which Ruminococcaceae loss in severe patients inflates the relative abundance of surviving C. scindens strains; strain heterogeneity, because not all C. scindens strains carry the bai operon, so species-level presence is not functional presence; and ecotype mixing, because pooled analysis averages across patient subgroups with different microbiomes. ₄ Compositional-aware DA (ANCOM-BC / MaAsLin2 / LinDA) run within patient ecotypes was the proposed fix, and the digest records it as proposed rather than as completed validation. ₄

A confound-free within-substudy analysis later contradicted the paradox-resolution reading rather than the enrichment itself: C. scindens is genuinely CD↑ under within-substudy confound-free analysis, with pooled CLR-Δ (centered-log-ratio difference) +1.18, FDR 1e-8, and 4/4 sign concordance, and the earlier within-ecotype non-significant call was a leakage self-selection artifact, since a leave-one-out refit showed C. scindens CD↑ in both E1 and E3 once the ecotype definition excluded the test species. ₄ Both the pooled preliminary test and the within-substudy analysis therefore report CD enrichment; what the revised analysis overturns is the reading of that enrichment as a compositional or ecotype-mixing false call. ₄

The design that produced this correction is a confound-free within-IBD-substudy CD-versus-nonIBD meta-analysis that resolves both the feature-leakage and substudy-confounding pitfalls of pooled-cohort DA; it yielded six actionable Tier-A targets — H. hathewayi (4.0), M. gnavus (3.8), E. coli (3.6), E. lenta (3.3), F. plautii (3.3), E. bolteae (2.8). ₅ This supports selection-on-outcome-leakage.

The identifiable contrast is CD-versus-nonIBD within IBD sub-studies: four curatedMetagenomicData studies have ≥ 10 CD and ≥ 10 nonIBD samples — HallAB_2017 89/73, LiJ_2014 76/10, IjazUZ_2017 56/38, NielsenHB_2014 21/248, aggregating to 242 CD / 369 nonIBD — and within each sub-study CLR-Δ is computed on the same samples processed the same way by the same group, then combined across sub-studies by inverse-variance weighted meta-analysis. _5,6

The design makes two guarantees explicit. The clustering axis and the DA axis are disjoint — samples are partitioned by ecotype and CLR-Δ is computed across sub-studies within a partition — which eliminates feature leakage. ₅ And the within-substudy contrast has no study-level variation, which removes between-study variation from the differential-abundance comparison. ₅

The repair pipeline that produced the replacement analysis ran NB04b (bootstrap confidence intervals, leakage-bound sensitivity, leave-one-out refit, Jaccard permutation null, ecotype stability) → NB04c (confound-free within-IBD-substudy CD-versus-nonIBD meta plus linda implemented in pure Python) → NB04d (rigor-controlled stopping rule) → NB04e (within-ecotype × within-substudy CD-versus-nonIBD meta), and it recorded curatedMetagenomicData substudy-nesting unidentifiability, feature leakage in cluster-stratified DA, and the LME-unidentifiability documentation as its generalizable lessons. ₅

Across the four IBD sub-studies (242 CD / 369 nonIBD pooled), the curated battery — which the report labels a 14-species battery although its displayed table does not list that many species rows — produced the expected canonical CD signature, pathobionts up and protective commensals down, with strong sign concordance: ₅

• Higher in CD: Mediterraneibacter gnavus +5.13 (FDR ~0, 4/4), Eggerthella lenta +2.30 (FDR 4e-9, 4/4), Escherichia coli +1.43 (FDR 2e-4, 3/4), Clostridium scindens +1.18 (FDR 1e-8, 4/4), Enterocloster bolteae +1.09 (FDR 3e-6, 4/4), Hungatella hathewayi +0.92 (FDR 5e-4, 3/4), Bilophila wadsworthia +0.07 (not significant, 3/4). ₅
• Lower in CD: Lachnospira eligens −1.01 (FDR 4e-3, 3/4), Roseburia intestinalis −1.14 (FDR 2e-3, 4/4), Akkermansia muciniphila −1.30 (FDR 3e-3, 3/4), Faecalibacterium prausnitzii −1.67 (FDR 9e-8, 4/4), Roseburia hominis −1.77 (FDR 4e-7, 4/4), Coprococcus eutactus −3.09 (FDR 4e-15, 4/4). ₅

This contradicts the earlier within-ecotype calls for several species. The within-ecotype analysis had called F. prausnitzii, R. hominis and L. eligens CD↑ within both E1 and E3 — the “Simpson’s paradox” of its original section 5 — whereas the confound-free analysis shows them CD↓, consistent with their classical protective-commensal role, and the earlier C. scindens “n.s.” call within both ecotypes is likewise contradicted. ₅ The report attributes those reversals to a compound artifact of feature leakage plus the pooled-cohort substudy × diagnosis confound, and notes that both compositional-bias-aware DA methods tried on the within-ecotype subsets, CLR-MW and LinDA, share the bias, so evidence counts from within-ecotype methods alone do not resolve it. ₅

The disagreement is sign-level for the engraftment pathobionts: confound-free within-substudy CD effects were M. gnavus +5.13 (FDR 0, 4/4) against retracted within-ecotype values of E1 CD↓ −2.7 and E3 CD↑ +1.6; E. lenta +2.30 (FDR 4e-9, 4/4) against E1 CD↓ −2.4 and E3 CD↓ −0.8; E. coli +1.43 (FDR 2e-4, 3/4) against E1 CD↓ −1.8 and E3 CD↓ −0.9; E. bolteae +1.09 (FDR 3e-6, 4/4) against E1 CD↓ −1.6 and E3 n.s. +0.6; H. hathewayi +0.92 (FDR 5e-4, 3/4) against E1 CD↓ −1.4 and E3 n.s. −0.1; K. oxytoca fell below the prevalence filter. ₅

The stated mechanism is exactly the nesting this page describes: the within-ecotype DA was systematically reversing their sign because the HC samples in each ecotype came from entirely different source studies than the CD samples in the same ecotype. ₅

The pitfalls digest records the downstream reclassification: an independent-design cross-check using the confound-free within-IBD-substudy CD-versus-nonIBD contrast found that most of the earlier Tier-A E1 candidates had negative within-study effects, making them ecotype markers rather than CD drivers. ₆

The report abstracts these into methodology lessons beyond its own findings: curatedMetagenomicData substudy-nesting unidentifiability, feature leakage in cluster-stratified DA, LOSO ARI as an honest stability metric, pool ≠ flux for pathway versus metabolite DA, multi-line cross-corroboration as a portable rigor pattern, and cohort batch dominating clustering for absolute-intensity multi-cohort modalities. ₅

A classifier result in the same project is a further instance: is_ibd dominated classifier AUC and the performance collapsed on a single-cohort test. ₅ The minimal classifier predicted E1 for 19/22 UC Davis patients; in the training cohort IBD samples split ~58 % E1 / ~40 % E3 / ~2 % E0 / ~0 % E2, so its dominant learned rule was “is_ibd = 1 → E1,” and applied to the all-CD UC Davis cohort where is_ibd is constant that rule collapses to the marginal mode. ₅ The extended classifier’s training subset was 702 E1 / 959 E3 / 3 E0 / 11 E2 — effectively an E1-versus-E3 binary problem — and severity markers did not separate the two reliably. ₅

The report abstracts a transferable diagnostic from this: the gap between OvR-AUC (one-vs-rest area under the curve, computed for each class against all others) and per-patient agreement flags classifier-utility overstatement when a cohort-axis variable dominates the features. ₅

Site aggregation can pull a site-local association toward null

Aggregation across sites can move a reported signal in the opposite direction from aggregation across studies. The serology analysis crossed six EU-deduplicated serology axes (ANCA, ASCA, CBir1, IgA-ASCA, IgG-ASCA, OmpC) with 8 species — 6 actionable Tier-A core plus 2 module-anchor commensals — for 48 tests, using partial Pearson r residualized on site dummy variables with BH-FDR. ₅

The cohort-aggregate signal was dominated by ANCA × M. gnavus (+0.31), but the site-stratified breakdown showed most associations concentrated in 1 of 2 productive sites (Harvard or CCHMC) and diluted when cohort-aggregated; with Harvard 33 + CCHMC 31 + Emory 3 subjects, single-site partial r reached |r|≈0.40–0.46 individually, while the top observed aggregate was |r|=0.31 at FDR=0.40, so the strict |r|>0.40 plus FDR<0.10 plan threshold was not met and strict H3e was not supported. ₅ Site-stratified single-site r reached +0.46 at Harvard while cohort aggregation pulled toward null, and the report records a structural single-cohort caveat. ₅ This refines the page’s claim: pooling does not always inflate an association — here cohort aggregation pulled site-local correlations toward null, and those site-local correlations themselves remain unvalidated because the planned significance threshold was not met. ₅

Bile-acid signals were weaker cross-cohort than within-cohort in the same project: primary bile acids replicated at 57 % (4/7) and secondary bile acids at 40 % (2/5). ₅ The report treats the paired sample-level substrate–product signature for F. plautii, E. lenta and E. bolteae measured within HMP2-paired samples as the more robust evidence stream than cohort-aggregate cross-cohort DA, because bile-acid pool sizes are more susceptible to dietary and antibiotic confounding than within-cohort species × metabolite correlation; it states explicitly that this paired finding does not require Franzosa replication, and that the cohort-level result documents only that bile-acid-pool DA is more cohort-dependent than within-cohort correlation. ₅

Design confounds inside a single study

Holding study constant does not remove confounding when the design itself is unbalanced. The Harvard Forest Long-Term Warming — DNA vs RNA Functional Response design has 42 biosamples in a factorial layout of treatment (control vs heated) × horizon (organic vs mineral) × incubation (direct vs lab-incubated), with a DNA cohort of n=28 (mineral direct + organic incubated) and an RNA cohort of n=39; coverage is unbalanced because the underlying pipeline did not produce metagenomes from organic-direct or some mineral-direct samples, which drives the sensitivity caveat on its first hypothesis. ₇

The hypothesis that RNA composition shifts more than DNA under warming was not supported once a horizon × incubation confound was removed. In the paired n=25 subset (samples with both DNA and RNA), PERMANOVA (permutational multivariate analysis of variance) pseudo-F for treatment was 3.12 in DNA (R²=11.9%, p=0.020) versus 0.76 in RNA (R²=3.2%, p=0.60), seemingly contradicting the hypothesis in the opposite direction; but the paired subset confounds horizon with incubation, because every organic sample is incubated and every mineral sample is direct. ₇

A sensitivity analysis using direct samples only, with no incubation confound, gave treatment R²=12.7% for DNA × mineral (n=14, F=1.75, p=0.081), 11.4% for RNA × mineral (n=11, F=1.15, p=0.254), and 10.0% for RNA × organic (n=14, F=1.33, p=0.190). ₇

DNA and RNA pools therefore showed comparable treatment R² of 10-13%, and the premise that the transcript pool is more sensitive was not supported in this dataset — a null result that depends entirely on removing the design confound. ₇ This feeds multi-omics-integration.

The residual confound is acknowledged rather than removed: organic-horizon DNA samples are all lab-incubated, since direct organic DNA was not produced by the pipeline, making it impossible to factor incubation cleanly out of the DNA pool’s organic-horizon analysis. ₇

Cross-site reproducibility of the project’s methanotrophy-up signature remains untested here; the report names SPRUCE peatland study nmdc:sty-11-33fbta56 and Alaskan permafrost-thaw study nmdc:sty-11-db67n062 as comparison datasets for that test. ₇

Time can be confounded with material in the same way. In the SSO Subsurface Community Ecology — Spatial Structure, Functional Gradients, and Hydrogeological Drivers, sediment cores were collected once per well (Feb-Mar 2023) with no within-well temporal replication for sediment, so the 18-month offset between sediment (2023) and groundwater (2024) confounds material type with time, and cross-material temporal comparisons must be read with that caveat. ₈ This refines subsurface-hydrogeological-zonation.

The same report states the offset concretely: sediment cores (Feb-Mar 2023) and groundwater (Sep 2024) were collected 18 months apart, so seasonal or plume dynamics could confound the comparison. ₈

Partial control is also a limitation in its own right: in the Contamination Gradient vs Functional Potential in ENIGMA Communities, site structure is represented by coarse location_prefix effects rather than full hierarchical or random-effects modeling, so residual site structure remains a candidate explanation for reported associations. ₉

Conditioning on batch changes what an explained-variance figure means

A batch term in the model changes what an explained-variance figure refers to. In the soil metal analysis, db-RDA — distance-based redundancy analysis, a constrained ordination of community distance matrices on measured predictors — reported R² = 0.799 with p = 0.005 from a permutation test with 999 permutations, described as metal concentrations explaining 80% of variance in community COG (clusters of orthologous groups) profiles after conditioning on batch/project effects. ₁₀

That figure refines rather than contradicts this page’s claim, provided its scope is stated: the db-RDA conditions on project accession before fitting the metal concentration predictors, and the report’s own caveat is conditional — if batch effects (project) account for the majority of community COG variance, then the reported R² describes how well metals explain the residual variance, not the total variance. Which case holds is not determined, because the unconditional R² of metals alone, without first removing project effects, is not reported and may be substantially lower; the report states that it must be reported alongside the conditional figure. ₁₀

Batch-aware validation can also leave nothing predictable. In the global soil functional-genomics analysis, all three predictive model families had negative out-of-sample R² — Soil & Climate R² = −0.205 ± 0.197, Geochemical R² = −0.331 ± 0.071, Industrial R² = −0.221 ± 0.042 — meaning they predict worse than the training mean. ₁₁

The report declines to read that as a clean biological null, and this refines the page’s own use of negative out-of-study R² as evidence of non-portability: three causes have not been distinguished — train/test distributional shift due to spatial autocorrelation within GroupKFold folds, high-leverage outlier samples in test folds that dominate the mean squared error, and genuine unpredictability of functional gene counts from the measured predictors at global scale. Only the third supports a strong null-result interpretation; the first two would indicate a modelling failure rather than a biological null. ₁₁

Single-study dominance and pseudo-replication inflate significance

A nominally cross-study dataset can be one study. Study-level analysis of the 131-sample H1 dataset showed that 95% of samples (125/131) come from one NMDC study, nmdc:sty-11-r2h77870, with the remaining 6 samples from a second study. ₄ The project report names that second study as nmdc:sty-11-547rwq94 and reports that excluding its 6 samples leaves the arginine association at r = −0.264, p = 0.019, n = 78, so the arginine signal is not driven by the minor study; on that basis it treats cross-study LC-MS protocol heterogeneity as effectively absent as a confounder for leucine, whose samples all come from the dominant study. ₁₂

Within that composition, the leucine signal is reported robust to soil-only stratification (q = 0.022 unchanged) and to study-level blocking, with the 95%-single-study structure stated alongside the claim rather than resolved. ₄ The soil-only re-run, performed on the 125 Soil-ecosystem samples after excluding the 6 Unknown-ecosystem samples in the H1 dataset, reports leucine remaining FDR-significant at r = −0.390, q = 0.022, n = 62 — identical to the full-cohort result because the leucine samples are all Soil — arginine losing FDR significance at r = −0.264, q = 0.117, n = 78 while remaining directionally consistent, and the binomial sign test unchanged at 11/13 negative, p = 0.011. ₁₂

The same report keeps the heterogeneity caveat pointing in both directions: its 175 metabolomics samples derive from multiple NMDC studies, yet because 125/131 H1 samples originate from a single study, cross-study LC-MS protocol heterogeneity is largely absent as a confounder for leucine and arginine and the second study (n = 6) does not materially affect the results; it proposes modelling study-level random effects to improve power in future analyses with broader multi-study coverage. ₁₂

Single-study dominance can also masquerade as geography. The global metal-resistance biogeography report lists expedition-level clustering as an untested check, proposing to test whether the Atacama (OR=9.8) and USA clusters are single-study artefacts by checking distinct sample_accession prefixes per hotspot. ₁₃

Non-independence also inflates significance directly. In the ENIGMA Carbon Census — Tiered Knowledge Census of 83 Enrichment Compounds, each metagenome was treated as an independent observation in a rank test over compositional, zero-inflated relative abundances; with thousands of non-independent samples this inflates significance massively, and all 83 genera reach q<0.05, many at q~1e-70. ₁₄ This supports pooled-run-pseudoreplication-and-metadata-label-noise.

The report accordingly treats its environmental enrichment ranking as exploratory and proposes replacing the per-sample rank test with a study-aware mixed model, or a sample-level permutation respecting study structure, to obtain defensible enrichment statistics. ₁₄ This refines confirmatory-exploratory-ecological-association-discordance.

Metadata, classifier, and join limitations

NMDC did not populate DNA-extraction kit, size fractionation or filtration, host-depletion method, or library-preparation fields, and has_filtration was all false in the processing metadata. ₁ These unmeasured factors can remain residual explanations for apparent environmental differences even after measured sequencing variables are included. ₁ Classifier and reference-database compatibility must therefore be audited before comparing contamination or taxonomic fractions across studies; the GOTTCHA2-specific measurement scale supports this qualification. ₁

The analysis used workflow_run_id rather than the biosample level because 1,067 of 2,759 runs were pooled from multiple biosamples, avoiding pseudo-replication from biosample-level joins. ₁ Metadata linkage connected workflow runs to biosamples and studies through NMDC tables and linked 99%+ of classified runs. ₁ For pooled runs, environmental and collection metadata came from a representative biosample selected by MIN(biosample_id); heterogeneous pooled biosamples could therefore introduce predictor label noise. ₁ The report characterizes this as a conservative bias that can weaken associations rather than manufacture them. ₁ This refines pooled-run-pseudoreplication-and-metadata-label-noise and environment-embedding-geography. ₁

Missing batch metadata is not confined to that collection: the lignin community-enrichment project reports no technical metadata on extraction or library-preparation batches, so batch effects there cannot be formally assessed at all. ₁₅

Incomplete statistical reporting limits what a reader can audit in the same way: a comparison of community-weighted mean values between EVO09 and EVO17 is reported as a Mann–Whitney test over 52 versus 81 samples with p = 0.0012, while the U statistic itself is not reported. ₁₆

The ecotype reanalysis further supports controlling unequal observation counts: environmental species had a 21% NaN partial-correlation rate, compared with 7% for human-associated species, so environmental species were disproportionately removed by missing-correlation filtering. ₂ Because this filtering would, if anything, favor finding a stronger environmental signal, its failure to do so does not support clinical sampling bias as the explanation for the null comparison.

The original ecotype analysis refines this metadata concern from another angle: AlphaEarth embeddings were available for only 28.4% of genomes, and host-associated geographic coordinates may represent collection sites rather than the organisms’ actual microenvironments. ₃ It therefore suggests the hypothesis that incomplete or ecologically misaligned environmental metadata can attenuate environmental associations even when study or batch is not the primary explanation. ₃

Implications for environmental inference

Random cross-validation should not be treated as evidence of environmental generalization when samples from the same study or batch can occur in both training and test partitions. ₁ Study-held-out GroupKFold or within-study contrasts are more appropriate tests of whether an association survives batch changes. ₁ Absent extraction, filtration, size-fractionation, host-depletion, and library-preparation metadata should be treated as unresolved confounders, not as evidence that protocol effects are absent. ₁

The ecotype result refines this implication by showing that correcting the environmental-versus-human comparison for genome-level classification and using a consistent within-method analysis still produced no stronger environmental association: environmental species had median partial correlation 0.051, while human-associated species had 0.084, with p=0.83 for Environmental > Human-associated. ₂ Embedding similarity therefore need not represent ecologically relevant variation; the report notes that environmental embeddings had a 3.4x geographic-signal ratio versus 2.0x for human-associated embeddings, while the captured environmental variation did not strongly predict bacterial gene content. ₂

The original ecotype analysis supports prioritizing functional subsets over whole-genome association: it proposes testing COG functional categories, including V-Defense and L-Mobile, and comparing gene content within identified ecotype clusters. ₃ This refines the implication that a weak genome-wide signal should not be interpreted as absence of ecological adaptation, particularly when the environmental representation is incomplete or the relevant effects are locus-specific. ₃

Tensions

The cross-study analysis shows strong environmental differences, including environment-model R²=0.35 under random cross-validation, while study-held-out validation gives R²=−0.30 and detection AUC=0.56. ₁ The within-NEON analysis nevertheless gives positive predictive performance, with R²=+0.17 ± 0.06. ₁ These findings are not mutually exclusive: fine-scale environmental structure may be real while the pooled cross-study association is not portable because study, batch, and unmeasured protocol variables are entangled. ₁ Constant protocol and batch reduce confounding but do not remove possible sub-batch confounding from sampling campaigns. ₁

The ecotype reanalysis refines this tension rather than resolving it: environmental species did not show stronger environment–gene-content correlations than human-associated species, despite a confirmed clinical sampling bias. ₂ Its absolute median correlation differed from the original analysis because genome inclusion and downsampling differed, so the reanalysis cannot be used to compare absolute effect sizes across methods; its within-method group comparison remains valid. ₂

The original ecotype analysis introduces a related scale and interpretation tension: it reports a median environment partial correlation of 0.0025 versus 0.0143 for phylogeny and no significant environmental effect in 90.7% of species, whereas the ecotype reanalysis reports an all-species median partial correlation of 0.081 versus 0.003 in the original analysis. _3,2 These values must not be reconciled by averaging: the reanalysis attributes the discrepancy to different genome inclusion and downsampling, and the original analysis also had only 28.4% embedding coverage. _3,2 The evidence therefore supports a robust caution about portability and measurement coverage, but does not establish whether the remaining weak signal is primarily ecological, methodological, or confounded.

The discoveries digest states that cross-study LC-MS protocol heterogeneity is effectively absent as a confounder, contrary to the naive expectation for NMDC metabolomics, and that the leucine BQH signal is internally consistent within one analytical protocol — while the same cohort includes 6 samples from a second study. ₄ These statements are retained side by side rather than merged: what the digest reports is internal consistency within one analytical protocol, and 125 of the 131 samples come from a single study. ₄

The IBD project is also internally inconsistent about which serology correlation belongs to which aggregation level: a future-directions line describes a multi-cohort meta-analysis as resolving a “cohort-aggregate vs site-stratified r=0.46 vs +0.31 gap,” while its results section describes +0.46 as the Harvard site-specific value and +0.31 as the site-adjusted cohort value. ₅ The labeling contradiction is recorded here rather than resolved, and the proposed multi-cohort analysis is the report’s own stated route to settling it. ₅

Open Directions

• Recover DNA-extraction kit, filtration, size-fractionation, host-depletion, and library-preparation records from laboratory information systems, then fit hierarchical or batch-stratified models to test which protocol variables explain residual eukaryotic-fraction variation. ₁
• Reanalyze the 2,759 runs with classifier-compatible reference databases and a common taxonomic estimator, testing whether environmental associations persist after database effects are controlled. ₁
• Use additional NMDC studies with overlapping biome categories and study-held-out GroupKFold to test whether prediction remains above the original cross-study R²=−0.30. ₁
• Reanalyze the 1,186-run NEON subset with sampling-campaign identifiers, restricted permutations, or mixed-effects models, testing whether the within-study R²=+0.17 ± 0.06 survives sub-batch control. ₁
• Separate pooled and unpooled workflow runs, retain biosample-level protocol and environmental labels, and compare run-level models with measurement-error or partial-pooling models to quantify how the 1,067 pooled runs affect environmental estimates. ₁
• Expand the analysis to studies with recorded sequencing depth and laboratory metadata, then use GroupKFold by study to ask whether protocol-aware models generalize beyond the dominant NEON soil study. ₁
• Perform a paired downsampled-versus-full-genome extraction using the same species and environment labels to determine why the original and reanalysis median partial correlations were 0.003 and 0.081, respectively, without treating the 27x difference as a biological effect. ₂
• Add genome count as a covariate and repeat the environmental-versus-human comparison with structured ENVO terms and functional gene subsets, testing whether sampling depth, ontology resolution, or whole-genome Jaccard distances mask environment-specific associations. ₂
• Reanalyze the 172-species ecotype dataset with alternative embedding distances and direct environmental metadata, then test COG categories such as V-Defense and L-Mobile and compare gene content within ecotype clusters to determine whether subset-level effects are hidden by whole-genome distances. ₃
• Import a many-study raw-read resource — the report proposes the SPF corpus, described as 136K metagenomes across thousands of studies — or compute eukaryotic fraction across a broader collection, breaking the batch confound and enabling a properly powered cross-study correlate analysis. ₁
• Apply explicit batch correction (ComBat, SVA, RUV, or quantile normalization) to the m/z-bridge metabolite panel before clustering, then recompute cross-cohort LOSO ARI against the 0.113 taxonomic baseline to test whether any diagnosis-driven structure survives removal of the cohort batch effect. _4,5
• Replace the per-sample rank test behind the soil–freshwater genus enrichment with a study-aware mixed model or a sample-level permutation respecting study structure, closing the pseudo-replication gap that currently drives all 83 genera below q<0.05. ₁₄
• Repeat the warming comparison with a balanced direct-sample design, or with the named NMDC warming studies as replicates, to test whether the comparable DNA and RNA treatment R² of 10-13% holds once organic-direct metagenomes exist. ₇
• Add within-well temporal replication for sediment, or resample sediment and groundwater in the same season, to separate material type from the 18-month collection offset. ₈
• Refit the contamination functional-potential models with full hierarchical or random-effects site terms in place of coarse location_prefix effects, testing whether site structure absorbs the reported associations. ₉
• Re-run the serology × species screen with site as a stratum rather than only as a residualized covariate, using the Harvard and CCHMC subsets separately, to test whether the single-site |r|≈0.40–0.46 associations survive multiple testing within site. ₅
• Re-run the global soil functional-gene models with spatial blocking to decompose the negative out-of-sample R² into distributional shift, outlier leverage, and true unpredictability, since only the last of the three licenses a strong null-result reading. ₁₁
• Report the unconditional metals-only db-RDA R² alongside the project-conditioned R² = 0.799 so that the share of total, rather than residual, community COG variance attributable to metals is known. ₁₀
• Add available abiotic variables (pH, temperature, total organic carbon) to the H1 metabolite–taxonomy analysis in a partial-correlation or mixed-effects model once NMDC abiotic data is populated for the metabolomics-overlap samples, and replicate across additional NMDC studies with paired taxonomy and metabolomics, to test whether the leucine and arginine signals generalise beyond an effectively single-study cohort. ₁₂
• Check distinct sample_accession prefixes per metal-resistance hotspot to test whether the Atacama (OR=9.8) and USA clusters are expedition-level single-study artefacts. ₁₃
• Apply the design-consistent within-cohort case-versus-control contrast to other pooled public inventories (SRA study inventories, Qiita studies) whose case and control cohorts were collected by different groups, instead of attempting random-effects adjustment of the pooled contrast. ₆