The Permian–Triassic paleogeography is derived from PANALESIS (Vérard, 2019; Vérard, 2021), a global plate tectonic model providing maps every 10 Myr from 888 Ma (Tonian) to the present. The PANALESIS paleogeography for the PTB, which we use as a fixed boundary condition in our climate simulations, has proven to be in good agreement with geochemical and paleontological records (Chablais et al., 2011; Peyrotty et al., 2020; Bucur et al., 2020; Le Houedec et al., 2024), particularly in locating elements within the intertropical zone. The location of island arcs, continental ribbons, and even parts of Pangea, such as South China, remains, however, subject to uncertainties, with latitudes potentially varying up to ca. 10 ° . The raw PANALESIS map was adapted to the climate model horizontal resolution of 2.8 ° ( ∼ 300 km). Specifically, seaways narrower than a few pixels were enlarged, while the smallest ones, along with epicontinental seas and lakes, were closed. The resulting topography used in the simulations is shown in Figure 1A.

The study by Nowak et al. (2020) was based on a dataset of plant macrofossil assemblages from Nowak et al. (2019), which is a compilation of previously published and unpublished plant fossil collections. This dataset is used in the present paper to determine major biomes at the stage level, spanning from the Wuchiapingian to the Anisian. Each plant genus is associated with the major biome(s) it could potentially occur in, possibly exclusively (i.e., when a plant genus is known to be limited to habitats aligned with a certain biome; see Supplementary Table S2 in Nowak et al. (2020) for details).

The records are aggregated within a moving 100-km radius, where each genus is assigned a weight corresponding to the inverse of the number of major biome(s) it is associated with. For each aggregated assemblage, a vote for all possible biomes is then tallied based on the weights of all present genera. The final major biome is determined based on either i) the exclusive major biome if a characteristic representative is present, ii) the major biome whose combined weight is > 50 % in the area, or iii) the major biome manually assigned by the authors considering taxa, location, and neighborhood (if the combined weight for each calculated biome is < 50 % ). The resulting major biome distribution, thus, reflects the distribution of all plant macrofossil collections and corresponds to the area with comparable, climatically controlled plant and animal assemblages (Nowak et al., 2020). This method helps reduce the uncertainty in location and short-scale variability. However, some of the manual biome assignments are tenuous and mainly based on surrounding data points. In these cases, if a single calculated biome is indicated as the most likely one (as opposed to multiple candidates with similar combined weights, each being the case in approximately half of the manual assignments or a quarter of all aggregated assemblages), we use the latter in the present study. This approach aims to keep data points independent without introducing new interpretations.

The classification of major biomes used in this study comprises six categories, which are mostly adapted from those in Nowak et al. (2020), which, in turn, followed the set of biomes introduced by Ziegler (1990) as far as they could be applied to the fossil dataset at hand. The resulting simplified major biomes are briefly described as follows: the tropical everwet major biome includes various vegetation types developing under constantly hot and humid conditions, typically near the equator, but it can extend up to mid-latitudes. The tropical summerwet major biome represents intermediate vegetation between tropical everwet and desert, found in middle to low latitudes with marked seasonality and wet summers. The desert major biome includes both subtropical and mid-latitude deserts but is generalized to all latitudes in this study and is marked by water deficiency. The warm-to-cool temperate major biome encompasses vegetation ranging from evergreen (i.e., plants that keep their needles or leaves all year) to deciduous (i.e., plants that shed their leaves in autumn), affected by seasonal changes in climatic conditions. The cold temperate major biome is associated with areas with low evaporation, where the short growing season is mainly controlled by temperature and sunshine. The tundra major biome has a very short growing season and was not recorded by Nowak et al. (2020).

The simulated vegetation distribution is different between the three attractors but also varies along each stable branch, together with the atmospheric CO 2 content. The range of vegetation pattern on each attractor is determined by applying the asynchronous coupling procedure between the MITgcm and BIOME 4 models, as described in Supplementary Appendix SA at various positions along the branches. This includes a common value of 320 ppm and, in general, at the edges of each branch (see Figure 1B). For the warm state, which has a narrow stable branch, only two positions are selected, at 320 ppm and near the upper edge of the branch (328 ppm). The two models, the climate model MITgcm (Marshall et al., 1997a; Marshall et al., 1997b; Adcroft et al., 2004; Marshall et al., 2004) and the vegetation model BIOME 4 (Haxeltine and Prentice, 1996; Kaplan, 2001; Kaplan et al., 2003), have been described in the Methods section of Ragon et al. (2024). The convergence criteria between MITgcm and BIOME 4 are 1) the global SAT does not change between two iterations (within the uncertainty) and 2) the land surface fraction where albedo varies between two iterations is smaller than 10%. The corresponding values and surface imbalance over the ocean at each iteration of the coupling procedure are reported in Supplementary Table S6.

