Plant Biology ISSN 1435-8603

RESEARCH PAPER

Molecular data and ecological niche modelling reveal the phylogeographic pattern of Cotinus coggygria (Anacardiaceae) in China’s warm-temperate zone W. Wang†, C. Y. Tian†, Y. H. Li & Y. Li College of Forestry, Henan Agricultural University, Zhengzhou, China

Keywords China’s warm-temperate zone; cpDNA; ecological niche model; yellow River. Correspondence Y. Li, Henan Agricultural University, Zhengzhou 450002, China. E-mail: [email protected] † These authors contributed equally to this work.

Editor T. Peeters Received: 5 September 2013; Accepted: 20 December 2013

ABSTRACT The phylogeography of common and widespread species helps to elucidate the history of local flora and vegetation. In this study, we selected Cotinus coggygria, a species widely distributed in China’s warm-temperate zone. One chloroplast DNA (cpDNA) region and ecological niche modelling were used to examine the phylogeographic pattern of C. coggygria. The cpDNA data revealed two phylogeographic groups (Southern and Northern) corresponding to the geographic regions. Divergence time analyses revealed that divergence of the two groups occurred at approximately 147,000 years before the present (BP), which coincided with the formation of the downstream area of the Yellow River, indicating that the Yellow River was a weak phylogeographic divide for C. coggygria. The molecular data and ecological niche modelling also indicated that C. coggyria did not experience population expansion after glaciations. This study thus supports the fact that Pleistocene glacial cycles only slightly affected C. coggygria, which survived in situ and occupied multiple localised glacial refugia during glaciations. This finding is contrary to the hypothesis of large-scale range habitat contraction and retreat into a few main refugia.

doi:10.1111/plb.12157

INTRODUCTION China’s warm-temperate zone generally refers to the area between 32°30′–42°30′ N and 103°30′–124°10′ E, which comprises an area of approximately 7 9 105 km2 (Shangguan et al. 2009). The typical vegetation in this region is deciduous broadleaved forest (DBLF; Gao et al. 2001). This warmtemperate zone has not been strongly influenced by massive Quaternary glaciers since the Tertiary Period. Thus, the majority of deciduous broadleaved species among the vegetation are preserved in this region, in contrast to the large-scale extinction of broadleaved species in Europe and North America (Zhu 1997). Although no massive ice sheet directly covered China’s warm-temperate zone during the Quaternary glacial period, the climate in this region became cooler and drier. Thus, both the evolution and distribution of the warm-temperate DBLF are still affected by climate fluctuations. Phylogeographic analyses are now widely used to reconstruct the history of species using molecular data and ecological niche modelling (Knowles et al. 2007; Qi et al. 2012). However, our knowledge on phylogeographic histories of organisms occurring in China’s warm-temperate zone, as well as causal correlations with climatic fluctuations, has been limited due to finite phylogeographic studies in this region, particularly for plants such as Cercidiphyllum (Qi et al. 2012), Platycarya strobilacea (Chen et al. 2012), Pteroceltis tatarinowii (Li et al. 2012), Bupleurum longiradiatum (Zhao et al. 2013) and Forsythia suspensa (Fu et al., unpublished). Interestingly, these studies

have shown remarkable differences in the phylogeographic history of the species involved, and two hypotheses may be inferred from the findings. The first hypothesis is that climatic changes during the Pleistocene deeply impacted China’s warmtemperate zone, which may have caused plants in this region to retreat into a few main refugia during glaciations. The second hypothesis is that Pleistocene glacial cycles only slightly affected these species, which survived in situ and caused the plants to occupy multiple localised glacial refugia during glaciations. Unlike the biodiversity hotspots in China (e.g. Qinghai–Tibetan Plateau, Himalaya–Hengduan Mountains and subtropical China), China’s warm-temperate zone harbours relatively low biodiversity. However, this region has the highest population density in the country and has well developed agriculture. Dramatic climate change is adversely affecting the global ecological system and agricultural cultivation. Thus, increased attention must be given to China’s warm-temperate zone and to the influences of Quaternary climate change in this area based on more phylogeographic studies. Cotinus coggygria Scop. (Anacardiaceae) is a deciduous shrub or small tree endemic to China’s warm-temperate zone, found in Hebei, Henan, Hubei, Shaanxi, Shandong and Shanxi Provinces, at elevations of 330–2400 m a.s.l. As a typical component of DBLF, C. coggygria is an ideal candidate for investigating the influences of Quaternary climate change in China’s warm-temperate zone. In this study, one chloroplast DNA (cpDNA) region and ecological niche models were used to examine the phylogeographic pattern of C. coggygria. This study aims to elucidate

