East palearctic treefrog past and present habitat suitability using ecological niche models

Focal species

There are currently four recognized species of Dryophytes in the east palearctic: D. japonicus, D. immaculatus, D. flaviventris and D. suweonensis. Dryophytes japonicus is split into two clades (hereafter Clades A and B) with the Clade A being found in Japan northeast of the Chugoku-Kansai border (Dufresnes et al., 2016) and Clade B, likely associated with the name Hyla stepheni (Dufresnes et al., 2016; Borzée et al., 2020b), is found widely across mainland northeast Asia and Japan southwest of the Chugoku-Kansai border (Dufresnes et al., 2016). Dryophytes immaculatus, D. flaviventris, and D. suweonensis are combined into the Dryophytes immaculatus group, with D. immaculatus found in mainland China, D. suweonensis found in the Republic of Korea and the Democratic People’s Republic of Korea, north of the Chilgap mountains, and D. flaviventris found in the Republic of Korea south of the Chilgap mountains and north of the Mangyeong river (Borzée et al., 2020b). The species within the D. immaculatus clade are generally restricted to the low elevation alluvial plains around the Yellow sea (Borzée et al., 2020b) while D. japonicus is found at variable elevation, breeding in a variety of water bodies (Tsuji et al., 2011; Roh, Borzée & Jang, 2014; Borzée et al., 2019b), and human alteration of the landscape resulted in hybridization between the two groups (Borzée et al., 2020a).

Climatic suitability modelling and projection to past climates

We used maximum entropy (MaxEnt ver. 3.4.4; Elith et al., 2011; Phillips et al., 2017) modelling to predict present and past distributions of Dryophytes in the east palearctic, separating the Dryophytes japonicus group into two clades (called A and B in Dufresnes et al., 2016) and the Dryophytes immaculatus group into three species (D. flaviventris, D. immaculatus and D. suweonensis; Borzée et al., 2020b). Occurrences were obtained from the Global Biodiversity Information Facility (GBIF.org, 2021) and non-focal surveys and opportunistic observations conducted by the authors (Supplemental File). In filtering GBIF data, we removed occurrences with more than 4,000 m accuracy, which was the cell size of the environmental layers used to train our models. We omitted occurrences by the National Institute of Ecology Korea (NIEK) since the methodology used by the institute (gridded sampling along transects) differed from other data sources and led to uneven sampling which overestimated common species (D. japonicus) while underestimating rare species (D. flaviventris and D. suweonensis). To further reduce sampling bias, we reduced the occurrences of D. japonicus in the Republic of Korea to 20% using the “Subset” tool in ArcMap 10.8.1 (Esri, Redlands, CA, USA). We finally removed occurrence points of D. japonicus along the Chugoku-Kansai border since the lineage could ambiguously belong to either clade A or B of the species.

When running MaxEnt, we used seven environmental layers (Table 1) with a 2.5 arc-minute (0.0417 decimal degree or ~4 km) resolution. This resolution was used as it was the smallest scale available for the Last Glacial Maximum downscaled climate reconstruction. We used seven independent bioclimatic variables (WorldClim 1.4; Hijmans et al., 2005) with low multicollinearity (Pearson’s r < 0.8) which were previously determined to influence present and past distributions of D. japonicus in the east palearctic (Dufresnes et al., 2016). We opted to only use climatic variables and exclude terrain variables that potentially correlate to present-day distributions (Borzée et al., 2020b) as current topography (e.g., elevation) may not correlate to the same climatic conditions across geologic time scales (Peterson et al., 2011; Jarnevich et al., 2015; Bobrowski et al., 2021). For example, other species have been predicted to shift their elevational distribution as a result of past climate change (Yousefi et al., 2015). For model replication, we used the cross-validate run type option with five replicates (5 k-fold), which uses 80% of occurrences for training and 20% for testing. We selected the option to remove duplicate presence points within the same cell of the environmental layers. To address overfitting (Radosavljevic & Anderson, 2014), we created models at three levels of background points (5,000, 10,000, 15,000) and four levels of randomization multipliers (0.5, 1.0, 1.5, 2.0). We finally projected our trained models to past climates of the last glacial maximum (LGM; 22 kya; WorldClim 1.4; Hijmans et al., 2005) and last interglacial (LIG; 130 kya; WorldClim 1.4; Otto-Bliesner et al., 2006). For our climate projected models, we ran multivariate environmental similarity surface (MESS) analysis (Elith, Kearney & Phillips, 2010) to identify areas of novel past climates with no modern analog where a projected climate suitability model may require more careful interpretation. For each model, we calculated the area under the curve (AUC) for training and test data, true skill statistic (TSS; Allouche, Tsoar & Kadmon, 2006; Lobo, Jiménez-Valverde & Real, 2008), and percent omission of test data for minimum and 10 percent omission thresholds. TSS was calculated from MaxEnt output test and background model values by (from Allouche, Tsoar & Kadmon, 2006):

