Major Forest Changes in Subtropical China since the Last Ice Age

: In the subtropical zone of southern China, there was a considerable conversion of forests from deciduous to evergreen broadleaf in the early Holocene. However, the exact timing of this vegetation change and its relationship to climate are still unclear. We examined a high-resolution pollen record collected in the mid-subtropical zone and then performed a correlation with regional data to reconstruct the history of forest ecosystems since the last deglaciation. Our data show that the expansion of the evergreen plant component already occurred at low elevations during the last deglaciation. The subtropical mountain landscape was not recolonized by evergreen forests until the mid-Holocene at about 8.1 ka BP. Based on fossil pollen reconstruction and climate model simulation, we conclude that the primary increase in evergreen components of subtropical ecosystems was triggered by postglacial temperature increase, and that a complete conversion from deciduous to evergreen forest ecosystems did not occur until Holocene winter temperatures and seasonal temperature contrast reached a threshold suitable for the growth and persistence of evergreen tree species.


Introduction
Evergreen broadleaved forests (EBLFs) are now the dominant biome in the subtropical zone in east China with abundant endemic plant species and high biodiversity [1]. The EBLF extends from about 22-23 • N to Yangtze River at 31-33 • N [2,3]. During the Last Glacial Maximum (LGM), subtropical China was covered by temperate deciduous (TEDE) and cool mixed (COMX) forests in a cooler climate [4,5]. Reconstructed LGM biomes suggest that tropical forests were absent from mainland southern China, or their range was reduced to very limited areas [6]. EBLF and warm mixed forest had retreated to tropical latitudes [7]. Therefore, it is likely that the EBLF was restricted to the southernmost part of China, between 6 and 8 • N south of its present range [6][7][8].
However, exactly when and how the EBLF recolonized subtropical regions since the last deglaciation is poorly known [5,9]. Moreover, the climatic interpretation of forest conversion remains controversial. Although long-term global warming since the last deglaciation is generally considered to be the main cause of forest changes and shifts in zonal vegetation boundaries [10,11], some have claimed that summer rainfall might have played a dominant role [10]. The mechanism of low latitude vegetation response to global Forests 2021, 12, 1314 2 of 10 climate change driven by ice sheet mass and/or insolation is also debated [12,13]. In this study, we examine a high-resolution pollen record of SZY peat bog in the mid-subtropical zone. This pollen record allowed us to reconstruct the timing of the conversion from the LGM deciduous biome to the now dominant EBLF. Climate thresholds for EBLFs were derived from its modern distribution using a species distribution model and past climate variables from transient climate simulations. A reconstruction of biome turnover and climate data helped us identify the key driving climate variables and the effects of postglacial global warming on EBLF's recolonization of evergreen ecosystems in subtropical China.

Regional Setting
Study core SZY (119.034 • E, 26.777 • N; 1007 m) was collected from a mountain peat bog in Jiufeng Mountains, Fujian Province, southeastern China ( Figure 1). The summits of nearby mountains are between 1300 and 1600 m. The area is located in the southern part of the mid-subtropical zone, which is strongly influenced by the East Asian Summer Monsoon (EASM). The local mean annual temperature is about 16 • C, and the mean January temperature and daily minimum are 6-7 • C and −4 • C, respectively, at about 1000 m elevation. The average altitude constant is about 0.59 • C/100 m. Annual precipitation varies between 1700 and 1900 mm, with most precipitation falling between late March and July. Evergreen broadleaved forest is the main vegetation type, which can be divided into the following five communities: coniferous/broadleaved mixed forest, warm coniferous forest, evergreen broadleaved forest, bamboo forest and shrubland. In the regions below 1000 m altitude, there are few natural forests left due to heavy human impact. Dominant arboreal plants at the genus level in the local natural or secondary forests of the mountains are Pinus, Castanopsis, Cyclobalanopsis, Cunninghamia, Altingia, Cryptomeria, Cinnamomum, Machilus, Schima and Eleocarpus. The common shrub taxa include Eurya, Neolitsea, Lindera, Ardisia, Rhododendron, Prunus, etc. The temperate deciduous and conifer elements, such as Fagus, Quercus, Carpinus, Acer, Tsuga, Rhododendron, Symplocos, are generally found at higher altitudes of 1200 m.

