Ttyn Y. Petrenko, Kirill A. Korznikov,＊, Dmitry E. Kislov, Ndezhd G. Belyev,Pvel V. Krestov

a Botanical Garden-Institute FEB RAS, Makovskogo St. 142, Vladivostok, 690024, Russia

b Institute of Geography RAS, Staromonetniy 24-4, Moscow, 119017, Russia

Keywords:Climate change Bioclimatic modeling Species distribution model Korean pine Northeast Asia

ABSTRACT Background:Pinus koraiensis Siebold&Zucc.(Korean pine)is a key species of the mixed cold temperate forests of Northeast Asia.Current climate change can significantly worsen the quality of P.koraiensis habitats and therefore lead to a large-scale structural and functional transformation of the East Asian mixed forests. We built a species distribution model (SDM) for P. koraiensis using the random forest classifier – a versatile machine learning algorithm, to discover overlap areas of potential species occurrence in the climate condition of the Last Glacial Maximum (～21,000 year before present) and in the projected future climates (2070 year), from which possible permanent refugia for P. koraiensis were identified.Results: Using the random forest supervised learning algorithm, we developed models of the modern distribution of P.koraiensis in accordance with the five selected bioclimatic variables(Kira’s warmth and coldness indices,the index of continentality, the rain precipitation index, and the snow precipitation index). In addition to current climatic conditions, we performed this analysis for the climate of the Last Glacial Maximum and for the future projected climate(2070)under scenarios RCP2.6 and RCP8.5.Among the predictors,the rain index appears to be the most significant.The land area estimates with high suitability for P.koraiensis was 303,785 km2 under current climatic conditions, 586,499 km2 for the Last Glacial Maximum, and 337,573 km2 for the future (2070) period under the RCP2.6 scenario, and 397,764 km2 under the RCP8.5 scenario.Conclusions: Most of the potential range of P.koraiensis during the Last Glacial Maximum was located outside the current distribution area of the species.The climatically suitable P.koraiensis habitats will likely disappear in the western part of its modern range. In the southern part of the range, which includes glacial refugia, the areas of continuous distribution of the P. koraiensis populations since the end of the Pleistocene are expected to be fragmented, but some localities in the north of the Korean Peninsula, northeast China, southern Primorye (Russia),and central Honshu (Japan) with suitable climatic conditions for the species will support the existence of populations.

1. Background

Pinus koraiensis Siebold & Zucc., Pinaceae (Korean pine) is a key species of mixed deciduous broadleaved-coniferous forests (DBCFs) in the cold-temperate zone of Northeast Asia (NA) (Kolesnikov, 1956;Nakamura and Krestov, 2005; Krestov et al., 2006; Nakamura et al.,2007). The DBCFs with the dominance or presence of P. koraiensis formerly covered a vast territory in the region.Currently,approximately 70%of this forest type has been severely destroyed or transformed,and intact ecosystems remain only fragmentary,mostly in the territory of the Russian Far East (Dom‵enech et al., 2019). Some physio-ecological features of P.koraiensis,such as the currently observed trends of decreasing radial increments in populations located near the southwestern border of the species’ range, are associated with changes in climatic conditions(Lyu et al., 2017). Against the background of a changing climate, dramatic human-induced transformation of the DBCFs, including P.koraiensis,forms a generally unfavorable environment for the sustainable existence of the populations of their consorts.

Species distribution modeling (SDM) based on machine learning algorithms in the last two decades has firmly strengthened the arsenal of methods for analyzing ecological niches (Guisan and Zimmermann,2000; Elith et al., 2006). SDM at regional and global levels involves bioclimatic indices as predictors of habitat suitability for a particular species since the climate controls the distribution of species on a wide spatial scale. SDMs are considered to represent and accurately reflect complex climatic conditions when the potential distribution of the species under evaluation is likely (Araújo and Guisan,2006).