$$Sensitivity={\displaystyle \frac{a}{a+c}}$$

$$Specificity={\displaystyle \frac{d}{b+d}}$$

$$TSS=Sensitivity+Specificity-1$$where:

a = number of occurrences (test points) for which presence was correctly predicted by the model.

b = number of occurrences (background points) for which the species was not found but the model predicted presence.

c = number of occurrences (test points) for which the species was found but the model predicted absence.

d = number of occurrences (background points) for which absence was correctly predicted by the model.

We additionally calculated the difference between training and test AUC for each model. We selected best models for each clade based on these fit statistics, choosing models with low difference between training and test AUC, minimum test omission closest to 0% and 10 percent test omission closest to 10% (or below). To calculate the area of suitability for each clade, we thresholded the models at the maximum test sensitivity plus specificity threshold (true scale statistic or TSS; Allouche, Tsoar & Kadmon, 2006), a threshold which maximizes true presences and absences while minimizing false presences and absences. Finally, to visualize areas of potential climatic refugia, we used the “Cell Statistics” tool (specifying “sum” as the overlay statistics) in ArcMap 10.8.1 (ESRI, Redlands, CA, USA) to create stability models for each clade by calculating the sum of cloglog outcomes of our suitability models across time periods (Loera, Ickert-Bond & Sosa, 2017).

Model evaluation

Models for all clades produced AUC values greater than 0.9 (Table 2), which is considered excellent fit (Dolgener et al., 2014). Models for the D. japonicus Clade B (TSS = 0.8077), D. immaculatus (TSS = 0.7946) and D. flaviventris (TSS = 0.8329) had very good presence-absence prediction capability while models for the D. japonicus Clade A (TSS = 0.9475) and D. suweonensis (TSS = 0.9127) had almost perfect presence-absence prediction capability (Table 2). See Figs. S1–S5 for all tested models and statistics.

Present and past distributions and areas

The D. japonicus Clade A had a similar suitable range (defined here as the thresholded area of suitability) during the last interglacial (LIG) to its current suitable range (Fig. 1). During the last glacial maximum (LGM), its suitable range shifted south and contracted in the northern part of its current suitable range, with a disjunct suitable range in the Yellow Sea Basin and mainland China. The D. japonicus Clade B also shifted south during the LGM with a LIG suitable range similar to its current suitable range in mainland northeast Asia and southern Japan (Fig. 2). Both D. japonicus clades experienced expansion from the LIG to the LGM, and then subsequent detraction in the present (Table 3). Dryophytes immaculatus experienced suitable range contraction from the LIG to the LGM, but an expansion from the LGM to the present (Fig. 3). Dryophytes flaviventris experienced a suitable range detraction across all three time periods, becoming the clade with the smallest suitable range in the present at 10,760 km². The clade’s suitable range in the LIG covered large patches in China and the Korean peninsula and had a large area of suitability in the Yellow Sea Basin during the LGM (Fig. 4). Dryophytes suweonensis experienced a marked increase in suitable range from the LIG to the LGM (29,480 km² to 816,268 km²) mostly in the Yellow Sea Basin (Fig. 5), but in the present its suitable range is similar to its LIG suitable range.

For D. japonicus Clade A, the highest climatic stability areas were in lowland and coastal areas of central Japan (Fig. 1), while the highest climatic stability for D. japonicus Clade B occurred in the Korean peninsula and southern Japan (Fig. 2). For D. immaculatus, high climatic stability was found along the Yangtze River (Fig. 3). High climatic stability for D. flaviventris was limited to coastal lowlands in the southwestern corner of South Korea (Fig. 4). Climatic stability for D. suweonensis occurred in the lowlands of the western Korean peninsula, with a small area of stability along the eastern coast of the Democratic People’s Republic of Korea around the city of Wonsan (Fig. 5).

MESS analysis was similar for all clades, with a small area of novel climatic conditions during the LGM in the north and a larger area during the LIG mostly in present-day China (Figs. 1–5, S6). During the LGM, the main novel limiting climatic variables were Bio1, Bio2 and Bio3. During the LIG, the main novel limiting climatic variables for all clades were Bio2, Bio5 and Bio15.

Variable importance