Coring, Sampling and Chronology
The SZY core was taken in 2011 using a Russian Corer. In this study, the upper 273 cm of the peat core (425 cm in total) was investigated. The uppermost sediment (0-32 cm) consists of marsh clay with abundant plant debris and high water content. The main part of the core below 32 cm consists of homogeneous brown-black peat and dark clay. For the upper 273 cm of core, a total of 180 samples (including 132 previously published [5] and 48 newly added) were taken at 1-2 cm intervals, corresponding to a temporal resolution of about 80 years between samples. The age-depth model was constructed using radiocarbon dating with accelerator mass spectrometry (AMS 14 C). Details of radiocarbon samples can be referred to [5]. Based on the calibrated age, the interval extends to 15,000 cal yr BP at 274 cm. The depth-age plot [5] shows a higher sedimentation rate during the Holocene (2.2 cm/100 yr) and a lower rate during the last deglaciation (0.83 cm/100 yr).

Pollen Analysis
The pollen preparation method basically followed the traditional laboratory treatment techniques [14] with some improvements such as the use of hydrofluoric acid. Each sample was weighed to 2-3 g of sediment, and one tablet of trace Lycopodium spores (27,637 per tablet) was used for calculation of pollen concentration (grains/g). The main treatment procedure included the hydrochloric acid (10%) to remove carbonates, hydrofluoric acid (40%) to remove silicates and potassium hydroxide to remove organic matter. The palynomorphs were separated from mineral sediment using standard heavy liquid flotation techniques [15]. In general, more than 400 pollen grains were counted, and pollen percentage calculation was based on sum of total pollen excluding aquatic taxa and all spores.

Biomization
In order to better illustrate temporal vegetation changes, the pollen taxa were grouped as plant functional types (PFTs). Based on which, the biomes were reconstructed by applying the biomization technique. We basically applied the existing schemes of plant functional types (PFTs) previously used for China [7], except for the exclusion of PFTs from high latitude arctic and boreal zones, and some moderate modifications to fit the regional biomes. In this study, only the following six biomes occurring in the study area were reconstructed: subtropical evergreen broadleaved forest (EBLF), tropical rain forest (TRFO), temperate deciduous broadleaved forest (TEDE), cold mixed forest (CLMX), cool mixed forest (COMX) and alpine shrub and meadow (ALPM). We named the biome of subtropical (warm-temperate) evergreen broadleaved forest as 'EBLF', which is approximately the same with the warm evergreen and mixed forest (WAMF) used in [7] and adopted the biome of ALPM proposed by [4] to identify subtropical mountain vegetation.

Species Distribution Modelling
By associating modern presence data and climate data, the species distribution model (SDM) [16] predicts the climatic envelope of the study species or vegetation type and projects it to different time periods to obtain potential distributions. Using this method, the potential distributions of EBLF at the time slices of LGM and the present were predicted. The presence data of EBLF was extracted from a digitalized vegetation map (Supplementary Figure S1). To avoid sampling bias, we divided the distribution range of presence data into several 1 • × 1 • grids, and randomly picked up only one point from each grid to obtain a filtered subset of presence data (n = 129) using the GridSample function of the R 'dismo' package. Seven key parameters among 19 bioclimatic variables representative of climate conditions from the WorldClim database were used for model calibration (Supplementary  Table S1). For the LGM climate conditions, we relied on the results estimated using an Earth System Model, MIROC-ESM [17]. The distribution of EBLF at each time slice was finally predicted using one of the most popular SDM methods, the maximum entropy algorithm (Maxent) within the R 'dismo' package with default settings based on the filtered presence data and climate data. The resulting habitat suitability, the value of which ranges from 0 (unsuitable) to 1, was used to indicate the potential distribution of EBLF. By comparing the predicted distribution model for the present time and the known distribution of EBLF, the model was validated with a threshold-independent measure of the area under the receiver operating characteristic curve (AUC). The highest value of 1 (0 being the lowest) represents a complete match between the modeled distribution and the known distribution.

Transient Data Analysis
TRACE21 is the first state-of-the-art transient simulation of global climate change over the last 21 ka at the T31 spatial resolution [18]. Two climatic variables including mean winter (December-February) temperatures (corresponding to bioclimatic variable 11 in WorldClim) and temperature seasonality (July-January) were used for comparison. We first extracted the temperature values of the grid where SZY is located, and then uniformly correct the simulated values by substituting in the bias caused by the altitudinal lapse-rate of modern climate observation in the study site.