Global climate change leads to a shift in the boundaries of species ranges (Chen et al., 2011) and changes the structure and functioning of ecosystems(Vill′en-Per′ez et al.,2020). The use of climatic niche models allows prediction of the change in climatic conditions suitable for a particular taxon in accordance with various scenarios of the predicted climate (Pearson and Dawson, 2003; Hijmans and Graham, 2006). This opens up the possibility of adaptive management in areas such as forestry and agriculture as well as ex situ species conservation (Janowiak et al.,2014;Schelhaas et al.,2015).The use of retrospective models of climatic niches often supports research on molecular phylogeny and genosystematics and is a working tool for paleobiology(Sakaguchi et al.,2010).

Modern SDMs, in their essence, are statistical models fitted to observed data. These include, but are not limited to generalized linear models (GLMs),generalized additive models (GAMs),ensembles of randomized decision trees,as well as such state-of-the-art generalizations as boosting and stacking approaches(Beery et al.,2021).Being equivalent to logistic regression – a widely used algorithm in various SDMs, a single-layer perceptron is another example of using neural networks in the SDM field(Hegel et al.,2010).The latest advances in neural networks of deep architectures in various applied problems of human beings inspire extending of deep learning algorithms to SDMs.It was shown that convolutional neural networks of deep architecture in the case of multi-species distribution modeling can overcome classical approaches,such as MaxEnt and logistic regression model (Hegel et al., 2010).Another group of versatile machine learning methods is presented by the support vector machine (SVM) concept and its generalizations, such as kernel SVM(Evgeniou et al.,2001).These algorithms can build decision regions with a predefined level of confidence, excluding data points belonging to non-confident areas. The kernel trick allows handling the SVM approach of non-linear decision boundaries (Suykens, 2001). The main drawback of SVM approach is its training complexity, which requires a considerable amount of training time in cases when the dataset is presented by millions of data points(Mor′e,1978).The comparison study presented in Kaky et al.(2020)revealed that random forest(RF)classifier and MaxEnt can show in certain cases comparable performance metrics,whereas other algorithms, such as GLMs and SVM, have relatively poor accuracy. Because of simplicity and the ability to handle complicated decision regions in factor spaces,we chose the RF algorithm as a primary tool for developing our model.

This article focuses on modeling the potential distribution of P. koraiensis using data on the growth of the species in its natural habitats(not in the conditions of arboriculture or plantations)and the WorldClim 1.4 climate dataset (Hijmans et al., 2005). With the aid of the past climate data MIROC-ESM (Watanabe et al., 2011; Kawamiya et al., 2020), we reconstructed the spatial distribution of areas with climatic conditions suitable for P. koraiensis in the periods of the Last Glacial Maximum(LGM,～21,000 years before present).We also built prognostic models of the potential distribution of P. koraiensis for 2070 year according to the climate change scenarios RCP2.6(representative concentration pathway)(van Vuuren et al., 2011) and RCP8.5 (Riahi et al., 2011). Both prognostic and retrospective SDMs were developed using ensembles of decision trees(RF algorithm). The decision tree is an efficient algorithm for handling complicated decision boundaries in multi-factor spaces(Rokach and Maimon, 2008). Another advantage of decision trees and their ensembles for SDMs is the ability to estimate climatic variable importance and impacts on the observed species distribution. In our study, we focused on building SDM for P. koraiensis using the RF classifier – a versatile machine learning algorithm, to overlap ranges of potential species occurrence in the climate condition of the LGM and in the projected future climates,which could be considered permanent refugia for P.koraiensis.

2. Materials and methods

2.1. Range and natural condition of P. koraiensis

P.koraiensis occupies the south of the mainland Far East of Russia,the northeastern part of China,the mountainous regions of the northern part of Korean Peninsula and the central part of Honshu Island of Japan(Sochava,1967;Okitsu and Momohara,1997;Krestov,2003;Nakamura and Krestov,2005;Aizawa et al.,2012;Su et al.,2020).Artificial forest stands of P. koraiensis exist in northern Japan, northeastern China, the southern part of the Korean Peninsula and the southern part of Sakhalin Island of Russia.The mixed forests with P.koraiensis occur in areas with complex, mainly mountainous terrain. The climate of the region is influenced by the East Asian monsoon; winters are cold and nearly snowless, and summers are relatively cool and humid (Krestov, 2003).The distribution area of P.koraiensis corresponds to the zone of northern temperate mixed forests (Nakamura and Krestov, 2005; Krestov et al.,2006).