The biomes resulting from the simulations (28 biomes; see Kaplan, 2001) are grouped in the same major biomes described in Section 2.2 to facilitate the comparison with the geological records provided by Nowak et al. (2020) and eliminate perturbations coming from small differences between similar biomes. The corresponding classification is described in Table 1.

The present-day “grassland” is replaced by the herbaceous non-graminoid “forbland” because graminoids only appeared during the Cretaceous–Cenozoic (Gradstein and Kerp, 2012). There is consequently no direct equivalence between simulated grasslands and any fossil assemblage or fossil-based biome from the Permian and Triassic, but we can infer likely correspondences based on the climatic conditions they represent. Both tropical savanna and tropical forbland biomes have been included in the tropical summerwet major biome. Even if they slightly differ from the broadleaf forest category, their seasonality has been considered the main argument for their classification. The classification of barren, a desertic biome mainly associated with polar regions in present-day vegetation, is distinguished by the latitude of formation: the desert major biome is found below 40 ° latitude, while tundra is located above.

The criteria for temperate forbland biome in BIOME 4 correspond to desert-like conditions, more than everwet or seasonally wet. However, our simulations show that this biome is mostly formed in high latitudes, ≳ 40 ° (see for example gray areas in Figure 2). In polar regions, reduced evaporation gives rise to cold temperate vegetation, allowing it to thrive instead of forming deserts. For that reason, temperate forbland is classified as a cold temperate major biome (case 1 in Table 1). For completeness, we repeat the same analysis with an alternative interpretation, whereby temperate forbland is classified as a desert (case 2 in Table 1) in Section 3.4.

The paleontological record provides local information on vegetation, while the model simulates broader areas, posing a challenge for the direct comparison between the two. Visual comparisons (see example in Figure 2) can offer valuable insights into how well the simulation aligns with the records; however, this approach is not quantitative and does not always provide a clear answer.

Each plant fossil assemblage record is assigned to the model-grid cell corresponding to its location. Although the single outcrops are considered independent, their spatial resolution is sometimes higher than that of the model [i.e., records fall within the same grid cell of ∼ (300 km)²]. Because typical statistical fluctuations decrease with the square root of the number of observations, each distance is weighted based on the number of similar individual records in the same cell (see Equation 2). This results in the merging of plant fossil assemblages, representing the same biome within the same cell, thereby increasing their confidence. The merging contributes to a reduction in sample sizes n : from 97 to 56 (Wuchiapingian), 59 to 43 (Changhsingian), 66 to 46 (Induan), 88 to 42 (Olenekian), and 202 to 64 (Anisian). Note that the weight is applied only for the mean distance calculation, whereas for the median, only one record is considered among the identical ones within the same grid cell (see Section 2.5).

The similarity between the plant macrofossil record at a given stage and the simulated vegetation distribution of an attractor is estimated as follows. For each reported assemblage of fossil plants, we compute the smallest geodesic distance d to the nearest region where the same major biome has been predicted by the simulation, given by d = R ⋅ arccos sin ϕ M sin ϕ F + cos ϕ M cos ϕ F cos λ F − λ M , (1) where R = 6371 km is the Earth radius, ϕ is the latitude, λ is the longitude, and indices refer to the position of the fossil record ( F ) and the nearest point on the map where the considered major biome is simulated ( M ) . Figure 3 shows an example with a simulation and few paleobotanical records, where arrows correspond to the distance between the record and the nearest cell with the same major biome. We have checked that if the model predicts the same major biome at the position where it has been observed, then the distance is 0 (see the two records without arrow in Figure 3). When a major biome present in the fossil record does not appear in a simulation, the distance is arbitrarily set to the longest possible distance, i.e., the distance between the two poles, π ⋅ R ∼ 2 ⋅ 1 0 7 m.

We define the weighted mean distance d ̄ of the ensemble as follows: d ̄ = ∑ i = 1 n 1 n i   d i ∑ i = 1 n 1 n i , (2) where d i is given by Equation 1 for each observation i, n is the sample size, n i is the number of co-located identical records, and 1 n i is the weight applied to each point. The associated variance V is as follows: V = ∑ i = 1 n 1 n i   d i − d ̄ 2 ∑ i = 1 n 1 n i , (3)

The distances d i are calculated for each combination of fossil records (for each stage) and simulation (for each attractor and CO 2 value), together forming a sample. For a given stage, eight samples corresponding to the eight simulations are available. The one minimizing the mean distance d ̄ (or median distance) is assumed to better represent the vegetation distribution recorded at that stage. The difference between this sample and the others can either be due to i) a real difference between the populations from which the samples are drawn or ii) sampling fluctuations, meaning that geological records do not permit discrimination between the two simulations. In the first case, the simulation minimizing d ̄ is significantly better at reproducing the recorded vegetation, while in the second case, the data do not allow for the clear determination of which attractor minimizes the distance. To distinguish between these two cases, statistical tests are employed. These tests determine whether the difference between two samples, A and B , is significant.