Results
The high resolution and well dated record of SZY shows that the vegetation of the last 15,000 years passed through four main phases characterized by four distinct pollen zones ( Figure 2). During the last deglaciation between 15 and 11.5 ka BP, regional vegetation was dominated by cold-tolerant deciduous forest. The constituent deciduous taxa are Fagus, Quercus, Alnus, Betula, Corylus, Carpinus and others. Evergreen taxa such as Cyclobalanopsis and Castanopsis were also present in this pollen zone, but the percentages are generally lower than the subsequent Holocene values. Conifers such as Tsuga and Abies and the shrub Rhododendron (Ericaceae) were more abundant during the time interval earlier than 11.5 ka BP. Herbaceous taxa including Artemisia and Poaceae fluctuated and sometimes reached 5% (Figure 2). During the Holocene, at least three zones of pollen assemblage can be identified. The lower zone is from 11.5 to 8.2 ka BP when Pinus and Taxodiaceae (possibly Chinese watercypress Glyptostrobus) and wetland sedge increased rapidly. Alnus was also relatively abundant, the species of which may belong to the wetland environment taxon in the region (i.e., A. trabeculosa). The middle zone dated to 8.1-0.8 ka BP corresponds to the pollen  During the Holocene, at least three zones of pollen assemblage can be identified. The lower zone is from 11.5 to 8.2 ka BP when Pinus and Taxodiaceae (possibly Chinese water-cypress Glyptostrobus) and wetland sedge increased rapidly. Alnus was also relatively abundant, the species of which may belong to the wetland environment taxon in the region (i.e., A. trabeculosa). The middle zone dated to 8.1-0.8 ka BP corresponds to the pollen assemblages marked by evergreen broadleaved taxa mainly including Cyclobalanopsis, Castanopsis and Altingia. The wetland sedge pollen of Cyperaceae was extremely abundant in the lower part of this zone, which was followed by an increase in Poaceae and fern spores. The upper most zone corresponding to the last~0.8 ka BP shows an abrupt increase in Poaceae, Artemisia, Dicranopteris and Pinus, which are considered as the pollen taxa associated with human activities [19]. The total non-arboreal pollen percentage is up to 52% in the uppermost part of the record.
The biome reconstruction result shows a strong dominance of temperate deciduous broadleaved forest (TEDE) between 15 and 11.5 ka BP, and the biome with the second highest score was the cool mixed forest (COMX). An obvious decrease in all biomes except alpine shrub and meadow (ALPM) happened at around 11.5-12.5 ka BP. During the transition from early to mid-Holocene (ca. 11.5-8.2 ka BP), there was a gradual increase in subtropical (warm-temperate) evergreen broadleaved mixed forest (EBLF), although its scores were still lower than those of TEDE.
The onset of the EBLF increase after the obvious biome change at ca. 11.5 ka BP corresponds to the cold Younger Dryas (YD), suggesting that the EBLF's recolonization of the subtropical regions resumed and accelerated after the YD termination. After 8.1 ka BP, EBLF became dominant over TEDE, with several fluctuations of abundance (e.g., at around 3.2 ka BP) until the end of the Holocene (Figure 3).  The species distribution model was used to predict the potential distributions of evergreen broadleaved forest at the LGM and in the present (Figure 4). The predicted present-day distribution of EBLF agrees well with the known distribution (AUC = 0.94), indicating that the model is reliable for drawing further conclusions. The simulated LGM distribution shows that EBLF migrated southward into northern Indochina in the western part but found potentially suitable conditions (in situ refugia) in the central east (i.e., some large river basins between the Wuyi and Xuefeng mountain ranges). The species distribution model was used to predict the potential distributions of evergreen broadleaved forest at the LGM and in the present (Figure 4). The predicted present-day distribution of EBLF agrees well with the known distribution (AUC = 0.94), indicating that the model is reliable for drawing further conclusions. The simulated LGM distribution shows that EBLF migrated southward into northern Indochina in the western (CLMX), cool mixed forest (COMX) and alpine shrub and meadow (ALPM). The species distribution model was used to predict the potential distributions of evergreen broadleaved forest at the LGM and in the present (Figure 4). The predicted present-day distribution of EBLF agrees well with the known distribution (AUC = 0.94), indicating that the model is reliable for drawing further conclusions. The simulated LGM distribution shows that EBLF migrated southward into northern Indochina in the western part but found potentially suitable conditions (in situ refugia) in the central east (i.e., some large river basins between the Wuyi and Xuefeng mountain ranges). The quantile result of critical climate thresholds for the simulated present distribution is shown below. The suitability value ranges from 0 (unsuitable) to 1. The green to blue distribution ranges represent the predicted most probable climatic-suitable habitats for the EBLF.