Plant Biology © 2014 German Botanical Society and The Royal Botanical Society of the Netherlands

1

Phylogeography of Cotinus coggygria

Wang, Tian, Li & Li

(i) the genetic structure of C. coggygria populations in China as revealed by the cpDNA data, and (ii) how this species responded to the climate oscillations of the Quaternary Period. A

C

MATERIAL AND METHODS Population sampling In this study, 105 individual plants from 15 wild populations of C. coggygria were sampled, representing almost the entire range of the species in China (Table 1, Fig. 1). Each population included two to ten individuals that were collected at least 10 m apart to increase the likelihood of sampling inter-individual variation. Healthy leaves were collected and dried in silica gel until total DNA was extracted. DNA isolation, amplification and sequencing Total genomic DNA was isolated from leaves dried in silica gel using the Plant Genomic DNA kit (Tiangen Biotech, Beijing, China) according to the protocol of the manufacturer. After preliminary screening of 12 fragments, we chose one intergenic spacer (trnL-F) for the full survey because it contained the most polymorphic sites. The primers for trnL-F described in Taberlet et al. (1991) were used to amplify the sequences. Polymerase chain reaction (PCR) was then performed in a 30-ll reaction volume containing 30 ng genomic DNA, 0.2 mM of each dNTP, 0.3 lM of each primer, 3 ll Taq buffer and 1 U Taq polymerase (Takara, Dalian, Liaoning, China). The PCR protocols involved initial denaturation for 4 min at 94 °C, followed by 35 cycles of denaturation for 40 s at 94 °C, annealing for 45 s at 50 °C and extension for 90 s at 72 °C, with a final extension step of 8 min at 72 °C. The PCR products were

B Fig. 1. Geographic distribution and genealogical relationships of cpDNA haplotypes recovered from C. coggygria populations in China. A: Distribution ranges of C. coggygria (red lines) in China. B: Sampling localities of 15 populations of C. coggygria and the geographic distribution of 13 haplotypes (H1–H13) detected (for population numbers see Table 1), black dot A to B represented the downstream area of the Yellow River. C: Statistical parsimony network of haplotypes H1–H13. The area of circles corresponds to the frequency of each haplotype. Each solid line represents one mutational step interconnecting two haplotypes, for which parsimony is supported at the 95% level.

purified with an E.Z.N.A. Gel Extraction Kit (Omega Bio-Tek, Winooski, VT, USA) and then sequenced on an ABI 3730 DNA Sequence Analyzer at BGI (Beijing, China).

Table 1. Details of population locations, sample size and cpDNA variation of Cotinus coggygria sampled in China. p: nucleotide diversity, h: haplotype diversity, Nph: number of private haplotypes. population no. and code Southern 1. HBWD 2. HNLJ 3. SXLJ 4. HNSM 5. SXWL 6. HNJL 7. SDYM regional mean Northern 8. SXTB 9. SXTT 10. SDBD 11. SXHM 12. HBWZ 13. SXLK 14. SXTL 15. HBTG regional mean (not including P9) total

2

locations

Lat.(N)/Long.(E)

Wudang Mt., Hubei Laojun Mt., Henan Laojun Mt., Shaanxi Song Mt., Henan Wulaofeng, Shanxi Jiulian Mt., Henan Yuan Mt., Shandong

32°24′/111°00′ 33°45′/111°38′ 34°20′/110°15′ 34°28′/113°05′ 34°50′/110°35′ 35°35′/113°35′ 36°28′/117°51′

Taibai Mt., Shaanxi Tiantai Mt., Shaanxi Baodugu, Shandong Hua Mt., Shaanxi Wuzhi Mt., Hebei Lingkong Mt., Shanxi Tianlong Mt., Shanxi Tiangui Mt., Hebei

