Dear editor, this document contains a point by point reply to the issues raised by the reviewers (in red); changes made in the MS are indicated in blue. All original comments by the reviewers are left in black font. You can see the colors in the attached Supplement PDF file. See also the Supplement PDF file on author reply for the Specific comments.

Referee comment from Anonymous Referee #1

Summary

In this research article, Kunio Kaiho presents novel findings on the future development of metazoan diversity in superterranean, subterranean, surface-water, and deep-water habitats based on diversity changes in the past. By incorporating seven different environmental drivers, the author projects the com-plete extinction of metazoans within the next 700 million years, which is 300–400 million years earlier than previously estimated.

General comments

Overall, the manuscript is well written and provides novel insights into an important field of research. The language is almost perfect, clear, and easy to follow. However, there are a few general points that should be addressed before final publication of the article.

Neither the Introduction nor the Discussion provides much context regarding previous research efforts. While the Introduction nicely explains the different environmental drivers incorporated into the current study, it is unclear what previous research entailed and what the current study adds to it. These aspects should be included in the revised manuscript.

Author Reply: I agree with the comments. I revised words marked in blue in the attached manuscript (Kaiho Revise Marked 1).

Added lines 38-43 in the Introduction. Estimates for the end of Earth’s biosphere published after 2010 vary widely. Projections based on surface temperature range from 1.0 to 5.0 Gyr (O’Malley-James et al., 2012; Rushby, 2013; Leconte et al., 2013; Wolf and Toon, 2015), while scenarios based on CO₂ depletion yield estimates of 0.84–1.08 Gyr (Rushby, 2015; Ozaki and Reinhard, 2021). Mello and Friaça (2019) suggest that biosphere collapse is unlikely before 1.5 Gyr based on thermal constraints. However, a decline in atmospheric oxygen to 1% PAL within 1.08 ± 0.14 Gyr (1σ) as predicted by Ozaki and Reinhard (2021) may lead to the earlier extinction of metazoan life.

Added lines 62-63 in the Introduction. This study builds upon previous models by integrating anthropogenic crises (Waters et al., 2011; Ceballos et al., 2015; Waters et al., 2016; Kaiho, 2022, 2023), cyclical climate rhythms, and abrupt climate events. It considers seven key factors—anthropogenic crises, long-term warming, cyclical climate rhythms, abrupt events, C₃ plant collapse, C₄ plant decline, and oxygen depletion—to project metazoan extinction over the next 1.5 Gyr. Projections are grounded in temperature and oxygen modeling, thermal tolerance limits, and observed metazoan diversity trends.

Similarly, the Discussion repeats the major results of the current study without discussing them in the context of previous findings. For example, it is repeatedly mentioned throughout the manuscript that the current study projects metazoan extinction to occur 300–400 million years earlier than previous estimates, but these previous estimates are not further specified. What differences between previous studies and the current study may cause these different results? Why are the results of the current study more/similarly realistic? These questions should be addressed in the Discussion.

Author Reply: Added 3.2 (Results) and 4.1 (Discussion) sections.

3.2  Metazoan lifetime estimation under four scenarios

When only the long-term warming trend driven by the gradual increase in solar luminosity is considered (Mello and Friaça, 2019), metazoans are projected to go extinct at approximately 1.3 Gyr, based on the intersection of the black dashed line with the upper boundary of GATEU90 in Figure 1. Under this scenario, surface-dwelling metazoans are expected to go extinct slightly earlier, at 1.2 Gyr, based on the same trend intersecting the upper boundary of GATES90.

When long-term cyclical fluctuations between icehouse and greenhouse phases are incorporated into the model, extinction is projected to occur at 1.2 Gyr, corresponding to the intersection of the orange line with the top of GATEU90. In this case, surface metazoans are expected to disappear by 1.0 Gyr, as indicated by the intersection with GATES90.