Substantial Forest Transformation during the Holocene
Our data show that the vegetation in subtropical mountains around 1000 m of elevation during the last deglaciation was dominated by deciduous forest (TEDE) and conifers (COMX). The presence of evergreen components might reflect a kind of mixed forest in which deciduous trees predominated, but there were also a significant number of evergreen taxa. It is also possible that some of the pollen of the evergreen plants was transported from lower altitudes by updraft wind, since the evergreen forest belt in the lower mountains was not very far from the study area during the last deglacial. However, considering the small size of our study site, we suppose that the pollen assemblage of our record largely reflects local vegetation composition [20]. The main transition from a deciduous to an evergreen biome occurred rapidly over a narrow time span at ca. 8.1 ka BP (Figure 3). However, this timing of biome turnover raised controversy as pollen records from other latitudinal positions show different timing. In Taiwan Island, the expansion of subtropical (and/or warm temperate) evergreen forest took place much earlier at about 11.5 ka BP [4]. Fossil pollen data over eastern China show that the timing of vegetation change is not simultaneous, but rather follows a time-transgressive pattern from the south (10.3 ka BP) to the north (6 ka BP) [11].
Effectively, the complete absence of the evergreen broadleaved forest during the LGM in the subtropical zone of the Chinese mainland inferred by previous biome reconstruction [6] is under debate as phylogeographic studies have suggested that species of the EBLF survived the LGM locally in numerous northern refugia [21]. Our pollen record shows that the evergreen forest may have occupied the lower altitude regions even during the glacial periods, taking into consideration the regional altitudinal air temperature lapse rate, which may have steepened during the LGM [22]. Our simulated distribution of evergreen forest shows a shrinkage of evergreen forest distribution, which also persisted in nearby lowland areas as refugia such as the basins between the Wuyi and Xuefeng mountain ranges in southeastern China. Such model simulations imply that the local mountain evergreen forest experienced a complex migration upward from lowland during the early to mid-Holocene, agreeing well with the refugia theory proposed by phylogeographic studies [21,23], although there may have also been some northward shift of constituent evergreen species from the south to the north [21].

