Different processes shape the patterns of divergence in the nuclear and chloroplast genomes of a relict tree species in East Asia

Abstract Isolation by spatial distance (IBD), environment (IBE), and historical climatic instability (IBI) are three common processes assessed in phylogeographic and/or landscape genetic studies. However, the relative contributions of these three processes with respect to spatial genetic patterns have seldom been compared. Moreover, whether the relative contribution differs in different regions or when assessed using different genetic markers has rarely been reported. Lindera obtusiloba has been found to have two sister genetic clades of chloroplast (cpDNA) and nuclear microsatellite (nSSR), both of which show discontinuous distribution in northern and southern East Asia. In this study, we used the Mantel test and multiple matrix regression with randomization (MMRR) to determine the relative contributions of IBD, IBE, and IBI with respect to L. obtusiloba populations. Independent Mantel tests and MMRR calculations were conducted for two genetic data sets (cpDNA and nSSR) and for different regions (the overall species range, and northern and southern subregions of the range). We found a significant IBI pattern in nSSR divergence for all assessed regions, whereas no clear IBI pattern was detected with respect to cpDNA. In contrast, significant (or marginal) divergent IBD patterns were detected for cpDNA in all regions, whereas although a significant IBE was apparent with respect to the overall range, the effect was not detected in the two subregions. The differences identified in nSSR and cpDNA population divergence may be related to differences in the heredity and ploidy of the markers. Compared with the southern region, the northern region showed less significant correlation patterns, which may be related to the shorter population history and restricted population range. The findings of this study serve to illustrate that comparing between markers or regions can contribute to gaining a better understanding the population histories of different genomes or within different regions of a species' range.