The variables contributing the most to the distribution of the D. japonicus Clade A was Bio19 (precipitation of coldest quarter) at 49.69% (Table 4), while Bio15 (precipitation seasonality) had the highest permutation importance at 84.36% (Table 5). For the D. japonicus Clade B, elevation was the highest for both percent contribution and permutation importance. Elevation was also the highest variable for D. flaviventris, with a contribution of 68.19% and permutation importance of 67.30%. The variable with highest contribution (37.36%) and permutation importance (48.88%) for D. immaculatus was Bio1 (annual mean temperature). For D. suweonensis, elevation had the highest percent contribution at 36.56%, while Bio12 (annual precipitation) had the highest permutation at 76.68%.

Habitat suitability and stability

The highest stability areas for the two D. japonicus clades match their current respective ranges. For the D. japonicus Clade A, this range is in central Japan. The highest stability areas for the D. japonicus Clade B are in the Korean peninsula, between Vladivostok and Khabarovsk along the Ussuri River, and in southern Japan. Comparison with previous research (Dufresnes et al., 2016) suggests a possible correlation between high stability areas and a greater number of haplotypes of Dryophytes japonicus clades in the Korean peninsula and southern Japan. Future research may find a significant correlation between climatic stability and genetic diversity for the two clades. There is also moderately stable area in mainland China west of the Yellow Sea, which makes it surprising that the clade is not found in this region. This may simply be because the species did not disperse to the area while it was accessible, or competition with D. immaculatus or another species (Borzée et al., 2016; Borzée, Kim & Jang, 2016) may have prevented expansion.

High stability areas for D. immaculatus follow the Yellow River Basin where the species is most prevalent. Dryophytes flaviventris has high stability areas in its current range as well as an area south of its range (South Jeolla province of Rep. of Korea). As the species is declining and facing continued habitat destruction, translocation to the South Jeolla province of Rep. of Korea may aid conservation efforts. Additionally, since there is high stability in two coastal areas of Shandong and Hebei, molecular tools or call analysis may be needed to confirm the species in those areas. The highest stability areas for D. suweonensis also match with its current known range with the exception of a high stability area along the east coast of the Democratic People’s Republic of Korea. However, as the species has not been recorded in the area and is geographically isolated from the known range, it is unlikely the species occurs there.

During the LGM, all clades had suitable area in the then-exposed Yellow Sea Basin. This would have allowed for contact and created a situation where admixture was possible. As previous research has estimated the divergence time between D. immaculatus and the D. suweonensis group at around 1.02 mya (Borzée et al., 2020b), it is possible that a minor glaciation during this time (Ehlers, Gibbard & Hughes, 2018) saw the clades splitting from a shared range in the Yellow Sea Basin. However, the merging of all the rivers into the Yellow Sea Basin may have created a barrier that was too difficult for Hylids to cross as observed by the Vistula River segregation of Hyla arborea and Hyla orientalis in Poland (Stöck et al., 2012) and in the Middle East where the Dead Sea Rift has divided Hyla savignyi and Hyla felixarabica (Dufresnes et al., 2019).

Range expansions for the D. japonicus clades, D. flaviventris and D. suweonensis during the LGM are in line with previous studies which have seen LGM expansion in a temperate hylid Hyla sarda (Bisconti et al., 2011). Cooler temperatures that would allow temperate species expansion might have also limited the LGM range of D. immaculatus which is currently found at lower latitudes and more subtropical climates.

Variable responses and species ecology

Variables having the highest percent contribution and/or permutation importance in ecological models indicate potential preference or limiting factors for species. For the D. japonicus Clade A, these variables included precipitation of coldest quarter and precipitation seasonality. Since the region is generally marked by dry winters, this shows a preference towards wetter winters and less precipitation variability between seasons (Fig. 6). This is interesting because a similar preference is not seen in Clade B. Instead, elevation is the variable with the highest contribution and importance, with the species found at low elevations mostly between 0 and 400 m (Fig. 6). Elevation was similarly important for D. flaviventris and D. suweonensis, whose cloglog responses both peak around sea level (0 m). This fits with the species’ habitat preferences for alluvial plains and makes the species particularly vulnerable to sea level rise. For D. suweonensis, the high permutation importance of annual precipitation indicates a preference for moderate annual rainfall with a tolerance for higher rainfall better than lower rainfall and peaking at 1,240 mm. Finally, for D. immaculatus, there is a preference for higher annual mean temperature, which relates largely to its range as the clade present at the lowest latitudes.

Ecological niche modeling also allows for among-species comparisons in responses to environmental variables. Interestingly, annual mean temperature for all clades peaked between 10 to 16 °C indicating preference for warmer climates, with the D. japonicus Clade B tolerating cooler temperatures likely in its northern range in Russia. Maximum temperature of warmest month peaks below 35 °C for all clades likely because of critical thermal limits. A study of another hylid, Dryophytes versicolor, showed a decrease in tadpole speed past about 32 °C with a sharp decrease occurring around 34 °C (Katzenberger et al., 2014).