Forecasts of whether marine species can evolve fast enough to survive climate change rely largely on single-stage tolerance assays, tacitly assuming that tolerance measured at one developmental stage reflects adaptive capacity across the organism’s life cycle. We test this assumption in a common marine bivalve (Mercenaria mercenaria) and find that it fails in sign, not degree: larval and post-metamorphic tolerance to end-century ocean acidification are strongly negatively correlated (r = −0.53, 95% CI [−0.69, −0.33]) across family-level, pedigree-conditioned and animal model estimates, with pervasive rank reversals among the top performers of each stage. A multivariate analysis of the among-family covariance matrix G shows why: under acidification, 57% of the among-family variance collapses onto a single dominant axis that loads positively on every larval trait and negatively on every post-metamorphic trait, whereas under ambient conditions the same families show no cross-stage covariance at all. Acidification therefore does not merely depress performance—it reorganizes the genetic architecture into a cross-stage conflict. Under ramping RCP 8.5 forcing this conflict geometry dictates the direction—not merely the magnitude—of adaptation: only balanced multi-stage selection preserves family diversity, whereas single-stage selection systematically erodes the opposite axis. In complex life cycles, resilience is a geometric property of ontogeny that forecasts cannot compress into a single-stage without inverting the sign of the answer.

Fig 1. Three competing models show that climate tolerance can integrate, decouple or reverse across metamorphosis.

This figure is a conceptual schematic, not data: the points and lines are generated from a simulation in which 51 hypothetical families are drawn from bivariate normal distributions with imposed cross-stage correlations of +0.53, 0.00 and −0.53; the trajectories in d–f are schematic curves illustrating the qualitative expectation under each architecture. It illustrates the three hypotheses the experiment is designed to distinguish and the trajectory signatures each would produce; the corresponding empirical results are shown in Figs 2 and 6. a–c, Three alternative cross-stage architectures for family-level larval and post-metamorphic OA tolerance: positive covariance (co-adaptation), zero covariance (modularity) and negative covariance (conflict). Lines show the simulated relationships; shaded bands indicate 95% bootstrap intervals on the simulated draws. d–f, Expected trajectory signatures under uniform larval-acting selection. Positive covariance improves both stages, modularity improves the selected stage without changing the other, and conflict improves larval performance while eroding post-metamorphic performance, producing a lower net life-cycle composite. The numerical data underlying this figure are provided in https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.g001

Metamorphosis reverses family rankings for climate tolerance

We first tested the strongest prediction of the conflict model: the families that look climate-resilient before metamorphosis should not be the same families that look resilient afterward. For each family, we calculated stage-specific composite tolerance as weighted geometric means of OA/ambient trait-retention ratios, so that tolerance reflects stage-wide performance rather than a single endpoint (Materials and methods, Eqs. 2 and 3; S3 Table). The cross-stage answer was strongly negative: larval and post-metamorphic tolerance were inversely associated across families (Pearson r = −0.530, 95% CI [−0.691, −0.332], P = 6.2 × 10⁻⁵; Spearman ρ = −0.548, P = 3.1 × 10⁻⁵; Fig 2a). The reversal was visible at the level that matters for selection and screening. Of the 13 families in the top larval quartile, 10 fell below the adult median; of the 13 families in the top adult quartile, 11 fell below the larval median; no family occupied the top quartile at both stages (Fig 2b, 2c). Thus the first empirical result is not simply that OA affects two stages differently, but that metamorphosis changes which genetic families look resilient.

The reversal is robust to pedigree structure and tolerance definitions

We next asked whether the inversion was a structural property of the pedigree rather than a consequence of one scoring rule or a few influential families. A 5,000-shuffle permutation null placed the observed cross-stage correlation in the extreme negative tail (P < 1/5,000; Fig 2d; S4 Table). Removing sire and dam structure left a negative residual association (r = −0.486, P = 3.0 × 10⁻⁴; Fig 2e), and leave-one-sire/dam analyses never changed the sign (S3d Fig). The result also survived five alternative tolerance definitions, including equal-weight geometric means, arithmetic means of ratio terms, and adult composites omitting gonadal index or MO₂ (S3 Fig; S5 Table). Quantitative-genetic summaries pointed in the same direction: collapsing to the 12 shared sires gave r = −0.66 (95% CI [−0.90, −0.14], P = 0.019), dam-level summaries remained negative, and a bivariate pedigree animal model supported a negative cross-stage genetic correlation (S4 Fig; S6 Table). The sign is therefore not a family-outlier effect, not a composite-weight artifact and not created by ignoring the partial pedigree.