| BACKG ROU N D
Disentangling the various contributions of geographical, environmental, and historical climatic dynamics to population genetic structures could provide significant insights into the species diversification (Lee & Mitchell-olds, 2011;Meirmans, 2012;Wang & Bradburd, 2014;Wright, 1943). Under the isolation-by-distance (IBD) scenario, gene flow is considered to be negatively correlated with geographical distance, and accordingly genetic differentiation between populations is expected to be positively correlated with distance (Hedrick, 2005;Meirmans, 2012;Wright, 1943). In contrast, the isolation-by-environment (IBE) hypothesis proposed that it is differences in environmental variables between sites that can explain most of the observed variance in gene flow and implies that local adaptation may occur and restrict the gene flow between sites (Manthey & Moyle, 2015;Sexton, Hangartner, & Hoffmann, 2014;Wang & Bradburd, 2014). In addition to IBD and IBE, historical climatic fluctuations, particularly glacial-interglacial cycles during the Pleistocene era (2.58-0.0117 million years ago, hereafter, Ma), have markedly influenced species distribution and demographic history (Qiu, Fu, & Comes, 2011;Soltis, Morris, McLachlan, Manos, & Soltis, 2006). These fluctuations in species distribution may subsequently lead to an instability in habitat connectivity and thereby contribute to enhancing population genetic differentiation. Isolation by instability (IBI) denotes a negative correlation between the extent of the historical instability of habitat connectivity and gene flow. Although both IBE and IBI can be calculated based on environment data, these processes differ in important respects, with IBE being characterized by the influence of differences in habitat environment on gene flow, whereas IBI is characterized by the influence of long-term habitat instability due to climate change and its resistance on dispersal between habitats. Previously, a number of studies have compared IBD and IBE to evaluate the contribution of local adaptation to present-day patterns in spatial genetic (Pelletier & Carstens, 2018;Sexton et al., 2014;Shafer & Wolf, 2013;Wang, Glor, & Losos, 2013 (Sexton et al., 2014), Quercus lobata Née (Gugger, Ikegami, & Sork, 2013) Hypsiboas lundii Burmeister, 1856 (Vasconcellos et al., 2019)].
Owing to its extensive physiographical heterogeneity, in conjunction with frequent historical changes in climate and sea level, East Asia currently harbors the world's highest diversity of temperate flora (Qian & Ricklefs, 2000). It has been suggested that Cenozoic (66-23 Ma) relict genera in East Asia can be subdivided into two refugial regions (Milne & Abbott, 2002), namely the northern refugial region (NEA), comprising Japan, Korea, and adjacent areas of Northeast China, and the southern region (SEA), comprising south and southeastern China with extensions into Himalayan areas. These two refugial regions are divided by a low-precipitation zone, referred to as the "aridity belt" (Figure 1), which developed during the Eocene  and has persisted to the present day (Chen, Deng, Zhou, & Sun, 2018;Donoghue, Bell, & Li, 2001;Milne & Abbott, 2002). The arid belt is believed to represent a climatic barrier that restricts gene flow between the NEA and SEA regions (Bai, Wang, & Zhang, 2016;Milne & Abbott, 2002;Ye et al., 2017). Although studies over the past two decades have begun to focus on the spatial genetic patterns of temperate plant species in both the NEA and SEA regions and their responses to historical climate change (Chen et al., 2018;Harrison, Yu, Takahara, & Prentice, 2001;Qiu et al., 2011;Tang et al., 2018), the relative roles of spatial arrangement (IBD), local adaptation (IBE), and historical climate fluctuations (IBI) in shaping spatial genetic patterns in East Asia have rarely been thoroughly examined (but see He, Wang, Li, & Yi, 2016;Zhang, Wang, Comes, Peng, & Qiu, 2016 (Pauwels et al., 2012), Juglans mandshurica Maxim. (Bai, Liao, & Zhang, 2010)].
Nonetheless, the examination of IBD-, IBE-, and IBI-related patterns in nuclear markers tends to be considerably more frequent than that using chloroplast markers. We assume that this disparity in usage can be attributed to the fact that IBD, IBE, and IBI have been studied to a greater extent in landscape genetics (in which chloroplast markers tend to be assessed less frequently) than in phylogeographic (Manel, Schwartz, Luikart, & Taberlet, 2003), which reflects the slight difference between the two fields (Rissler, 2016;Wang, 2010). Despite the disparate status of these markers, we maintain that comparisons among IBD, IBE, and IBI patterns based on both markers types would be necessary for gaining a more comprehensive understanding of the processes occurring in different genomes in the same cell, which in turn can contribute disentangling the complex processes associated with the origins of biodiversity.
As its designation indicates, the SEA is located in a lower latitude region than the NEA, and is also characterized by a higher degree of physiographical heterogeneity. Both pollen records and genetic data indicate that the SEA region has had the potential to sustain more refugia for temperate flora than the NEA region during glacial periods (Qian & Ricklefs, 2004;Qiu et al., 2011). Hence, we presume that IBD, IBE, and IBI processes would have differed to certain extents in the NEA and SEA regions. Although some independent studies have examined these processes in the SEA or entire region (He et al., 2016;Zhang et al., 2016), the regional differences between the NEA and SEA regions have yet to be compared. Moreover, even within the same region, different species with differing traits may show contrasting responses to IBD, IBE, and IBI. Accordingly, from the perspective of assessing the relative roles of these different processes, studies of widespread species with sister clades in the NEA and SEA regions would contribute to minimizing species difference and maximizing regional differences.
Lindera obtusiloba Blume is a drought-sensitive deciduous shrub or small tree that has a widespread distribution in both the NEA and SEA regions (Figure 1). In a previous study on this species (Ye et al., 2017), both chloroplast DNA sequence (cpDNA) and nuclear microsatellite (nSSR) data indicated the occurrence of two clearly discontinuous genetic clades in the NEA and SEA regions ( Figure 1). Over the entire distribution range of L. obtusiloba, the greatest divergence in the populations of this species occurs between the NEA and SEA populations. Both the most recent common ancestor of cpDNA haplotypes and population divergence estimates based on nSSRs have dated the divergence between NEA and SEA populations to the Pliocene (5.33-2.58 Ma;Ye et al., 2017). Notable in this regard is the fact that the genetic break between NEA and SEA populations of L. obtusiloba is considerably sharper than that reported for other species with a similar distribution, for example, Asian butternuts (Bai et al., 2016), Acer mono Maxim. (Guo et al., 2014), Kalopanax septemlobus Koidz. (Sakaguchi et al., 2012), and Prinsepia Royle (Ma, Wang, Tian, & Sun, 2019). We suspect that this abrupt genetic break between L. obtusiloba populations could be attributed to the drought sensitivity of the species and that the arid belt located between NEA and SEA regions may represent a strong climatic barrier to the spread of L. obtusiloba.
The findings of both field and common garden studies have indicated that the NEA and SEA populations of L. obtusiloba can be distinguished with respect to seed weight and multiple phenological traits, including bud break and fruit maturation (Dai, 2019). In spite of such trait differences, however, the two genetic clades are sister clades that are characterized by similar climatic niches (Tian 2019, unpublished data). Hence, L. obtusiloba would appear to provide an ideal system for comparative analysis of evolutionary processes occurring in the NEA and SEA regions. In this regard, although Ye et al. (2017) have previously provided a detailed description of the spatial genetic patterns of L. obtusiloba in the NEA and SEA regions, there has to date been no comparative examination of the relative contributions of IBD, IBE, and IBI to present-day distributions. Moreover, Ye et al. (2017) found that compared with SEA populations of L. obtusiloba, those in the NEA harbors a higher genetic diversity of nSSRs, although not of cpDNA, thereby implying the differential history of the two marker types. Accordingly, on the basis of a combination of cpDNA and nSSR data and ecological niche modeling (ENM), we examined IBD, IBE, and IBI patterns in both types of markers with respect to the overall range of L. obtusiloba and also that in the two divided subregions (NEA and SEA). Our goal was to identify potential differences in the evolutionary processes occurring in the different regions and for different types of genetic markers.