Holocene Climatic Changes as a Cause of Subtropical Forest Transformation
The rapid EBLF/TEDE biome change at about 8.1 ka BP suggests that the rate of vegetation change increased as temperature increased and crossed a threshold that delineates the climatic range of the two biomes. Mean minimum temperature and seasonal temperature contrast are commonly considered as effective climate indices for defining the upper limit for the evergreen broadleaved forest [24], as photosynthetic efficiency leads to a reduction in evergreen broadleaved species [25]. From the analysis of Holocene climate changes in the region ( Figure 5), we further confirm that the seasonal temperature difference and/or the coldest climate condition must be the most important climatic factors controlling the switch between evergreen/deciduous trees, and that it is much less related to the change of rainfall amount, because both the moisture proxy from stalagmites [26] and the simulated precipitation [27] show a maximum amount centered around 10 ka BP, which contradicts the timing of our pollen records. . Temperature and seasonal difference as the key factors of evergreen forest turnover during the Holocene. Temperature indices including brGDGT-based mean annual temperature reconstruction [27], TRACE21 [18] modeled winter temperature (Twin) and seasonal mean temperature difference between July and January. The values of the two variables were revised from the regional output according to the local air temperature (SZY).
The local air temperature was estimated by the proxy of branched glycerol dialkyl glycerol tetraether (brGDGT) [28], which shows synchronous changes with the reconstructed vegetation changes from the same site of SZY. There was a rapid increase in the mean annual temperature (MAT) from 13 to 16 °C at about 8 ka BP. The temperature at about 16 °C corresponds to the most appropriate MAT range (14.5-18.5 °C) for the present evergreen forest based on the simulated distribution (Figure 4). This suggests that a local forest change from a deciduous-dominated forest to an evergreen-dominated forest may have been triggered by the temperature change up to the threshold. Note also that the . Temperature and seasonal difference as the key factors of evergreen forest turnover during the Holocene. Temperature indices including brGDGT-based mean annual temperature reconstruction [27], TRACE21 [18] modeled winter temperature (Twin) and seasonal mean temperature difference between July and January. The values of the two variables were revised from the regional output according to the local air temperature (SZY).
The local air temperature was estimated by the proxy of branched glycerol dialkyl glycerol tetraether (brGDGT) [28], which shows synchronous changes with the reconstructed vegetation changes from the same site of SZY. There was a rapid increase in the  Figure 4). This suggests that a local forest change from a deciduous-dominated forest to an evergreen-dominated forest may have been triggered by the temperature change up to the threshold. Note also that the subsequent climate fluctuations after 8 ka BP to the threshold MAT value of about 14.5 • C (e.g., at 3.2 ka BP) repeatedly caused a switch between evergreen/deciduous biomes. This further supports the significant role of temperature threshold in the biome turnover. The vegetation changes in the studied pollen record are closely related to the output of the TRACE21 climate model [18], e.g., the trends in mean winter temperature (Twin) and the seasonal temperature difference between July and January (STdif) from the last deglaciation to the Holocene ( Figure 5). The modeled winter temperature in the subtropical climate zone shows a gradual change from 2 to 5 • C in the early Holocene at 11.5-8.1 ka BP (more abrupt since 9 ka BP). Similarly, there was a gradual decrease in the seasonal temperature difference at 11.5-8.1 ka BP from 24.5 to less than 23 • C ( Figure 5). Note that these values were corrected to the altitude of the SZY record based on the local average elevation lapse rate. These values at 8.1 ka BP agree very well with the present-day range of the most suitable temperature for evergreen trees (Twin: 5.5-10.5 • C, STdif: 14-23.5 • C). They are also in broad agreement with the meteorological data delineating the modern evergreen broadleaved forest in the subtropical mountains [24,25]. Therefore, there is no doubt that evergreen trees became dominant in the forest when the winter and seasonal thermal conditions reached their thresholds. After the YD at 11.5 ka BP, progressive winter warming and decreasing seasonality led to an increase in evergreen broadleaved elements in the local forest, and the subsequent onset of the Holocene optimum at~8.1 ka BP forced the general forest turnover.
According to our record, we suggest that the forest change from deciduous to evergreen composition was initiated after the YD in the subtropical region. The greatest extent of EBLF in our study area corresponds to the Holocene optimum, which begins at~8 ka BP and shows a close relationship with the increase in the winter temperature and the decrease in the seasonal temperature difference to thresholds. Our data are consistent with winter climate dynamics influencing the strength of the Asian summer monsoon [29]. This study provides an improved understanding of the ecologically sensitive region at the interface between the warm temperate and tropical bioclimatic zones. Ice sheet melting exerts an important influence on low latitude climate through associated changes in the global sea level [30], and the sea surface temperature of Indo-Pacific warm pools shows a progressive warming since the beginning of the Holocene [31]. The observed amplitude of subtropical ecosystem changes during global climate warming over the last postglacial period draws our attention to the period when one extended biome with new species replaced another. The biome turnover we identified occurred over a time span comparable to the few decades of ongoing global climate change. Thus, ongoing climate change will likely have a significant effect on ecosystems, especially upland species [32], although the impact is expected to be less significant in upland areas than in lowland areas [33]. Many endemic species, e.g., Fagus, Cryptomeria, Tsuga, Pseudotsuga, Quercus, in subtropical Asia are near threatened, vulnerable, endangered or even critically endangered according to the IUCN Red List [34]. These species are all fragmented and scattered within the evergreen biome. Ongoing warming is likely to lead to their extinction, resulting in a depletion of biodiversity in Asian ecosystems, which deserve the establishment of new national parks and biosphere reserves.

Conclusions
This study shows that both ecosystem composition and their range have undergone major changes since the last ice age, especially in the last thousand years. The observed biome transformation in our study site appears to have taken place within a time span comparable to that of ongoing climate warming, providing important knowledge to better Forests 2021, 12, 1314 9 of 10 understand ecosystem transformations and manage their future changes. Many endemic evergreen tree species are now threatened with extinction, which also raises many questions regarding the managing of diversity in subtropical ecosystems.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/f12101314/s1, Figure S1: The distribution of the presence points extracted from digitalized EBLF distribution used for predicting the modern distribution of EBLF, Table S1: Estimates of relative contributions of the environmental variables to the Maxent model.