2.2. P. koraiensis occurrence data

The geographical coordinates of P. koraiensis localities (presence points)were extracted from the herbarium labels from the collections of the V.L.Komarov Botanical Institute of the Russian Academy of Sciences(LE), the Botanical Garden-Institute of the Far Eastern Branch of the Russian Academy of Sciences (VBGI), Moscow State University (MW),and the Federal Research Center for Biodiversity (VLA). Some distribution data were obtained from the Global Biodiversity Information Facility(GBIF) (https://www.gbif.org/), as well as from our own observations.We used only those herbarium samples, the coordinates of which could be identified with an accuracy of at least 2 km.

The GBIF presence points were subjected to additional filtering. We excluded from the analysis all localities for which the location coordinates were rounded to one hundredth of a degree because they were judged less accurate. We assumed that coordinate data with decimal representation to only the hundredths place would not offer the geographic level of detail necessary for use in our analysis. Data with coordinates outside the natural distribution area (forest plantations on the islands of Hokkaido and Sakhalin)or urbanized territories(gardens,parks)were also excluded.

This resulting dataset was then searched with an algorithm to remove the presence points located closer to each other than at a distance of 2 km(essentially duplicate points).The algorithm was implemented using the Python GeoPy, package https://geopy.readthedocs.io/en/stable/. After the removal of these points,we calculated the Average Nearest Neighbor,index for the remaining data(implemented in the ArcMap 10.8 program)https://desktop.arcgis.com/en/arcmap/latest/tools/spatial-statisti cs-toolbox/average-nearest-neighbor.htm,which compares the observed average distance between all points with species presence with the expected distance for a set of evenly distributed points.If the index is less than 1,the pattern exhibits clustering;if the index is greater than 1,the trend is toward dispersion or competition. Thus, we managed to avoid pronounced effects of data imbalance when one region could be represented by a disproportionately large number of points of presence. As a result,207 unique presence points of the species were used in modeling(Fig. 1a; Table S1). The number of pseudo-absence points randomly placed throughout the simulation area was estimated to be 10 times larger than the number of presence points.

Fig. 1. Potential range of Pinus koraiensis according to the SDM; (a)current climate condition, black dots indicate species presence points used in modeling; (b)Last Glacial Maximum (～21,000 years BP) climate condition; (c) 2070 year RCP2.6 climate condition; and (d) 2070 year RCP8.5 climate condition.

2.3. Climate data

To model the distribution of P. koraiensis, we selected the five most informative bioclimatic indices from a set of more than 30 indices provided by Nakamura et al. (2007) and Noce et al. (2020). The highest contributions to the model were shown by Kira’s warmth index(WKI,the sum of the average monthly temperatures above+5°C);Kira’s coldness index (CKI, absolute value of the sum of the average monthly temperatures below +5°C); the index of continentality (IC, the difference between the year maximum and minimum average monthly temperatures);the rain precipitation index (Pp, the amount of precipitation for the period with positive average monthly temperatures); and the snow precipitation index (Pn, the amount of precipitation for the period with negative average monthly temperatures) (Kira, 1977; Nakamura et al.,2007).

The index values were calculated from the data on average monthly temperatures and precipitation totals in WorldClime v.1.4 (Hijmans et al.,2005)with a spatial resolution of 30 arcseconds(～0.0083°),which were extracted from the source data files using the GDAL, https://gdal.org/. Similar data presented in the framework of the MIROC-ESM climate model (Watanabe et al., 2011; Kawamiya et al., 2020) were used to reconstruct the climatic situation during the LGM and forecast the climatic situation for the year 2070. The prognostic modeling was performed in accordance with two scenarios of global climate change: a restrained(optimistic) RCP2.6(implies an increase in the average planetary temperature by 0.3°C–1.7°C by 2100) (van Vuuren et al., 2011)and a pessimistic RCP8.5(increase by 2.6°C–4.8°C)(Riahi et al.,2011).

The test for multicollinearity of the selected bioclimatic indices was performed by checking the completeness of the rank of the data matrix using the NumPy,package https://numpy.org/and by means of omcdiag function from the mctest R-package (Imdadullah et al., 2016). It was found that determinant of the covariance matrix for selected subset of climatic variables significantly differs from zero(at least,at significance level 0.05).

Feature selection was implemented according to a recursive feature elimination algorithm provided by Scikit-learn package.Having applied it to the original subset of climatic variables multiple times,we ended up with slightly different subsets of optimal climatic variables (including widely used BIO1-BIO19 climatic indices)each time.Finally,we decided to use predefined, but sub-optimal composition of climatic variables characterized by a trade-off between clearness of interpretation and resulting model performance.Therefore,we concluded that it is better to choose a few predictor variables having clear biological interpretation rather than select slightly optimal subset of variables, which have implicit and unclear effects on species distribution(Santini et al.,2021).

2.4. Modeling and statistical processing

The identification of nonlinear relationships between the distribution of P. koraiensis and bioclimatic parameters was performed using the RF machine learning method implemented in the Python programming language in the Scikit-learn package(Pedregosa et al.,2011).Scikit-learn is a general purpose machine learning package aimed at fast prototyping,validating and deploying supervised and unsupervised learning models.It allows to formulate the process of building SDM at a high level of abstraction in Python programming language and efficiently performs all necessary steps of developing supervised learning models, including feature engineering and feature selection stages, cross validation and model testing. Additional parameters for the RandomForestClassifier during the grid search cross-validation procedure were found to be equal to their values used in similar models (Korznikov et al., 2019). In particular,the optimal number of random trees was found to be equal to 100, and the maximum tree depth was limited to 10. The remaining parameters of the RF corresponded to their default values accepted in Scikit-Learn.

The built model was evaluated using the continuous Boyce index(CBI) (Boyce et al., 2002), which is calculated using only the points of presence of the species,based on 100 iterations by randomly dividing the original set of spatial data into a training set(75%of the points)and a test set(25%of the points).The use of CBI to assess the quality of the model is preferable to using ROC AUC,since it is based solely on empirical data on the placement of localities of the species, without referring to pseudo-presence points (Lobo et al., 2008). The calculations were performed using the SciPy, (stats module) and the NumPy, https://www.s cipy.org/ of the Python environment. SciPy and the older, related NumPy are general purpose packages for efficient manipulating arrays and performing fast scientific computations in Python-based computational environments.It includes a large number of algorithms and special functions that are not limited to the machine learning field.

The contribution to the final model of each of the five predictors was evaluated by the feature_importances_attribute implemented for the RandomForestClassifier from the Scikit-learn package.

Climatic variables are not the only ones influencing species distribution. Human activities have a profound disturbance on nature and impact the final accuracy of species distribution models(Sanderson et al.,2002). One of the ways to handle this is to introduce the Human Influence Index and include it as a predictor variable.This index is composed of eight factors (human population density, railroads, roads, navigable rivers, coastlines, nighttime lights, urban footprint, land cover). However, our training data points fall almost entirely to the territory of the Russian Far East, which is characterized by relatively small human population density. Since values of Human Influence Index at training data points are expected to be almost constant,we decided to exclude this index from our model.

The result of using the trained classifier is a probability map(from 0–the presence is improbable,to 1–the maximum probability of presence)of habitat suitability for P. koraiensis. For practical purposes, such as calculating the area of the territory that a species can potentially occupy,we represented the probability maps in binary form,namely,“the species is absent” or “the species is present”. The transition to binary representation requires finding the optimal value of the threshold, the excess of which by the probability of presence will be identified with the “presence”of the species at a given point.To estimate the threshold value of the probability of the obtained probability map, we considered the problem of optimizing the average value of the maxSSS indicator (Liu et al.,2013),calculated on the basis of random partitions of the original spatial data set into 100 training and test set pairs for evaluation. The similar optimization problem was stated against true skill statistic and F1 score metric. To validate obtained optimal value of the threshold, we followed an expert-based approach (Liu et al., 2015; Konowalik and Nosol, 2021). As a result of computational experiments, we concluded that optimal values of maxSSS yields to distribution maps of P.koraiensis which are consistent with an expert-based assessments.The analysis was conducted using Scikit-learn packages.

In summary,the process of building the SDM using RF classifier was consisted in the following stages: 1) collecting P. koraiensis occurrences data; 2) data preprocessing (removing duplicates, local point density balancing,generation of pseudo-absence points);3)expert-based feature selection stage; 4) grid search for the best set of model parameters (the number of trees,tree depth,available criterion for trees constructing);5)finding the best threshold value (by maximizing maxSSS and expertbased approach); 6) applying the model to the past, current, future climatic data and results interpreting.