| Genetic differentiation
In the present study, we used a data set previously compiled by Ye  (Ye et al., 2017) are denoted by different symbols: NEA clades, red circle; SEA clades, red triangle. The different colored areas of the map indicate the distribution of annual precipitation. This indicates that the approximate range of the two Cenozoic relict flora refugial regions in East Asia (NEA and SEA) is characterized by annual precipitation greater than 1,000 mm (Milne & Abbott, 2002). The region of low precipitation dividing the NEA and SEA (precipitation <1,000 mm) is referred to as the "arid belt," which is suggested to constitute a climatic barrier for species migration between the NEA and SEA regions (Milne & Abbott, 2002). (c) The PCA results for all sampled populations, indicating bioclimatic differences between populations. The black triangles and red circles represent SEA and NEA clade populations, respectively the accession nos. KU645591-KU645758. Multiple population genetic distance indices were determined for both nSSRs and cpDNA, which constitute a comprehensive map of population genetic divergence (Appendix S2). For cpDNA, we calculated the population genetic distance indices N ST , G ST , and D xy using DnaSP version 6.12.03 (Rozas et al., 2017), and for nSSRs, the population genetic distance indices F ST and R ST were calculated using GenAlEx 6.5 (Peakall & Smouse, 2012) and SPAGeDi 1.5 (Hardy & Vekemans, 2002), respectively. All calculations for genetic distance indices were performed using the default settings in each program.

| Geographical distance
We represented the locations of the sampled populations in terms of longitude and latitude under the Beijing 1954 projection system.
Pairwise geographical distances between sampled populations were calculated by using the function "distm" in the package "geosphere" (Hijmans, 2015) of R version 3.6.0 (R Development Core Team, 2013). For distance computation, we applied the "distVincentyEllipsoid" function, which calculates the shortest circle distance between points.

| Environmental distance
As environmental variables, we used the data of 19 bioclimatic variables for the current period