33°57′/107°45′ 34°17′/107°11′ 35°00′/117°42′ 35°33′/110°06′ 36°30′/113°39′ 36°36′/112°05′ 37°42′/112°26′ 38°15′/113°44′

N

haplotype (no. of individuals)

p910

4 9 7 10 8 8 8

H5(2), H8(1), H13(1) H2(1), H4(1), H7(2),H9(2), H10(2), H12(1) H2(1), H3(1), H7 (2), H8(1), H9(1), H10(1) H2(1), H3(1), H6(2), H9(2), H10(3), H12(1) H6(3), H7(5) H2(1), H4(1), H7(1), H9(3), H10(2) H2(1), H7(1), H9(4), H11(2)

8 2 5 8 10 4 10 4

H1(3), H3(1), H6(1), H8(3) H2(2) H1(2), H9(1), H12(2) H1(3), H2(3),H3(1), H6(1) H2(3), H4(2), H9(1), H10(3), H12(1) H1(2), H2(1), H12(1) H2(5), H4(3), H9(2) H3(2), H7(1), H12(1)

105

h

Nph

5.10 4.95 5.45 5.17 3.56 5.03 5.37 4.95

0.833 0.917 0.952 0.889 0.536 0.857 0.750 0.819

2 0 0 0 0 0 1 0.43

5.98 0 5.32 3.85 3.04 3.32 2.51 4.43 4.06 4.86

0.786 0 0.800 0.786 0.844 0.833 0.689 0.833 0.796 0.900

0 0 0 0 0 0 0 0 0 3

3

Plant Biology © 2014 German Botanical Society and The Royal Botanical Society of the Netherlands

Phylogeography of Cotinus coggygria

Wang, Tian, Li & Li

Molecular data analysis

Past and current distribution inferences