The utilization of several methods (on means and medians), together with the classification of the temperate forbland biome into two different major biomes (cases 1 and 2 in Table 1), allows us to test the robustness of the results (see also Supplementary Appendix SB).

The simulated distribution of the major biomes corresponding to the attractors at various atmospheric CO 2 contents is shown in Figure 4. Tundra, barren, and land ice biomes, which have no correspondence in the macrofossil assemblages, are mostly present in polar regions in the cold and warm states, especially at the northern polar latitudes coinciding with the sea ice extent (Ragon et al., 2024).

The simulated major biomes in the cold and warm states (Figures 4D–H) are both dominated by desert at tropical and subtropical latitudes. Along 60 ° latitude, a band of warm-to-cool temperate vegetation is present, whereas in polar regions, cold temperate plants are favored and switch into tundra when moving poleward. The space occupied by the tundra in the northern polar region reduces in favor of cold temperate as the atmospheric CO 2 content increases along the stable branch. The changes in atmospheric CO 2 content mainly affect the biomes in the high latitudes of the northern hemisphere in both attractors. In the southern polar region, the vegetation is largely dominated by temperate forbland (i.e., either cold temperate or desert major biomes, as discussed in Section 2.3).

In the hot state, the warm-to-cool temperate major biome is also present but shifted poleward compared to the other attractors, so it dominates in the northern polar region. The southern polar region is occupied by both the warm-to-cool temperate major biome and temperate forbland, with the extent of the latter reducing as the atmospheric CO 2 content increases (see Figures 4A–C). At tropical and subtropical latitudes, desert areas are generally present, and plants adapted to tropical summerwet and everwet environments are well-distributed along the tropics and 50 ° latitude, respectively. Their extent increases together with the atmospheric CO 2 content along the hot-state branch. In general, the increase in atmospheric CO 2 content favors warmer and wetter major biomes: tropical everwet/summerwet over desert in mid-latitudes/tropical regions and warm-to-cool temperate over temperate forbland in the southern polar region.

Figure 5 shows, at the five considered stages, the mean distance d ̄ and associated standard error V (green circles and bars, respectively) for each simulation, along with the p -value. See also Supplementary Table S1.

For the two oldest stages, the Wuchiapingian and the Changhsingian (Lopingian), the vegetation distribution modeled for the cold state at 320 ppm provides the best match with the fossil records. However, it is not significantly different from the vegetation distribution modeled at other positions on the cold branch or from the warm state. In contrast, the vegetation simulated in the hot state is significantly different, and thus, it is not a good candidate to reproduce the records of these stages.

In the case of the Induan, the hot state at 308 ppm minimizes d ̄ but is not statistically different from any of the other simulated conditions, except for the upper branch of the hot state, which can be excluded. For both the Olenekian and Anisian, the lower part of the hot-state branch, i.e., at 308 ppm and 320 ppm, significantly differs from all the other conditions and better matches the records.

These results suggest a transition from the cold or warm state in the Lopingian to the hot state in the Olenekian, in concordance with the global warming recorded from the Lopingian to the Early Triassic (Retallack, 1999; Joachimski et al., 2012). The hot state persists during the Anisian, thus marking the stabilization of the climate system, as recorded, for example, in δ 13 C isotopic records (Payne et al., 2004). The unclear signal during the Induan might be the result of oscillations occurring at smaller temporal scales, as observed in the proxy of temperature δ 18 O (Romano et al., 2013).

The test on medians shows a trend comparable to what is observed for that on means, as shown in Figure 6; Supplementary Table S2. For both the Wuchiapingian and the Changhsingian, the cold state, not statistically different along the branch or from the warm state, has a smaller median and thus matches better with the data, while the hot state has a significantly larger median and thus can be excluded.

For the Induan, the median in the warm-state simulation with 328 ppm is lower than that for the other simulations and is statistically different from the whole hot state and the lower edge of the cold-state branch.

For the Olenekian and the Anisian, the hot state at 350 ppm minimizes the median but is not significantly different from the same attractor with other atmospheric content of CO 2 . For the Olenekian, the difference with the warm state at 328 ppm is also non-significant. However, in both cases, the hot state differs from the cold state, which can thus be excluded.

The classification of temperate forbland into cold temperate is questionable, as discussed in Section 2.3. In this study, we test the robustness of the results against the classification of this biome by including it in the desert major biome (case 2) instead of cold temperate (case 1).