Finally, we calculated response curves for each bioclimatic variable.Response curves in their essence are smooth estimations of modeled probability of species occurrence for a fixed value of a particular bioclimatic variable. Its interpretation is straightforward: higher values on response curves correspond to a higher probability of species occurrence and suitability of climatic conditions.

3. Results

The resulting model has a high predictive accuracy:the value of CBI is 0.887±0.096(the average value and the error of the average,n ＝100).Among the predictors, the rain index (Pp) appears to be the most significant.The contributions of the remaining four parameters were nearly equal(Table 1).Response curves,as a function of the dependence of the probability of the presence of a species on the values of climate indices,are shown in Fig. 2. The total area of a climatically suitable territory in our model is almost 500,000 km2less than the area estimated in the model of P. koraiensis distribution by Shitara et al. (2021) over the maximum sensitivity plus specificity logistic threshold value, although 100,000 km2greater than the area under high suitability values(0.5–1)in the cited model.

Within the framework of the obtained model,the threshold value of the probability of presence,found from the maxSSS value, was 0.263±0.171 (n ＝100). We consider the territory with a predicted probability higher than this value to be highly suitable for P.koraiensis because of the binary prediction estimation. This area, in accordance with the current climatic situation,is 303,785 km2(Fig.1a).The area of such a territory for the considered part of East Asia during the LGM is estimated to be approximately 586,499 km2(Fig. 1b). Most of this territory is currently located on the shelf of the Yellow Sea and the Sea of Japan.The area of the territory with highly suitable climatic conditions by 2070 within the framework of the RCP2.6 scenario will be 337,573 km2, and if the RCP8.5 scenario is implemented,it will be 397,764 km2.However,such a suitable area increase is in accordance with the binary predictor assessment-based maxSSS threshold value.The area of the territory with higher levels of suitability (>0.4) will be catastrophically reduced throughout NA(Fig.1c and d).