Ocean acidification imposes different filters before and after metamorphosis

Having established the inversion, we then asked how acidification imposes within-stage filters that could generate it. In larvae, OA reduced vial-level survival from 0.573 ± 0.004 s.e. to 0.432 ± 0.007 s.e. (β = −0.587 ± 0.028 s.e.; Wald χ² = 439.56, P < 0.001; S7 Table; Fig 3a, 3b). Shell area, shell thickness and shell hardness declined, whereas days to metamorphosis increased by 1.89 d (all P < 0.001; S8–S11 Tables; Fig 3c–3f). Importantly, acidification did not erase selectable family variance: family-level intraclass correlations (ICCs) were 0.187 for survival, 0.337 for shell area, 0.314 for thickness, 0.289 for hardness and 0.258 for metamorphosis (Fig 3g, 3h). Using the Lande–Arnold/Arnold–Wade selection-analysis framework [17,36], we then found that OA steepened directional selection toward larger first-shell area and faster metamorphosis (Fig 4a,4b): the standardized selection differential for shell area rose from +0.17 under ambient conditions to +0.28 under OA, and that for days to metamorphosis from −0.10 to −0.21 (treatment × trait interaction P < 10⁻⁵ for both; S12 Table), narrowing larval fitness around a high-calcification, rapid-development combination (Fig 4c,4d), with the differentials summarized in Fig 4e. This larval filter therefore does not merely depress performance; it redefines larval success under OA as rapid passage through a calcification bottleneck.

After metamorphosis, acidification acted on a different physiological accounting system. The juveniles posted a near-unity survival rate over the 90-day nursery, so at this stage acidification is expressed not as mortality but as a maintenance and reproduction cost, and gonadal index measured after two years of common-garden grow-out is the closest available proxy for post-metamorphic fitness. Shell growth declined from 1.395 ± 0.013 to 0.937 ± 0.023 mm month⁻¹, dissolution scores rose from 1.07 ± 0.05 to 2.01 ± 0.09, MO₂ increased by 45%, and gonadal index was still depressed by 26% 2 years later (S13–S16 Tables; Fig 5a–5d). Family-level variance again remained repeatable (ICC_family = 0.39 for growth, 0.30 for dissolution, 0.35 for MO₂ and 0.34 for gonadal index; Fig 5e, 5f), but the axis was no longer the larval one. Within the post-metamorphic stage, the two fitness-proximate traits moved together rather than in opposition: families that retained shell growth under acidification also retained gonadal output (r = +0.75, P = 2 × 10⁻¹⁰; Fig 5g). Growth nevertheless leaves 44% of the variance in reproductive retention unexplained, so it is an incomplete stand-in for post-metamorphic tolerance, and families with the highest metabolic cost under acidification tended to have the lowest delayed gonadal index, although this trend was not statistically significant (r = −0.27, P = 0.054; Fig 5h). The important contrast is therefore not within the post-metamorphic stage but between stages: family reaction norms make the switch visible, because the families that stand out as larval winners under acidification are not the families that stand out as post-metamorphic winners (S5 Fig). Metamorphosis reweights performance from rapid first-shell construction and developmental tempo to chronic maintenance, dissolution resistance, metabolism, and reproductive allocation.

Acidification reorganizes the multivariate architecture into a cross-stage conflict