| Resistance distance
To assess the influence of the historical instability of habitat connectivity on population genetic differentiation, we modeled the distributions of suitable habitats for the present, last glacial maximum (LGM; 21 thousand years ago, hereafter, kya), and last interglacial (LIG; 120 kya) periods based on ecological niche modeling (ENM) in MAXENT (Phillips & Dudík, 2008). The ENM procedure and the results obtained using this procedure have been described in detail in our previous study (Ye et al., 2017). In brief, 100 occurrence sites of L. obtusiloba showing no spatial redundancy were included in the ENM. The locations of these 100 sites and data sources can be found in Appendix S4. The modeling included only 11 bioclimatic variables that lacked strong correlation (r < .9), namely annual mean temperature, mean diurnal range, isothermality, temperature seasonality, maximum temperature of the warmest month, minimum temperature of the coldest month, mean temperature of the wettest quarter, annual precipitation, precipitation of the wettest month, precipitation of the driest month, and precipitation seasonality. Model validation was performed via 10 independent runs with 20% randomly selected data. High receiver operating characteristic (ROC) values (0.928 ± 0.025) have been reported for this type of modeling, which indicate its good accuracy. The established model was projected onto LGM (MIROC-ESM) and LIG climatic models (Shields et al., 2012;Watanabe et al., 2011).
The climatic suitability of each grid in each period was described in terms of the species occurrence possibility obtained from the ENM results for L. obtusiloba. To calculate resistance distance, the climatic suitability of each period was converted into a resistance value as climatic unsuitability ("1-suitability"), using the "create friction layer" function in the SDM toolbox (Brown, 2014;Brown, Bennett, & French, 2017

| IBD, IBE, and IBI analyses
To investigate the relationships between ecological factors (geographical, environmental, and historical habitat instability) and population genetic differentiation, we performed Mantel tests To compare the processes shaping the spatial genetic structure in the NEA and SEA regions, we performed the same Mantel test and MMRR analyses using only data pertaining to the NEA or SEA, respectively.
When performing analyses, we did not take into consideration nonlinear relationships. Although certain variables probably showed nonlinear relationships with genetic distance, their linear relationships were still significant in this study. Moreover, whereas a linear fit may not the best fit for these relationships, from the perspective of simplification, we only included linear relationships for variables when comparing the relative contributions of IBD, IBE, and IBI to spatial genetic patterns. Therefore, the use of a linear model may represent a more conservative test of the relationships.

| Environmental analyses
The results of PCA indicated that the first two coordinates (PC1  (Figure 1c).
ENM results (Ye et al., 2017) indicated that the area of habitats in the SEA with suitable climates was larger during the LGM than the present-day distributions, whereas the area in the NEA was more climatically limited. Furthermore, during the LIG, the area of suitable habitats in both NEA and SEA regions was smaller than that during the present and LGM periods.
Climate instability resistance values, which were obtained based converting the ENM results and summing over the three periods, showed an uneven spatial distribution (Figure 2). We found that the central SEA region and the Japanese Archipelago in the NEA had larger continuous ranges of long-term stable habitats than the remaining areas of species' range, which indicates a lower climate instability resistance for long-term gene flow in these two regions.
Interestingly, the population genetic differentiation of nSSRs was found to be generally lower in the central SEA region and Japanese Archipelago (Figure 2a), whereas no comparable pattern was observed with respect to cpDNA (Figure 2b).

| IBD, IBE, and IBI in cpDNA
Our analyses of the chloroplast genome revealed differences in the correlation patterns of the three genetic distance indices throughout the entire distribution range of L. obtusiloba. Using the Mantel test, we found that genetic distance D xy , which indicates divergence since the last common ancestor between populations, was significantly correlated with the three factors (Figure 4), with the highest correlation value observed for geographical distance (r = .4108, p = .0001). MMRR analysis revealed that D xy was only influenced by geographical distance (β = 0.0006, p = .0007; Table 1).

F I G U R E 2
The geographical distribution of population genetic distances for (a) the nSSR index R ST and (b) the cpDNA index D XY . The population genetic distance of each sampled populations is represented by the average R ST (or D XY ) with that of each of the other populations. A higher than the species-level average R ST (or D XY ) value is represented by a black circle, whereas a lower than average value is denoted by a while circle. The size of the circles is proportional to the difference from the average value. The background map in both (a) and (b) is the historical climatic instability resistance map, with 0 representing the lowest climate instability resistance due to the highest habitat suitability in each period [last interglacial (LIG), last glacial maximum (LGM), and the present), and 3 representing the highest climate instability resistance due to the lowest habitat suitability in each period  TA B L E 1 Multiple matrix regression with randomization (MMRR) analysis results showing the relative contribution of isolation by spatial distance (IBD), environment (IBE), and historical climatic instability (IBI) with respect to multiple population genetic distance indices of both nSSRs and cpDNA markers in the overall range (overall) of Lindera obtusiloba and two subregions (NEA and SEA) within this range

| Different processes underlying genetic differentiation in the nuclear and chloroplast genome
In this study, we identified differences in the processes underlying the divergent genetic patterns in nSSRs and cpDNA in L. obtusiloba.
With respect to nSSRs, we detected significant IBI and IBD patterns over the entire species range, with the contribution of IBI being found to be greater than that of IBD. Given that neither of the NEA and SEA subregions show a clear IBD pattern, we infer that the observed IBD pattern over the entire range may have been primarily influenced by the greatest population divergence between the NEA and SEA regions, which has been reported previously by Ye et al. (2017). Thus, climatic instability since the late Pleistocene is likely to have been the main driving force underlying the nuclear genetic patterns observed in L. obtusiloba.
Our observations indicate that climate instability resistance between populations would not appear to explain the genetic divergence of cpDNA. However, the genetic distance index D xy shows a clear IBD pattern, and the genetic differentiation of the indices N ST and G ST shows a clear IBE pattern over the entire range of the species. We suspect that local adaptation to different environments in the NEA and SEA regions ( Figure 2) may have driven chloroplast divergence in the population. Generally, the NEA region is considered to have greater seasonality than the SEA region with respect to temperature and precipitation (Bai et al., 2016;Qiu et al., 2011).
Although we were unable to detect any significant differences between the habitats of L. obtusiloba in NEA and SEA regions in relation to humidity (PC1) nor heat (PC2), the habitat types in the NEA are unique when these two climatic factors are considered together. We suspect that habitat in the NEA region is a marginal type of that in the SEA region (Figure 1c studies (Dai, 2019). Compared with IBD or IBI, the IBE patterns detected in the genome can be heterogeneous, as these patterns tend to be driven primarily by selection against immigration (Sexton et al., 2014;Wang & Bradburd, 2014). Therefore, analysis of a limited number of nSSR loci may not be sufficient to detect a significant IBE pattern in the nuclear genome.
The heredity of nuclear and chloroplast genomes may have contributed to the contrasting responses of these genomes to historical climate instability. The chloroplasts in our study species are maternally inherited and can only be transmitted to offspring by seed flow, whereas the nuclear genomes are biparentally inherited and can be transmitted via both seed and pollen flows. Dispersal of the seeds and pollen of L. obtusiloba is mediated by birds and insects, respectively (Smith, Hamel, Devall, & Nathan, 2004), and previous phylogeographic studies and fossils data indicate that climate change can affect the distribution of these vector species (Dong et al., 2017;Kharouba, Lewthwaite, Guralnick, Kerr, & Vellend, 2019;Wang, Zhu, Heller, Zhou, & Shi, 2020). Hence, it is to be assumed that dispersal based on these species would be influenced to a considerable extent by climate fluctuations (Guisan & Thuiller, 2005;Pearson & Dawson, 2003;Van der Putten, Macel, & Visser, 2010). Moreover, compared with birds, insects typically have a limited range of dispersal (Phillipsen et al., 2015;Phillipsen & Lytle, 2013). Accordingly, given that patterns in the nuclear genome can reflect the dispersal patterns of both seed and pollen, the effect of climate instability on bird and insect distribution may have resulted in a stronger IBI pattern in the nuclear genome compared with that of the chloroplast.
In addition to the contrasting modes of heredity, differences in ploidy may also have contributed to the different correlation patterns we detected for the two genomes. The nuclear genome is diploid, undergoes recombination, and generally has a larger effective population size than the chloroplast, which is haploid and does not recombine (Gaudeul, Gardner, Thomas, Ennos, & Hollingsworth, 2014;McCauley, 1994). Hence, the nuclear genome might be excepted to have a higher mutation rate and genetic diversity than the F I G U R E 4 Scatter plots showing the relationships between the population genetic distance index D XY (cpDNA) and (a) geographical distance, (b) environmental distance, and (c) climatic instability resistance distance. The statistical results of Mantel tests are shown in the upper right-hand corner of each plot chloroplast, which could reflect the greater effect the more substantial influence of recent historical climate change on the gene flow between populations (Bai et al., 2010). In addition, nSSRs, which are markers commonly used in phylogeographic studies, have a considerably higher mutation rate than cpDNA sequences, which may enhance the observed IBI pattern. In summary, differences in heredity and ploidy between the nuclear and chloroplast genomes or the differences between nSSR and cpDNA sequences could account for the stronger IBI pattern in nSSRs that we observed in the present study.  (Ye et al., 2017). The combination of climate instability and geological isolations would thus have led to a negative correlation between gene flow and resistance or geographical distance, thus resulting in significant IBI or IBD patterns.

| Genetic
Although molecular data for the NEA region indicate long-term persistence, most populations, particularly those on the mainland, were formed recently via postglacial migration from the southern margins of the NEA or Japanese Archipelago. Inferences with respect to NEA and SEA glacial refugia are also supported by the distribution of climatically stable habitats during the late Pleistocene ( Figure 2). These findings indicate that most NEA populations have a more recent population history than those in the SEA region. In addition, the southern clades in the SEA region have a considerably wider distribution than clades in the NEA region. Although the number of populations we sampled in the NEA (10 populations) was smaller than the number examined in the SEA region (18 populations), we do not believe this to be a prominent factor accounting for the less pronounce genetic patterns detected in the NEA. Relative to their respective range sizes, sampling densities in the two regions were similar, and thus, the influence of sampling effort on characterizing spatial variance in the two regions would be comparable. Insufficient time and/or space can restrict the effects of evolutionary forces (such as drift or selection) and thereby reduce the significance of detected patterns (Holderegger, Buehler, Gugerli, & Manel, 2010;Sexton et al., 2014). In the NEA region, the shorter evolutionary history and relatively homogeneous topographical have probably been important factors contributing to limited drift and geographical isolation. Moreover, postglacial expansion occurred in a restricted population range in the mainland area of the NEA region, which may also have contributed to limited population genetic divergence. Hence, the less prominent IBI or IBD patterns we identified in the NEA region may be related to the recent population history and restricted population range of L. obtusiloba in this region.

| CON CLUS IONS
In this study, we compared the relative contributions of IBD, IBE, and IBI to the genetic divergence of L. obtusiloba in East Asia, based on analyses of different types of genetic marker and different regions.
We found that different processes have probably driven the patterns of nuclear and chloroplast genetic divergence. The divergence of nSSRs appears to have been influenced primarily by the historical climatic instability (IBI) over either the entire species range or within the two distinct subregions (NEA and SEA) of this range. In contrast, we were unable to detect any significant IBI pattern with respect to the divergence of cpDNA. We envisage that changes in climatic change during the late Pleistocene have had a stronger influence on population divergence of the nuclear genome than that of the chloroplast genome, which can probably be attributed to differences in the inheritance mode and ploidy of these two genomes. Differences in the patterns of cpDNA divergence are assumed to be influenced by spatial distribution (IBD). In addition, local adaptation may have affected the divergence of the chloroplast between the northern and southern regions, whereas we detected no clear IBE pattern for nuclear divergence based on the set of nuclear loci analyzed in the present study.
We also identified convergence in the correlation patterns of both the nuclear and chloroplast genomes in the NEA and SEA regions. However, due to the shorter population history and more restricted population range of L. obtusiloba in the NEA, we found this region to be characterized by a less significant pattern of genetic divergence than the SEA region. Collectively, the findings of this study serve to emphasize that even within a single species, analysis of different genetic markers or different regions can reveal differences in IBD, IBE, or IBI patterns or differing extent of significance. Such findings will contribute to gaining a more complete understanding of the population history of different genomes or patterns of divergence in different regions of a species' range.

ACK N OWLED G M ENTS
We are grateful to Hai-Tao Yang, Shun Li, and Lu-Ye Shi for their constructive comments regarding the manuscript. We are also indebted to the two anonymous reviewers for their insightful comments and Wiley editing services for language assistance. This work was funded by the National Natural Science Foundation of China (31570381,

31770410) and the State Key Laboratory of Earth Surface Processes
and Resource Ecology (2017-ZY-08).

DATA AVA I L A B I L I T Y S TAT E M E N T
The chloroplast DNA sequence data have been deposited in GenBank with accession nos. KU645591-KU645758. Data relating to MaxEnt files, and microsatellite genotypes have been uploaded to DRYAD: https://datad ryad.org/stash /share /4ihCP pUOWD_zNdlb sTdxD wWD5u GcM7Q m0sI6 LcW-5aE.