Sequences were aligned with CLUSTAL_X version 1.81 (Thompson et al. 1997) and all indels were coded as substitutions following the method of Caicedo & Schaal (2004). Haplotype diversity (h) and nucleotide diversity (p) were calculated for each population (hS and pS) and at the species level (hT and pT) using DNASP version 5.0 (Librado & Rozas 2009). A network of the cpDNA haplotypes was constructed using TCS 1.21 (Clement et al. 2000), with the default parsimony connection limit of 95%, and each indel was treated as a single mutation event. To infer the most likely number of population genetic clusters (K), Bayesian analysis of population structure was performed as implemented in STRUCTURE version 2.3.1 (Pritchard et al. 2000). K ranged from 1 to 10, with ten replicates performed for each K and using a burn-in period of 2 9 105 and 5 9 104 Markov chain Monte Carlo (MCMC). The no-admixture model and independent allele frequencies were chosen for this analysis. The most likely number of clusters was estimated according to the model values (DK) based on the second-order rate of change, with respect to K, of the likelihood function (Evanno et al. 2005). Genetic relationships among 15 populations of C. coggygria, based on their pair-wise genetic distances, Kxy, were examine using the neighbour-joining method of the program MEGA 5.2 (Tamura et al. 2011). To quantify the genetic differentiation partitioned among all populations and regional groups, analyses of molecular variance (AMOVA) were carried out using ARLEQUIN version 3.5 (Excoffier & Lischer 2010), with 1000 random permutations, to test for significance of partitions. Isolation-by-distance (IBD) was assessed with the Mantel permutation procedure as implemented in IBD (Jensen et al. 2005, http://ibdws.sdsu. edu/) through testing the correlation between values of Slatkin’s similarity index, M (Slatkin 1993) against the geographic distance (Km). To relate the divergence of groups defined in STRUCTURE to the historic events involving the species, divergence times between groups was estimated. Net pair-wise divergence per base pair (dA) was calculated using MEGA 5.2 (Tamura et al. 2011) under the Kimura-2 model (Kimura 1980). Divergence time was calculated as T = dA/2l, where l is the rate of nucleotide substitution (Nei & Kumar 2000). An appropriate rate has not been calibrated for the chloroplast substitutions of C. coggygria; thus, the general substitution rates of the plastid sequence of 1.6 9 10 9 substitutions per site per year (Wolfe et al. 1987) was obtained as the approximate evolutionary rate. To infer demographic processes, the null hypothesis of the sudden expansion model in ARLEQUIN was tested through comparing the observed and expected distributions of pairwise sequence differences (i.e. mismatch distributions). The sum of squared deviations (SSD) and raggedness index (RAG) were used to test the goodness-of-fit of the observed mismatch distribution to the model expectation. In addition, tests of selective neutrality, namely Tajima’s D (Tajima 1989) and Fu’s Fs (Fu 1997) values, were used to infer potential population growth and expansion. If the expansion model was not rejected in the mismatch analysis and neutrality test, the relationship s = 2 lt (Rogers & Harpending 1992) was used to estimate the age of expansion (t), where l is the substitution rate for DNA sequences.

Ecological niche modelling was carried out in MAXENT version 3.3.3 (Phillips et al. 2006) to predict suitable past (i.e. Last Glacial Maximum [LGM]) and present climate envelopes for C. coggygria. Information on the geographic distribution of C. coggygria was based on the 29 populations included in this study, in which 14 points were obtained from the Chinese Virtual Herbarium (http://www.cvh.org.cn/cms_cvh) and the remaining 15 points were from the sampling sites. Current bioclimatic variables and LGM data were downloaded from the WorldClim database (http://www.worldclim.org/; Hijmans et al. 2005). LGM data were simulated using the Community Climate System Model (Collins et al. 2006), and 19 environmental parameters were used to model the potential distribution of the species. To test the performance of each model, 20% of the data in each run was randomly selected by Maxent and compared with the model output generated with the remaining data. The area under the receiver operating characteristic curves (area under curve [AUC]) was used to compare model performance (Phillips et al. 2006). RESULTS Genetic diversity and population structure One cpDNA region (trnL-F) was sequenced in C. coggygria from 105 individuals (15 populations). The region showed a length variation from 721 bp to 745 bp and aligned with a consensus length of 752 bp, which contained 14 nucleotide substitutions and four indels. Based on these polymorphisms, 13 haplotypes (H1 to H13) were identified among all the samples surveyed (Table 2). The 13 haplotype sequences were deposited in the GenBank database under accession numbers KF600590 to KF600602. Among the detected, the most widespread were H2 (in ten of 15 populations) and H9 (in eight of 15 populations). The geographic distribution of haplotypes H1 to H13 and their occurrence at each locality are shown in Fig. 1B and Table 1. Among the 15 populations, 14 were polymorphic, whereas only one (P9) exhibited only one haplotype because of its small sample size (i.e. two individuals). The statistical parsimony network of haplotypes H1 to H13 revealed that they were only one to eight mutational steps apart (Fig. 1C). The haplotype and nucleotide diversities of C. coggygria were hT = 0.900 and pT = 4.86 9 10 3, respectively. The highest nucleotide and haplotype diversities were found in P8 and P3, respectively (Table 1). Three private haplotypes of populations were found in C. coggygria, with P1 having the highest number of private haplotypes. In the Bayesian analysis of population structure (Fig. 2), the highest likelihood of the molecular data was obtained when the samples were clustered into two groups (K = 2; Fig. 3), namely, Southern group (P1 to P7) and Northern group (P8 to P15). The population-based NJ network (Fig. 4) confirmed this strong pattern of regional differentiation. Non-hierarchical AMOVA (Table 3) revealed a weak population genetic structure for cpDNA sequence variation at the entire species range scale (ΦST = 0.103, P < 0.05). This differentiation was partitioned between the two groups (ΦCT = 0.084), and 5.32% was explained among populations within each region (ΦSC = 0.058; Table 3). Significant IBD for cpDNA was detected at the entire

Plant Biology © 2014 German Botanical Society and The Royal Botanical Society of the Netherlands

3

Phylogeography of Cotinus coggygria

Wang, Tian, Li & Li

Table 2. Chloroplast DNA sequence polymorphisms detected in trnL-F of Cotinus coggygria, identifying 13 haplotypes (H1–H13). All sequences are relative to the reference haplotype H1. Numbers 1/0 in sequences denote presence/absence of length polymorphism, identified with a superscript letter (a, b, c, d). nucleotide position haplotype H1 H2 H3 H4 H5 H6 H7 H8 H9 H10 H11 H12 H13

91 G A G A A A A A A G G A A

130 A A G A A A A A A A G A A

174 a

1 1a 1a 1a 0 1a 1a 1a 1a 1a 1a 1a 1a

196

211

246

497

517

523

540

552

588

606

638

656

690

701

748

G G G G G T G G G G G G G

A A A A A G A A G A A A G

0 0 0 0 0 0 0 0 0 0 0 1b 0

C C C C A A A C A C C C A

C C C G C A C C C C C C C

G G G G G G A A G G G G G

T T T T T T T G T T T T T

C C C C C C C C T C C C C

0 1c 1c 1c 1c 1c 1c 1c 1c 1c 1c 1c 1c

A A A A A C A A A A A A A

0 0 0 0 0 0 0 0 1d 0 0 0 0

C A A A A A A A A C A A A

T T T T C T T T T T T T T

G G G G G G G G G G C G G

T T T T T T T G T T T T T

a: G; b: ATACGTAATACACTATGT; c: ACAAAC; d: GTATGTG.

Fig. 2. Estimated genetic structure for K = 2 obtained with the program STRUCTURE for 15 populations of C. coggygria based on cpDNA data. Each vertical bar represents an individual and its assignment proportion (not shown) into one of two (coloured) population clusters (K).

genetic distance (dA) between the Southern and the Northern groups was 0.00047, suggesting that divergence between the two groups occurred approximately 147,000 years before the present (BP). Tajima’s D and Fu Fs statistics for deviation from neutrality were not significant for the two subregions and the whole species dataset (all P > 0.05), indicating that no population expansion had occurred for C. coggygria. Similarly, the mismatch analysis for the whole species also showed no demographic expansions for significant SSD and RAG values. Unlike the two neutrality tests, the two subregions showed demographic expansions that were not significant for SSD and RAG values (Table 4).

A

B

Distribution inference by ecological niche modelling

Fig. 3. The uppermost hierarchical level of genetic structure in the cpDNA dataset of C. coggygria determined using STRUCTURE. (A) Values of log probability L(K) for the data as a function of K. (B) Values of △K based on the second-order rate of change, with respect to K, of the likelihood function.

species range scale (r = 0.234, P < 0.05), whereas no significant IBD was present in the two subregions (Southern: r = 0.234, P = 0.179; Northern: r = 0.085, P = 0.389). The 4

The inferred current and past (LGM) distributions of C. coggygria are illustrated in Fig. 5. The AUC values based on both training and test presence data for the present and LGM were higher than expected (0.995 and 0.993, 0.996 and 0.995, respectively), demonstrating good model performance. Compared with the two simulated distributions, the LGM distribution model predicted that the species did not contract into a few main refugia or migrate long distances into southern refugia. However, the species retreated at the eastern edge (Shandong Province) during the LGM (Fig. 5B). After glaciation, C. coggygria also did not experience significant habitat expansion (Fig. 5A).

Plant Biology © 2014 German Botanical Society and The Royal Botanical Society of the Netherlands

Phylogeography of Cotinus coggygria

Wang, Tian, Li & Li

Fig. 4. Neighbour-joining network illustrating the genetic relationships among 15 populations of C. coggygria based on their pair-wise genetic distances Kxy.

DISCUSSION Genetic diversity and spatial population genetic structure In our study, 13 haplotypes were identified in the C. coggygria populations analysed. The cpDNA variance among populations showed high levels of haplotype diversity (hT = 0.900) and nucleotide diversity (pT = 4.86 9 10 3). These results significantly exceeded the average value of 0.67 for the 170 plant species documented (Petit et al. 2005). Genetic diversity in plant species is usually affected by breeding system, distribution range and population size (Hamrick & Godt 1989, 1996; Hamrick et al. 1992; Nybom 2004). As a com-

Table 3. Non-hierarchical and hierarchical AMOVAs for cpDNA variation surveyed in populations of Cotinus coggygria in China. source of variation

df

non-hierarchical AMOVAs total 14 Southern 6 Northern 7 hierarchical AMOVAs among two groups 1 among populations 13 within populations 90

% total variance

Φ-statistic

P-value

10.27 7.01 4.63

ΦST = 0.103 ΦST = 0.070 ΦST = 0.046