Composite scores show that family rankings invert, but not whether the covariance structure itself is responsible, so we estimated the among-family covariance matrix G for all nine traits at once, separately within each environment (Materials and methods, §7.4; S17 and S18 Tables). Two quantities summarize such a matrix. Its leading eigenvector, g_max, is the single weighted combination of the nine traits along which families differ most; the effective number of dimensions, n_D, is the sum of the eigenvalues divided by the largest, and is 1 when all among-family variance lies on g_max alone and 9 when it is spread evenly across the traits. Under acidification, the architecture is close to one-dimensional: accounts for 57% of the among-family variance (95% CI 49%–66%; = 1.77) and exceeds its design-based null (P = 5 × 10⁻⁴; Fig 2f). Nineteen of the 20 cross-stage correlations are negative, with a mean of −0.375 (95% CI −0.50 to −0.20; P = 5 × 10⁻⁴; Fig 2g). Under ambient conditions, the same 51 families show no such structure: the mean cross-stage correlation is −0.001 (P = 0.49; Fig 2h), accounts for only 37% against a null expectation of 35%, and variance is spread over more dimensions ( = 2.70). Under acidification loads positively on all five larval traits and negatively on all four post-metamorphic traits, so the dominant axis of family variation is the cross-stage conflict itself, whereas the ambient loads on the post-metamorphic traits without a consistent sign (−0.22 for metabolic efficiency, +0.29 for gonadal index; Fig 2i), so it describes no cross-stage conflict. The contrast is well supported ( = −0.374, 95% CI −0.56 to −0.16), and the two matrices differ in orientation as well as dimensionality ( 42.5° apart; Krzanowski’s common-subspace statistic 2.28 of a possible 4; mean random-skewers vector correlation 0.66). Acidification does not simply depress performance at both stages—it converts an architecture in which the two stages vary independently into one in which they oppose each other.

The consequence for selection follows directly. Under acidification, the predicted multivariate response to selection on larval fitness, , is positive across the larval traits but negative across the post-metamorphic traits (summed post-metamorphic response −1.91, 95% CI −2.57 to −1.20; S6a Fig; S19 Table), and it points 100° away from the direction that post-metamorphic selection favors—more than a right angle, so the correlated response actively opposes adult-stage selection rather than merely failing to help it. The reciprocal holds (111°), whereas under ambient conditions the same calculation yields no cross-stage cost at all (+0.09; 88°). The conflict axis is not an artifact of nonadditive or maternal effects: among-sire covariances, which isolate the additive component, give a 16.6° from the family-level one (S6b Fig), explaining a similar share of variance (58%) and carrying the same opposing loadings (mean cross-stage r = −0.46). An unweighted principal-component decomposition of the nine tolerance ratios recovers the same axis as the a priori weighted composites (r = +0.95 and −0.70; S6c Fig), confirming that the weighting in Eqs. 2, 3 is not the source of the result.

Single-stage selection enriches hidden liabilities across the life cycle

Finally, we asked whether the observed sign of cross-stage covariance is consequential for adaptation, not merely descriptive of present-day family variance. We projected family frequencies for 20 generations using a relative-fitness family-frequency model with Wright–Fisher mean-fitness normalization (Materials and methods, Eqs. 5–8; S20 Table). Because these projections address evolutionary rather than purely plastic responses in a changing ocean [37,38], the model is not an abundance forecast; it isolates how the sign of cross-stage covariance changes evolutionary direction under the same family pool and the same forcing schedule. Under positive and near-zero scenario surfaces, larval-dominant, balanced, and adult-dominant selection all moved trajectories in broadly adaptive directions, differing mainly in which stage improved faster (Fig 6a, 6b, 6d, 6e, 6g, and 6h). Under the empirical conflict surface, however, stage weighting changed direction: larval-dominant selection improved larval tolerance while eroding adult tolerance, adult-dominant selection produced the inverse, and balanced weighting was the only regime that kept both axes net-positive while retaining the most effective family diversity (Fig 6c,6f,6i). Thus, the same family pool carries adaptive information in opposite directions depending on which stage is screened.