Incorporating average surface temperature anomalies associated with past mass extinction events further lowers the projected extinction time to 1.0 Gyr, based on the red circle's intersection with the upper boundary of GATEU90. However, actual complete extinction is expected to occur earlier, between 0.7 and 0.8 Gyr, due to compounded survival rate reductions (0.01–0.1) caused by food scarcity (FS) and oceanic anoxia. These stressors are expected to be triggered by elevated temperatures during abrupt extinction events (Events 8–10), which involve the collapse of surface metazoan populations and severe reductions in primary productivity caused by global-scale extreme warming (see Table A5).

By 0.7 Gyr, atmospheric oxygen and CO₂ levels are not anticipated to be the dominant extinction drivers. Instead, the primary cause of complete extinction is projected to be extreme surface warming, resulting from the combined effects of increased solar luminosity, long-term climatic oscillations, and large-scale volcanic activity. Therefore, the final extinction of all metazoan life is projected to occur at 0.7 ± 0.05 Gyr.

4  Discussion 4.1  Lifetime estimation

Previous studies have proposed varying estimates for the remaining lifespan of metazoan life on Earth. A 1.2 Gyr estimate has been cited as a plausible median value (Jebari and Sandberg, 2022), while 1.3 Gyr is projected based solely on the long-term warming trend driven by increasing solar luminosity (Fig. 1), and 1.1 Gyr corresponds to the point at which atmospheric oxygen is predicted to decline to 1% PAL (Ozaki and Reinhard, 2021).

In contrast, the present study estimates the total remaining metazoan lifespan at approximately 0.7–0.8 Gyr from now, which is 300 to 600 million years earlier than previous estimates. This discrepancy arises because earlier models do not account for the compounding effects of long-term cyclical icehouse–greenhouse climate phases, average surface temperature anomalies during mass extinction events, and cooling effects linked to mantle temperature decline.

The estimates from the current study are considered more realistic, as they incorporate all major drivers of surface temperature variability—including long-term trends, cyclical oscillations, and abrupt catastrophic events—providing a more comprehensive projection of the environmental conditions leading to metazoan extinction.

Author Reply: Revised to “Ultimately, a final extinction event—likely initiated by large-scale volcanism—will be primarily driven by extreme global warming.” in lines 16-17 (Abstract).

Added “, primarily driven by global warming” in line 704-705 (Conclusions). The sentence is “This scenario indicates that complete metazoan extinction is expected to occur within 0.7–0.8 Gyr, primarily driven by global warming.”

In addition, I think that some parts of the Methods section are difficult to follow. Firstly, this section uses many abbreviations, but not all of them are defined in the text itself, only in figure/table captions (e.g., PAL is only defined in the caption of Fig. 1). Secondly, many terms are unclear to the reader and require further explanation (e.g., what exactly are diversity rates and what is the difference between survival rates and survival area rates?). Thirdly, the argumentation is partly difficult to follow since the required explanations are either insufficient or provided later in the Results or Discussion section. I recommend adding further explanations and revising the structure of the manuscript where necessary. I give specific examples in the “Specific comments” section.

Author Reply: I agree with the comments. Revised the 2.7 section in Methods: Future metazoan diversity estimation as the following section (revised parts are marked in blue). Also revised Figures 5–7 based on the revised 2.7 section.

2.7  Future metazoan diversity estimation Future changes in metazoan diversity are influenced by the ongoing anthropogenic crisis ("event 0"), subsequent mass extinction events (events 1–11), C₃ and C₄ plant crises, and gradual oxygen depletion (see Table 2). The projected diversity of insects, terrestrial tetrapods, and marine metazoans is estimated before extinction events, immediately after, and following recovery using the equations below: Diversity loss due to extinction event: D_t = D_t-1 × SRC × FSR                                (12)

Diversity following recovery (after 50 Myr from the extinction event to before the next event):

D_t+1 = D_t + (D_t-1 – D_t)× RR                             (13)

Recovery Rate (RR):

RR = RRW × RRO × RRP                             (14)

In these equations, Dt​ represents metazoan diversity at time step t, corresponding to the level immediately following an extinction event. Dt−1​ denotes the diversity prior to the extinction event, while Dt+1​ reflects the diversity after the recovery phase, measured at the midpoint between extinction events. SRC is the Survival Rate associated with climate-driven crises, including mass extinctions and C₃–C₄ plant collapses. FSR is the Food Scarcity Rate, reflecting the impact of the collapse of plants and primary producers. The total Recovery Rate (RR) is calculated as the product of three components: RRW (recovery from gradual warming), RRO (recovery from progressive oxygen decline), and RRP (recovery from decreased primary productivity due to CO₂ reduction).