The identification of areas where P.koraiensis has persisted over time indicates the existence of long-term stable refugia in such areas (Tang et al.,2018).These areas should deserve the highest priority in the field of preserving the gene pool of this species,which,in these areas,should be represented by the most ancient populations. In order to visualize the geographical locations where P. koraiensis has always had favorable conditions,the obtained species distribution models of the species in the past(LGM), the present and future (2070) were superimposed on each other.The intersection of potential areas occupied by this species in different time sections showed the position of long-term stable refugia(Fig.3).

4. Discussion

4.1. Retrospective modeling: the LGM

The climatic conditions of the LGM maximum in the Far East were characterized by a colder and much drier climate compared to themodern climate(Ju et al.,2007).Due to the global decrease in sea level,a significant part of the modern marine shelves were the terrestrial areas represented by coastal, gently sloping plains (Clark and Mix, 2002).During the LGM,a significant part of the range of P.koraiensis was located outside the area of its modern distribution.The simulation results show that certain areas of the modern shelf of the Yellow Sea and the Sea of Japan were suitable for P.koraiensis.Similar results were obtained for P.koraiensis (Shitara et al., 2021) and other temperate trees species in NA(Sakaguchi et al.,2012;Kimura et al.,2014).

Table 1 Feature importance for the Pinus koraiensis SDM.