At generation 20, paired family-cluster bootstraps showed that balanced weighting exceeded larval-dominant selection with probability 0.81 and adult-dominant selection with probability 0.87 under the empirical conflict surface (Fig 6j–6l; Eq. 9). These are posterior-style probabilities from paired resamples, not significance tests: 0.81 means that in 81% of bootstrap resamples of the family pool balanced weighting retained more diversity than larval-dominant weighting, which is suggestive rather than decisive, whereas the corresponding probabilities under the positive and near-zero surfaces are 0.33 and 0.47—indistinguishable from a coin toss. The informative comparison is therefore between surfaces, not the absolute value of any one probability. In absolute terms, every regime loses most of the initial diversity by generation 20, because the model deliberately contains only selection: it has no mutation, no recombination among families and no immigration, and families are treated as indivisible units. It is a comparator of evolutionary direction under a common forcing schedule, not a forecast of standing variation, and the differences between regimes should be read as relative rather than absolute—balanced weighting retained ~80% more effective family diversity than larval-dominant weighting and about 3.7-fold more than adult-dominant weighting. Because recombination and mutation would replenish variation in all three regimes, incorporating them would be expected to raise every trajectory without reversing their order, but this remains to be tested. Continuous stage-weight scans placed the diversity maximum near balanced selection ( ≈ 0.45; Materials and methods, §8.4), and balanced weighting retained the most diversity in all 36 combinations of forcing intensity, ramp and horizon (S7 and S8 Figs; S21 and S22 Tables). These projections turn the empirical ranking reversal into an operational warning: single-stage selection does not merely omit information from the other stage; when metamorphosis reverses family rankings, it actively enriches hidden liabilities elsewhere in the life cycle.

1. Study organism and broodstock conditioning

Adult Mercenaria mercenaria (shell length 62–81 mm; n = 130) were collected from a single intertidal–subtidal flat in Long Island Sound, Connecticut, USA, in order to minimize confounding by population structure and to retain a common evolutionary background across the pedigree. Animals were transported to the hatchery and acclimated for 21 d in flow-through, UV-treated seawater (20 ± 0.5 °C; 30 ± 0.5 psu; pH_T = 8.05 ± 0.02), and fed ad libitum with Isochrysis galbana and Chaetoceros calcitrans at a 1:1 ratio (~5 × 10⁴ cells mL⁻¹). Shell length, total wet mass, and sex were recorded before crossing.

Spawning was induced by thermal stimulation and was carried out within a restricted temporal window to limit uncontrolled spawning-order effects. Broodstock were strip-spawned only when necessary to complete missing crosses. Embryos from each cross were rinsed through 40-µm mesh, checked for fertilization success after first cleavage, and retained only when cleavage exceeded 80%. This procedure generated the operational family units that were carried consistently through larval assays, juvenile assays, source-data worksheets and pedigree-aware models.

2. Quantitative-genetic family design

The pedigree was intentionally partial rather than fully factorial [30]. Each of 12 sires was crossed to a subset of the 6 dams, yielding 51 operational families (S1 Fig). This design sacrifices some balance relative to a complete diallel but preserves the paternal half-sib structure needed for quantitative-genetic inference while remaining practical for a nonmodel bivalve spawning programme. Operational families were mapped to an explicit sire–dam table and then carried through every analysis in one of three ways: as family identifiers, as nested pedigree terms within mixed models, and as the numerator-relationship structure used for the supplementary animal model [30,50]. Families sharing a sire are not independent evolutionary units, so retaining the pedigree in all downstream models prevents spurious certainty in family rankings and provides a conservative framework for evaluating cross-stage covariance.

Because the pedigree is partial and the phenotypic inputs used for the bivariate animal model are family-mean composite scores rather than repeated individual records, the quantitative-genetic analysis is interpreted as a formal support for covariance sign and mode rather than as a definitive estimator of all narrow-sense parameters. This boundary is made explicit in the main manuscript.

3. Seawater chemistry manipulation and characterization

Ambient and acidified systems were maintained independently and continuously monitored following established ocean acidification best-practice workflows [51]. Ambient seawater targeted present-day local conditions, whereas the acidified treatment targeted end-century carbonate chemistry under a high-emissions trajectory (pH_T ≈ 7.70, Ω_arag ≈ 1.0). Independent treatment loops were used to avoid pseudo-replication at the chemistry level and to preserve treatment-specific carbonate-system variance for downstream mixed-model inference.