These equations are applied sequentially across time steps from event 0 through event 16, encompassing extinction episodes, recovery phases that conclude at the midpoint between events, and the subsequent interval leading up to the next extinction event (see Table A5).

2.7.1  Survival Rate associated with Climate change (SRC)

Survival Area Rate (SAR) is rate of land and ocean area where metazoans survive in all land and ocean area (km2/km2). When extinction occurred in 0–10, 0–20, 0–30, 0–40, 0–50, 0–60, 0–70, 0–80, and 0–90° latitudes by warming, SAR values are defined as 0.83, 0.66, 0.50, 0.36, 0.24, 0.14, 0.06, 0.02, 0.00, respectively, under the same rate of land and ocean in those latitudes. The rates SAR are decided by only temperature 46 °C using Figure 2. SAR_S is Survival Area Rate for St and Sw metazoans, SAR_U and SAR_D are Survival Area Rate for U and D metazoans, respectively. These SAR are obtained from GATES, GATEU, and GATED in Figure 2.

The SRC (survival rate by climates) is calculated as:

SRC_T = 0.95 × SAR_S + 0.05 × SAR_U                                          (15)

SRC_M = 0.67 × SAR_S + 0.33 × SAR_D                                          (16)

Here, SRC_T and SRC_M represent the terrestrial and marine SRC, respectively. The coefficient 0.05 corresponds to the proportion of subterranean metazoan families among all terrestrial metazoan families (15 out of 315), based on mammalian lineage data (Recknagel & Trontelj, 2021; Benton, 2010). The remaining 0.95 represents superterranean (surface-dwelling) taxa. The coefficient 0.33 reflects the proportion of deep-sea fish families among all marine fish families, based on an estimated 6% of teleost species restricted to depths >200 m (Miller et al., 2022), with scaling applied via the species-genus-family extinction relationship from Kaiho (2022). The remaining 0.67 applies to surface-dwelling marine taxa. Equations (15) and (16) are used for modeling extinction scenarios during Events 5–16.

SAR_D is influenced by both temperature and dissolved oxygen levels. Elevated surface temperatures reduce oxygen concentrations in deep water, a process linked to deep-sea extinctions during the end-Permian and end-Cenomanian anoxia–euxinia events, despite elevated atmospheric O₂ (e.g., Sun et al., 2012; Kaiho et al., 2013, 2016a). High surface temperatures can cause extinction in both surface and deep-water taxa. Although deep water temperatures are lower than those at the surface, the greatest thermal anomalies occur in surface waters, while deep-water temperatures remain relatively constant throughout the water column. Consequently, SAR_D is assumed to approximate SAR_S.

Thus, for Events 5–6, 8–16, and all non-events after Event 8 (excluding the interval between Events 9 and 10), SAR_D is set equal to SAR_S, as these intervals are characterized by global surface temperatures comparable to or exceeding those of the end-Permian.

The dominant climate driver for mass extinction varies by event. Events –5 to 4 involve both warming and cooling phases, while Events 5–16 are exclusively warming-driven, corresponding to the yellow–orange shaded zone in Figure 1.

For Event 0 (the Anthropogenic Crisis), the maximum SRC values are set at 0.95 for insects, 0.70 for terrestrial tetrapods, and 0.90 for marine metazoans. These values reflect a worst-case scenario involving full-scale nuclear war, combined with moderate anthropogenic pollution, deforestation, and global warming (Kaiho, 2023). In the absence of nuclear conflict, SRC is assumed to be 1.0 for all groups. For Events 1–4 and 7, SRC values are 0.81 for insects, 0.63 for terrestrial tetrapods, and 0.74 for marine metazoans, based on average extinction percentages reported in Tables A3 and A4.