Fig.2. Response curves for probability of presence of Pinus koraiensis;mean value(solid)and 95%envelope,n ＝100.Abbreviations:WKI,Kira's warmth index;CKI,Kira's coldness index; IC, continentality index; Pp, rain precipitation index; and Pn, snow precipitation index.

The results of our retrospective modeling do not support the existing hypotheses about the possible survival of multiple northern P.koraiensis refugia during the LGM in remote areas from the sea, such as on the Sikhote-Aline Mountains (Potenko and Velikov, 1998) and the Greater Khingan Range (Bao et al., 2015). However, the results confirm the possibility of such refugia within the Manchurian-Korean Mountains system(the watershed of the Yellow and Japanese Seas basins),along the coast of the Sea of Japan north to latitude 43°N (the Black Mountains massif, the eastern tip of the Manchurian-Korean Mountains), in the north of the Korean Peninsula – the Changbaishan Plateau, and in the northern part of the Baekdudaegan Ridge on the Korean Peninsula(Aizawa et al., 2003;Kim et al., 2005;Tong et al.,2020) (Fig.1b).

In the context of the known data on the high level of genetic diversity and the presence of unique alleles in the populations of P.koraiensis in the Changbaishan region(Chung et al.,2017),our results allow us to conclude that P. koraiensis has survived in refugia in the Manchurian-Korean Mountains since the end of the Pleistocene (Aizawa et al., 2012). At the same time, paleobotanical data indicate the Holocene expansion of the species within most of the modern range of the species on the mainland(Makinienko et al.,2008;Belyanin and Belyanina,2019),which includes the Sikhote-Alin mountain systems, the Lesser Khingan (Xiao Hinggan)and the southern tip of the Bureya Massif(Fig.1).

Modern P. koraiensis refugia in the mountainous regions of Honshu Island(Japanese Alps)within the framework of our model correspond to areas with optimal climatic conditions for the species during the LGM.Isolated populations of P. koraiensis currently existing in these mountainous areas are considered to be relics of the Late Pleistocene(Momohara et al., 2016) and have different mitotypes compared to continental populations(Aizawa et al.,2012).In the mountains of central Honshu, there are local refugia for other relict coniferous tree species(Okitsu and Momohara,1997).

Fig. 3. Overlap of potential distribution models of Pinus koraiensis of the three time frames: the past (LGM), the present, and the future 2070. (a) under the scenario of RCP2.6; (b) under the scenario of RCP8.5.

Thus,the results of our retrospective modeling are consistent with the previously proposed assumptions about the location of the northernmost Late Pleistocene P.koraiensis refugia on the mainland in the Manchurian-Korean Mountains system.Taking into account the observed distance of P. koraiensis seeds dispersal by nutcrackers (Nucifraga caryocatactes(Linnaeus 1758))to the distance of 4–5 km and the maximum recorded distance of 15 km for the closely related species Pinus cembra L.(Hutchins et al., 1996), even without considering the possibility of long-distance seed drift, the rate of dispersal of the species from Pleistocene refugia to current distribution limits with an annual speed of 550–600 m?year-1looks quite realistic.In contrast,the results obtained do not support the hypothesis of multiple P. koraiensis refugia in the mainland part of the range,including on the Sikhote-Alin and the Lesser Khingan ridges.

4.2. Predictive modeling

The predicted climate change in NA by 2070 will not cause a significant change in the territory of the potential climatic niche of P.koraiensis relative to the modern niche.However,the changes of the configuration of the territory with optimal climatic conditions for P. koraiensis do not suggest a straight pattern, such as the “disappearance in the south and appearance in the north”. The simulation results demonstrate a high probability of preservation of climatic refugia in the modern centers of genetic diversity of the species in the southern part of the range–on the Korean Peninsula and Honshu Island (Fig. 1c and d). Nevertheless,climate change in these areas will and is already causing an increase in the timberline in the mountains and eliminating tree species at lower altitudes(Koo et al., 2015;Du et al.,2018).