pH was logged continuously with NBS-calibrated electrodes and independently checked by benchtop spectrophotometric measurements. Total alkalinity was measured regularly by open-cell titration, and carbonate-system variables were recalculated from pH and total alkalinity with CO2SYS [52], following best-practice carbonate chemistry reporting [53]. Temperature and salinity were recorded daily. Treatment separation and its stability over the exposure are shown in S2 Fig, and the full carbonate-system record is summarized in S1 Table. Treatment systems were serviced separately to avoid gas-exchange cross-talk, and water-exchange rates were adjusted to keep dissolved oxygen, temperature, and salinity tightly overlapping between treatments.

4. Larval experiment: Hierarchical design and phenotyping

The larval assay was designed to separate treatment effects from temporal and vessel effects. Two temporal blocks were staggered by 5 d. Each family was reared in its own independently plumbed 20-L polycarbonate tank within each treatment × block combination (51 families × 2 treatments × 2 blocks = 204 tanks), and each tank contained three 250-mL larval vials. Each vial received exactly N₀ = 180 larvae from a single family, counted individually under stereomicroscopy, giving 51 × 2 × 2 × 3 = 612 vials in total (S1 Data). Larval and post-metamorphic exposures were entirely independent of one another: no animal experienced acidification at both stages. Larvae used in the larval assay were terminated at that assay’s endpoint, and the juveniles entering the post-metamorphic experiment were drawn from separate ambient-reared cohorts of the same 51 families and then randomly assigned to the ambient or acidified nursery. Larval and post-metamorphic tolerance are therefore linked only through shared family (and hence pedigree) membership, which is precisely the level at which selection operates and at which cross-stage covariance is identifiable; no result in this paper depends on within-individual carry-over.

Larvae were reared under constant gentle aeration and fed every 48 h with I. galbana plus C. calcitrans with a 1:1 mixture at 10⁴ cells mL⁻¹. Vials were terminated when at least 50% of survivors reached the pediveliger stage, or at day 21 if that threshold had not yet been reached. All surviving larvae were counted, and dead individuals were identified by collapsed or darkened tissue. From each vial, up to 20 survivors were sampled for shell area (calibrated image analysis), shell thickness (cross-sectional SEM), nanoindentation hardness (Hysitron TI950; Berkovich probe; 500 µN peak load; 5-s hold) and days to metamorphosis. Image analysis, indentation, and SEM scoring were performed blind to family and treatment identity.

5. Adult experiment: Phenotyping and delayed reproductive signature

The adult phase was designed to capture both acute post-metamorphic maintenance costs and a delayed reproductive signature of earlier OA exposure. Therefore, our designation of “adult” stage includes reproductively active adults as well as post-metamorphic juveniles. Five juveniles per family were assigned to each treatment (n = 510 total individuals) and held individually in partitioned chambers nested within flow-through treatment tanks during a 90-d OA nursery phase.

Shell growth rate was calculated from repeated caliper measurements of maximum shell length over the nursery interval and is reported as mm per 30-d period. Shell dissolution and shell growth are not alternatives: in bivalves the periostracum-covered growing margin continues to accrete new carbonate while the older, abraded umbonal and lateral surfaces—where the periostracum is thin or worn—are exposed to undersaturated water and corrode. Net linear growth and localized surface dissolution therefore proceed simultaneously on the same shell, and both were scored. Dissolution was quantified on the external shell surface as a 0–4 ordinal scale (0 = intact periostracum and no visible etching; 1 = localized loss of surface luster; 2 = visible etching or chalky patches confined to the umbo; 3 = etching extending over the lateral surfaces with exposed prismatic layer; and 4 = extensive pitting or perforation of the outer layer) by two trained scorers working independently and blind to family and treatment identity (Cohen’s κ = 0.83); disagreements were resolved by consensus. Oxygen consumption (MO₂)—the mass-specific rate of oxygen uptake, used here as an index of routine maintenance metabolism—was measured by sealed-chamber respirometry using a YSI ProODO optical O₂ logger during 3-h incubations at days 30, 60, and 90; rates were standardized to dry mass. Nearly all juveniles survived the 90-day nursery (survival rate ≈ 1), so post-metamorphic fitness variation among families is expressed through growth, shell maintenance, metabolic cost and subsequent reproductive output rather than through differential mortality.