In events 5 and 6 where global average surface temperatures reach 38–39 °C in Figures 2a and 2c. The 38 °C and 39 °C correspond to GATES40 and 50 (SAR_S: 0.36 and 0.24) in Figures 2a and 2c, and GATEU00 and GATEU10 (SAR_U: 1.00 and 0.83) in Figure 2b. The both temperatures are lower than GATED showing 48°C (Fig 2d). Using equations 15 and 16, the SRC values are:

In Events 5 and 6, global mean surface temperatures reach 38–39 °C (Figures 2a, 2c). These correspond to GATES40 and GATES50 (SARS: 0.36 and 0.24), and to GATEU00 and GATEU10 (SARU: 1.00 and 0.83) in Figure 2b. Both values are below the stable GATED threshold of 48 °C (Figure 2d). Applying equations (15) and (16):

Event 5:

SRC_T = 0.95 × 0.36 + 0.05 × 1.00 = 0.39

SRC_M = 0.67 × 0.36 + 0.33 × 1.00 = 0.57

Event 6:

SRC_T = 0.95 × 0.24 + 0.05 × 0.83 = 0.27

SRC_M = 0.67 × 0.24 + 0.33 × 1.00 = 0.49

These values are listed in Table 2.

In Events 8–10, global surface temperatures rise to 43–45 °C (Figures 2a, 2c). At 43 °C, GATES90 is reached, resulting in complete extinction of surface-dwelling metazoans (SARS = 0). The corresponding SARU values, 0.50 and 0.30, are based on GATEU30 and GATEU45 (Figure 2b). GATED remains at 48 °C, so SARD = 0. Using equations (15) and (16):

Event 8:

SRC_T = 0.95 × 0 + 0.05 × 0.50 = 0.025

Events 9 and 10:

SRC_T = 0.95 × 0 + 0.05 × 0.30 = 0.015 

Events 8–10:                                 

SRC_M = 0.67 × 0 + 0.33 × 0 = 0

These SRC values are also summarized in Table 2.                   

2.7.2 Food Scarcity Rate (FSR)

In addition to direct climatic impacts, food scarcity significantly contributes to extinction risk. As plants and primary producers collapse, only organisms capable of surviving on bacterial biomass or sedimentary organic matter—along with their predators—will remain. Consequently, an additional decline in survival rate, quantified as the Food Scarcity Rate (FSR), is expected during abrupt extinction events such as Events 8–10.

These events are characterized by the extinction of superterranean and surface-water (SS) metazoans and severe reductions in primary productivity due to sunlight loss—conditions common in major mass extinction scenarios (see Fig. 1). A low FSR value of 0.01 reflects survival through alternative nutritional pathways, including hydrothermal vent ecosystems and bacterial-based underground food sources (Cosson and Soldati, 2008; Miroshnichenko, 2004; Kelley et al., 2005).

Conversely, a high FSR value of 0.1 represents an optimistic estimate, assuming evolutionary adaptation of primary producers to extreme temperatures, allowing some limited ecosystem function to persist. This range (0.01–0.1) is applied to adjust survival estimates in scenarios where abrupt collapse of food webs occurs due to light inhibition and temperature stress.

2.7.3  Recovery Rate by Warming(RRW)

The Recovery Rate by gradual Warming (RRW), applied outside of abrupt extinction events, is calculated using the same structure as Equations (15) and (16), but based on temperature anomalies associated with gradual climate changes. These temperature anomalies are derived from Figure 2.

The RRW is calculated separately for terrestrial and marine metazoans using the following equations:

RRW_T = 0.95 × SAR_S + 0.05 × SAR_U                                            (17)

RRW_M = 0.67 × SAR_S + 0.33 × SAR_D                                           (18)

Here, RRW_T and RRW_M denote the recovery rates for terrestrial and marine metazoans, respectively. The coefficients reflect the relative contributions of surface and subsurface habitats, consistent with SRC calculations in earlier sections.

2.7.4  Recovery Rate by Oxygen decline(RRO)

The Recovery Rate by Oxygen decline (RRO) is calculated separately for terrestrial and marine metazoans, reflecting projected atmospheric O₂ levels over the next 1.0 Gyr. The RRO varies across time intervals and is defined by the following equations:

For terrestrial metazoans (RRO_T):

During the next 0.5 Gyr: RRO_T = 1.0                                        (19)

During the next 0.5–0.8 Gyr: RRO_T = 2.67 – 3.33T                             (20)

During the next 0.8–1.0 Gyr: RRO_T = 0.0                                     (21)

For marine metazoans (RRO_M):

During the next 0.5 Gyr: RRO_M = 1.0                                        (22)

During the next 0.5–0.8 Gyr: RRO_M = 2 – 2T                                  (23)

During the next 0.8–1.0 Gyr: RRO_M = 1.84 – 1.80T                             (24)

In these equations, T is the numerical time variable in Gyr, ranging from 0.5 to 1.0. RRO_T and RRO_M represent the recovery rates for terrestrial and marine metazoans, respectively. These values are derived from oxygen level projections by Ozaki and Reinhard (2021).

Atmospheric oxygen levels are projected to steadily decline over the next 1.1 Gyr. Beginning at approximately 1.0 PAL around –0.1 Gyr (i.e., present time), O₂ levels are expected to drop to ~0.5 PAL (median: 0.3–0.7) by 0.5 Gyr. This decline continues to ~0.3 PAL (median: 0.07–0.5) at 0.8 Gyr, followed by a rapid collapse to ~0.01 PAL between 1.0 and 1.1 Gyr.

These projections inform the modeled RRO values and are used in conjunction with observed relationships between atmospheric O₂ levels and metazoan/plant diversity (Fig. 3) to assess the impact of oxygen depletion on biodiversity trajectories throughout Earth's future.

2.7.5  Recovery Rate by Primary productivity(RRP)

A molecular-level investigation of a C₃ plant's response to low CO₂ concentrations (100 ppm compared to the typical 380 ppm) revealed that reduced CO₂ levels lead to a significant decline in biomass productivity (Li et al., 2014). In the future, such low CO₂ conditions are projected to occur at approximately 0.5 Gyr (ranging from 0.4 to 0.65 Gyr), despite rising surface temperatures driven by increasing solar luminosity. This decline in atmospheric CO₂ is expected to cause a gradual reduction in net primary productivity (NPP), ultimately contributing to long-term decreases in metazoan diversity.

Approximately 40 million years after the extinction of C₃ plants, C₄ plants are expected to evolve into tree-like forms. This evolutionary transition is anticipated to support the diversification of metazoans that rely on such vegetation, paralleling the Devonian rise of terrestrial plant ecosystems. These recovering metazoan groups are expected to originate from species formerly dependent on C₃ plants. The Recovery Rate by primary Productivity (RRP) in this context is estimated to range from 0.3 to 0.8, reflecting the evolutionary and ecological challenges in re-establishing complex, tree-supported food webs from C₄ vegetation. This recovery event is centered at 0.4 Gyr, based on a temporal uncertainty of ±0.2 Gyr (Table 2).

Although oceanic primary producers are predominantly phytoplankton, both C₃ and C₄photosynthetic pathways coexist in marine environments (Reinfelder et al., 2000, 2004). To approximate the effect of terrestrial plant crises on marine metazoan diversity, equivalent reduction and recovery values are provisionally applied to marine systems, mirroring those used for terrestrial tetrapods under two modeled scenarios.

For Events 8–10, which involve abrupt primary productivity collapse due to sunlight reduction, RRP is estimated between 0.1 and 0.3, representing an intermediate range between C₃ and C₄ plant crises.

The C₄ plant crisis is expected to occur at approximately 0.97 ± 0.2 Gyr, aligning closely with Event 11. At this stage and beyond, RRP values decline to between 0.01 and 0.1, as NPP is assumed to approach zero. Under such conditions, only UD metazoans—those subsisting on bacteria, detritus, or residing in deep-sea hydrothermal ecosystems—would persist (Cosson and Soldati, 2008; Miroshnichenko, 2004; Kelley et al., 2005). This sharp reduction in RRP reflects the critical dependence of most metazoans on photosynthetically sustained food webs.