Both pessimistic and optimistic scenarios of climate change provide evidence of the disappearance of optimal conditions for P.koraiensis near the western edge of its distribution – in the southern part of the Bureya Massif and on Lesser Khingan. A less suitable climatic situation will develop in some areas of the south of the Russian Far East.

Considering possible changes in the range of P. koraiensis, it is worth noting that precipitation in the form of rain appears to be a major factor controlling the distribution.In this aspect,our model does not agree with the modelof Shitara etal.(2021),inwhichtheminimum temperatureof the coldest month is the most significant climatic parameter in the distribution of P.koraiensis.At the same time,the results of instrumental measurements show that the decrease in radial increments of P.koraiensis(Yu et al.,2011;Wang et al.,2013,2019;Cao etal.,2019;Ukhvatkina etal.,2021)andother coniferous tree species(Yu and Liu,2020)in NA is often associated with a lack of soil moisture, and this factor depends mostly on the amount of precipitation and the intensity of moisture evaporation.

4.3. Estimated distribution of P. koraiensis compared to other studies

Compared with the recently proposed species distribution model for P.koraiensis(Shitara et al.,2021),we did not use the standard Wordclim 1.4 indices (BIO1 – BIO19) but rather used the five most significant bioclimatic indices that take into account the heat supply of the territory and the distribution of precipitation in the form of snow and rain throughout the year(Noce et al.,2020).We also built our model on top of an RF classifier,which,in turn,could be considered a more flexible class of approaches than those based on the maximum entropy principle,especially in the case of complicated distributions in factor spaces (Mi et al.,2017;Santini et al.,2017).

The model presented in this study shows that the islands of Hokkaido and Sakhalin and the southern islands of the Kuril Archipelago are not suitable for P. koraiensis climatically and reflects the real situation. Our model,compared with the model of Shitara et al.(2021),is built on more exact distribution data obtained from the continental sector of the P.koraiensis distribution range and much closer describes the geographical range of species shown on vegetation maps of Russia and China(Sochava,1967;Su et al.,2020),being based on terrestrial observations and official records of the DBCF.

5. Conclusions

The distribution model of P. koraiensis, designed according to ecologically significant predictors, such as warmth and coldness Kira’s indices, continentality index, and rain and snow precipitation indices using 207 unique presence points,showed good agreement with both the species observations in nature and with the available cartographic materials built on the basis of forest revisions in the continental part of the specie’ range. The amount of precipitation in the period with positive average monthly temperatures appears to be the most significant predictor of the P. koraiensis distribution. Insufficient precipitation in the form of rain and aridization in the future can affect the distribution of the species in the driest parts of the range.These conclusions are consistent with the results of dendroecological studies of P.koraiensis in China and Russia, thus signaling a decrease in growth rates with a shortage of available moisture. At the same time, even with the realization of the pessimistic RCP8.5 scenario by 2070, climate change will not lead to a complete disappearance of suitable habitats in the Pleistocene refugia of the species. Efforts should be made for the total preservation of genetic resources of such populations in situ and ex situ within the framework of international cooperation of all countries in which P. koraiensis grows wild.

Funding

Not applicable.

Availability of data and materials

The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Authors’ contributions

Tatyana Y.Petrenko:data collection,experimental design and writing the original draft. Kirill A. Korznikov: data collection, experimental design,supervision,visualization,and writing the original draft.Dmitry E.Kislov:conceptualisation,statistical analysis,editing the original draft and visualization. Nadezhda G. Belyaeva: editing the original draft and visualization. Pavel V. Krestov: project administration, review and editing the original draft. The author(s) read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Acknowledgements

We are grateful to Victoria Chilcote and Mark Chilcote for linguistic help.

Appendix A. Supplementary data

Supplementary data to this article can be found online at https://do i.org/10.1016/j.fecs.2022.100015.