After the OA nursery phase, all juveniles were transferred back to ambient pond culture. Gonadal index was then quantified after two years of common grow-out as (gonad wet mass/total soft-tissue wet mass) × 100%, so reproductive status is interpreted as a delayed signature of earlier OA exposure rather than a trait measured concurrently with the 90-d maintenance phenotypes. This timing distinction is preserved throughout the analysis and the discussion.

6. Statistical analyses

For each trait, we first screened biologically justified random-effect structures and then estimated treatment effects using the most conservative supported structure. Candidate random terms followed the mating design and the experiment hierarchy: sire, dam, operational family, and block for larval traits; sire, dam, and operational family for adult traits. The objective was not to maximize model complexity but to preserve pedigree and block structure whenever it remained estimable and informative. All analyses were implemented in Python with pandas, numpy, scipy, statsmodels, openpyxl and matplotlib, following mixed-model logic commonly used in lme4 and glmmTMB workflows [54,55].

7. Cross-stage family association and composite tolerance

Cross-stage analyses were based on stage-specific composite retention scores rather than on raw trait deltas, so that resilience is captured as stage-wide performance under stress rather than success in any one endpoint. For each family , larval composite tolerance and adult composite tolerance were calculated as weighted geometric means of trait-retention ratios: Family composite scores, rank-reversal diagnostics, and robustness checks are summarized in S3–S5 Tables.

(2)

The multiplicative form preserves multiplicative penalties: a family cannot appear highly tolerant overall if it performs very poorly on any one key component trait. The equivalent log-form used as the regression target is

(3)

Default stage-specific weights were assigned a priori according to each trait’s expected proximity to fitness. Larval tolerance weighted survival most strongly (0.60), followed by metamorphosis rate (0.20), because these traits define passage through the pelagic demographic bottleneck; shell area, thickness, and hardness received smaller weights (0.08, 0.07, and 0.05) as mechanistic components of early calcification. Post-metamorphic tolerance weighted growth (0.45) and delayed gonadal index (0.35) most strongly as proxies for somatic performance and reproductive allocation, with dissolution resistance and metabolic efficiency retained as supporting maintenance-cost traits (0.10 each). Five alternative weighting schemes produced the same negative cross-stage association, indicating that the result was not driven by the default weighting choice (S5 Table; S3 Fig).

8. Family-frequency evolutionary projection

To project evolutionary trajectories under end-century OA, we implemented a relative-fitness family-frequency model [63] in which each of the N = 51 families occupies a point in larval–adult tolerance space and frequencies are updated by Wright–Fisher mean-fitness normalization [63]. Projection scenarios differed in two respects: the geometry of the cross-stage landscape (Model I, positive association; Model II, near-zero association; Model III, empirical conflict) and the relative weighting assigned to larval versus adult performance. Projection outputs are tabulated in S20–S22 Tables.

Initial frequencies are uniform, . The forcing ramp is

(5)

with = 0.22 and = 0.02 per generation as defaults. Under stage weights (, ) constrained to , the family-specific selection coefficient is

(6)

Here and are family f’s larval and post-metamorphic composite tolerance scores from Eqs. 2, 3, each standardized across the 51 families to zero mean and unit variance: a value of zero denotes a family of average tolerance for its stage, positive values a family more tolerant than average, and one unit is one among-family standard deviation. The stage weights and set how much each stage contributes to selection and are constrained to sum to one.

Relative fitness is and the family-frequency update follows Wright–Fisher mean-fitness normalization:

(7)

Population-level outcomes reported at each generation are the frequency-weighted means of the stage-specific composites, and —the quantities reported for each generation in the Results, so that 0 is the starting population mean and one unit is one among-family standard deviation—together with the effective family diversity :

(8)

S1 Fig. Experimental design and measurement timing. a, Life cycle of the hard clam, separating the planktonic larval stage from the benthic post-metamorphic stage. b, Schematic of broodstock sourcing, the partial half-sib pedigree, the larval vial-level assay, the post-metamorphic ocean acidification (OA) nursery, and the subsequent two-year common-garden ambient pond grow-out.

https://doi.org/10.1371/journal.pbio.3004079.s001

(TIF)

S2 Fig. Seawater chemistry fidelity and treatment stability.

a, Time series of pH_T. b, Aragonite saturation state Ω_arag. c, Temperature. d, Salinity. Dotted horizontal lines in a and b indicate treatment targets. Treatment separation was stable throughout the exposure, drift was minimal, and the OA manipulation was imposed under tightly controlled physicochemical conditions. The numerical data underlying this figure are provided at https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s002

(TIF)

S3 Fig. The negative cross-stage association is robust to tolerance definitions and pedigree groups.

a, Pearson r and Spearman ρ across the baseline and five alternative composite-scoring schemes; the negative sign is invariant to reasonable redefinitions. b, Equal-weight log-composite scatter. c, Threshold-quantile sweep showing the fraction of top-quantile families in each stage that fall below the opposite-stage median across quantile cut-offs. d, Pearson r after leave-one-out removal of each sire and each dam, relative to the full-data estimate (dashed vertical line). The numerical data underlying this figure are provided in https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s003

(TIF)

S4 Fig. Quantitative-genetic summaries support a negative cross-stage covariance.

a, Sire-level mean larval versus adult composite tolerance for the 12 shared sires (Pearson r = −0.66, 95% CI [−0.90, −0.14], P = 0.019). b, Summary of covariance-sign diagnostics across family-level, pedigree-conditioned, sire- and dam-level, and animal model estimates; every summary retains a negative sign. c, Restricted-likelihood profile across candidate rG values; the dashed horizontal line marks the 2Δ_REML = 3.84 threshold. The numerical data underlying this figure are provided in https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s004

(TIF)

S5 Fig. Family reaction norms reveal the stage-specific redistribution of OA performance.

a–d, Ambient-to-acidified reaction norms for 51 families across representative larval traits (a, survival; b, shell area) and adult traits (c, growth; d, delayed gonadal index). Families highlighted in red in a and b are the 10 top-quartile larval families that drop below the post-metamorphic median; families highlighted in red in c and d are the 11 top-quartile post-metamorphic families that drop below the larval median. These are the two sets ringed in red in Fig 2b, so the highlighted families differ between the left-hand (larval) and right-hand (post-metamorphic) panels. These reaction-norm panels show the stage-specific redistribution of performance that underlies the cross-stage ranking reversal. The numerical data underlying this figure are provided at https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s005

(TIF)

S6 Fig. Supporting detail for the cross-stage covariance analysis.

a, Predicted multivariate response to selection on larval fitness, , summed within each stage, under the acidified and the ambient matrix; vertical bars are 95% family-cluster bootstrap intervals (B= 3,000 resamples of the 51 families), and P(<0) is the bootstrap probability that the post-metamorphic sum is negative. b, Trait loadings on the leading eigenvector of the among-sire (paternal half-sib) covariance matrix against the family-level loadings; the two axes are 16.6° apart, so the conflict axis is present in the additive component. c, Family scores on an unweighted principal component of the nine trait-specific tolerance ratios against the a priori weighted composites, with least-squares fits (r = +0.95 and −0.70), showing that the weighting scheme of Eqs. 2 and 3 is not the source of the result. The numerical data underlying this figure are provided at https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s006

(TIF)

S7 Fig. Projected outcomes are sensitive to continuous stage weighting.

a, Generation-20 mean larval tolerance across the full larval-weight gradient. b, Generation-20 mean adult tolerance. c, Effective family diversity. Under the empirical conflict surface, the diversity maximum is located near balanced weighting ( ≈ 0.45) and drops off on both sides, confirming that the superiority of balanced weighting does not depend on the three anchor scenarios displayed in Fig 6. The numerical data underlying this figure are provided in https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s007

(TIF)

S8 Fig. Projection outcomes are robust to forcing design.

a–c, Generation-20 effective family diversity against generation horizon under a base forcing intensity of α₀ = 0.16 (a), 0.22 (b) and 0.28 (c), for each of the three stage weightings. Line color encodes the weighting; line style encodes the per-generation ramp β. Absolute values vary with forcing design, but balanced weighting retains the most family diversity in all 36 combinations of α₀, β, and horizon (S22 Table). The numerical data underlying this figure are provided in https://doi.org/10.6084/m9.figshare.34030317.

https://doi.org/10.1371/journal.pbio.3004079.s008

(TIF)

S3 Table. Composite tolerance scores, cross-stage family association and rank reversal.

Stage-specific OA tolerance is represented as a weighted geometric mean of the dominant demographic components for each life stage. Cross-stage ranking is based on these composite scores rather than on a single trait; robustness analyses are summarized in S4–S6 Tables.

https://doi.org/10.1371/journal.pbio.3004079.s011

(XLSX)

S6 Table. Bivariate pedigree animal model and sire/dam covariance summaries.

Model: y = Xβ + a + e, with a ~ N(0, A ⊗ G) and e ~ N(0, I ⊗ R). A is the numerator relationship matrix derived from the partial half-sib pedigree. Because phenotypes are family means on a partial pedigree, this table supports covariance-sign inference rather than serving as the sole primary effect-size estimate.

https://doi.org/10.1371/journal.pbio.3004079.s014

(XLSX)

S12 Table. Selection on larval shell area and developmental time.

Standardized selection differentials S′ for shell area and days to metamorphosis, with the retained quadratic term for days to metamorphosis, plotted in Fig 4e. Each value is an ordinary-least-squares regression of vial-level relative survival on one standardized trait, fitted separately within each treatment, with traits standardized across the pooled data and relative fitness computed within treatment; the multivariate selection gradients β used for the predicted responses are reported in S19 Table. Projection outputs are reported separately in S20–S22 Tables.

https://doi.org/10.1371/journal.pbio.3004079.s020

(XLSX)

S17 Table. Among-family covariance (correlation-scale) matrix G for nine traits, estimated separately within the acidified and the ambient environment.

Correlation-scale matrices for the five larval and four post-metamorphic traits described in Materials and methods §7.4, with eigenvalues, effective dimensionality, family-cluster bootstrap intervals and P values against the design-based null in which cross-stage covariance is set to zero. Traits were standardized to unit variance before the among-family decomposition.

https://doi.org/10.1371/journal.pbio.3004079.s025

(XLSX)

S18 Table. Trait loadings on the leading eigenvector of the among-family covariance matrix.

Loadings under acidification and under ambient conditions at the family-level, together with the among-sire (paternal half-sib) loadings used to check that the conflict axis is present in the additive component.

https://doi.org/10.1371/journal.pbio.3004079.s026

(XLSX)

S19 Table. Selection geometry on the estimated among-family covariance matrix.

Stage-specific standardized selection gradients, the predicted multivariate responses summed within each stage with family-cluster bootstrap intervals, the angles between predicted responses and stage-specific selection gradients, and evolvabilities. The matrix-comparison statistics contrasting the acidified and ambient matrices are given in S17 Table.

https://doi.org/10.1371/journal.pbio.3004079.s027

(XLSX)

S23 Table. Among-family variance components and the eigen-analysis repeated on the correlation scale.

Trait-by-trait among-family variance in each environment, with the leading-axis share, effective dimensionality and g_max loadings recomputed from the correlation matrices of S17 Table, and a cross-substitution of correlation structure and variance scale between the two environments. Because every trait is standardized to unit total variance within its environment, each variance component is that trait’s among-family share and their sum is the overall among-family variation.

https://doi.org/10.1371/journal.pbio.3004079.s031

(XLSX)

S2 Data. Adult individual-level source data.

Each row is one post-metamorphic individual. FamilyID, Sire, and Dam are merged from the partial half-sib pedigree scaffold so every record can be traced to the mating design and to the source data used for Figs 2, 5, and 6. Shell growth, dissolution, and MO₂ derive from the 90-day OA nursery phase; gonadal index was quantified after all juveniles were returned to ambient seawater pond culture for 2 years.

https://doi.org/10.1371/journal.pbio.3004079.s033

(XLSX)