Taxonomic status of Populuswulianensis and P.ningshanica (Salicaceae)

Abstract Species delimitation in the genus Populus is particularly challenging due to high levels of intraspecific polymorphism as well as frequent interspecific hybridisation and introgression. In this study, we aimed to examine the taxonomic status of Populusningshanica and P.wulianensis using an integrative taxonomy that considers multiple operational criteria. We carried out morphometric analyses of leaf traits and genetic examinations (including sequence variations at five barcoding DNAs and polymorphisms at 14 nuclear microsatellite SSR primers) at the population level between them and two closely related species P.adenopoda and P.davidiana. Results suggest that P.wulianensis belongs to the polymorphic species, P.adenopoda and should be considered as a synonym of the latter. P.ningshanica may have arisen as a result on the hybridisation between P.adenopoda and P.davidiana and therefore should be treated as P.×ningshanica. This study highlights the importance of the integrated evidence in taxonomic decisions of the disputed species.


Introduction
Species delimitation is essential to conserve and assess biodiversity (Agapow et al. 2004). Any incorrect species recognition may result in serious after-effects in related studies, for example, by an increase in species conservation (Wiens 2007) and un-der-or over-estimation of biodiversity (Douady 2007). Therefore, in addition to morphological traits, significant efforts have been made to delimit species based on DNA sequence variation (Wiens and Penkrot 2002;Sites and Marshall 2003;Kress et al. 2005;Bond and Stockman 2008;Fujita et al. 2012;Hendrixson et al. 2013) or other genetic polymorphisms that can assess gene flow and identify interspecific hybrids according to the biological species concept (Pérez-Losada et al. 2005). These molecular markers have been used to differentiate species, hybrids and even clones in the genus Populus (Salicaceae) (Hamzeh and Dayanandan 2004;Cervera et al. 2005;Hamzeh et al. 2006;Fladung and Buschbom 2009;Schroeder et al. 2012;Feng et al. 2013;Wan et al. 2013). Poplars are widely distributed in the Northern Hemisphere with an important ecological role in natural and artificial forests in both boreal and temperate regions (Dickmann et al. 2001). However, due to high levels of morphological variation and extensive inter-specific hybridisation, species delimitation within the genus is highly contentious (Eckenwalder 1996;Dickmann and Kuzovkina 2008). The number of the proposed species ranges from 22 to 85, plus hundreds of hybrids, varieties and cultivars (Dickmann and Stuart 1983;Fang et al. 1999). Numerous described species were doubted as being hybrids of the other independently evolving lineages (good species) or intra-specific variations of the polymorphic species. However, these ambiguous species have not been well examined.
In this study, we aimed to determine the taxonomic status of two species described from China: P. wulianensis S.B.Liang & X.W.Li and P. ningshanica C. Wang & Tung (Fang et al. 1999) based on morphometric analyses and genetic examinations at the population level as recently suggested for an integrated species delimitation (Liu 2016). P. wulianensis is restricted to eastern Shandong while P. ningshanica is distributed in southern Shaanxi and Northwest Hubei. Both are morphologically similar to P. davidiana Dode and P. adenopoda Maxim. of sect. Populus with widespread distributions in northern or middle to southern China. The key traits for their diagnosis are mainly based on leaf characters: blade and apex shape and margin incision (Fang et al. 1999). We firstly conducted morphometric analyses of leaf traits for representative populations of all four species. Then we examined genetic delimitations between them based on evidence from sequence variation of internal transcribed spacer (ITS) and four chloroplast DNA (cpDNA) and genetic polymorphisms from nuclear microsatellite loci (nSSR).

Sample collection
We sampled 163 individuals from 17 populations of four species (Table 1), including all recorded natural populations of both P. ningshanica and P. wulianensis. All individual trees were chosen with typical morphological leaf traits (Fang et al. 1999). Each tree was set apart by at least 50m in each population. Except for collecting specimens (SZ, herbarium of Sichuan University, Chengdu, China) for geometric morphometric analyses, we further selected healthy and fresh leaves from each tree and dried them immediately in silica gel for DNA extraction. We also used an Etrex GIS monitor (Garmin, Taiwan) to record latitude, longitude and altitude of each sampled population (Table 1; Fig. 1A).

Geometric morphometrics
Although we failed to find type specimens of P. ningshanica and P. wulianensis, we included the newly collected specimens from their type localities. A Canon 60D digital camera was used to photograph typical leaves of all specimens. We transformed every image into a vector diagram using TpsUtil version 1.64 (Rohlf 2013). Thirty-two homologous landmarks were assigned in order to quantify leaf blades shape in all specimens. Landmark positions of leaves included base, tip and margin. All landmarks were digitised for each individual using the software TpsDig version 2.22 (Rohlf 2015). We created a combined data file including all specimens. We implemented morphometrics analyses in MorphoJ version 1.01b (Klingenberg 2011), within which a principal component analysis of morphological variations was conducted and plotted.

Genetic analyses
We isolated the total genomic DNA from leaves of each individual, based on the hexadecyltrimethyl ammonium bromide (CTAB) method (Doyle 1987). We used a total of 14 SSRs primers (Suppl. material 1: Table S1) developed previously, based on the genome sequences of Populus euphratica and P. trichocarpa (Ma et al. 2013;Jiang et al. 2016) to genotype our samples. The PCRs were performed in a volume of 25 ml, which contained: 50-100 ng diluted genomic DNA, 0.5 mM of each dNTP, 0.5 µl of each primer, 2.5 µl 10 × Taq buffer and 0.5 units of Taq polymerase (Vazyme Biotech, Nanjing, China). The PCR programme used was: initially a single cycle at 95 °C for 5 min, followed by 36 cycles at 95 °C for 45s, 55 °C for 40 s and 72 °C for 80 s, with a final extension at 72 °C for 10 min. The PCR products at each locus were analysed on an ABI 3830xl DNA analyser (Applied Biosystems, Inc., Foster City CA) at Tsingke Biological Technology (Beijing, China). We used STRUCTURE version 2.3.4 (Falush et al. 2003) that allows a Bayesian hybrid mixture computation to identify genetic compositions of all sampled trees. We pre-assigned a number of genetic clusters (K) ranging from 1 to 10. All runs involved 1,000,000 Markov chain Monte Carlo repetitions after a burn-in period of 500,000 iterations. We used the long burn-in and run lengths as well as 10 replicates to ensure the reproducibility of STRUCTURE results (Gilbert et al. 2012). We estimated the posterior probability of K and Delta K (ΔK), the rate of change of Ln P (K) between successive K values (Evanno et al. 2005). We determined the most likely number of clusters. We also sequenced internal transcribed spacer (ITS) and four chloroplast DNA (cpDNA) fragments: matK, trnH-psbA, trnG-psbK and psbK-psbI for three to five individuals from each sampled population of four species used for nSSR genotyping. In addition, one individual of P. euphratica was sequenced as the outgroup. Primers, PCRs and sequencing followed Feng et al. (2013) (Suppl. material 1: Table S2). Sequences for each fragment were aligned and sequences from four cpDNAs were connected using MEGA 7.0 (Kumar et al. 2016). We constructed unrooted neighbour-joining (NJ) trees for both ITS and cpDNAs datasets by MEGA 7.0 (Kumar et al. 2016) respectively, using pairwise deletion and the P-distance model. Bootstrap values were estimated with 1000 random addition sequence replicates.

PCA analyses of geometric morphometric data
Geometric morphometric analyses of leaf traits yielded 30 principal components (PC), which accounted for all leaf variations. PC1 to PC3 were the only PCs that individually represented >5% of the variance (PC1=57.05%; PC2=12.69%; PC3=7.68%) and they together represented 77.43% of the variance. All other PCs accounted for <5% of the variance individually. The greatest amount of shape variance is observed across PC1 and PC2 (Fig. 1B). Across these two axes, individuals of P. davidiana and P. adenopoda were treated as a clear division, whereas individuals of P. wulianensis and P. ningshanica are clustered into one subgroup of the P. adenopoda group. All other PCs showed similar relationships.

Clustering analyses based on the SSR polymorphisms
We genotyped 14 nuclear SSR loci for 163 sampled individuals of four species. Using the method originally described by Pritchard et al. (Pritchard et al. 2000) and also the ΔK approach described by Evanno et al. (Evanno et al. 2005), we found the most likely number of Bayesian clusters was two (K = 2) ( Fig. 2A). When K = 2, individuals from P. davidiana clustered into one group and those from P. adenopoda into the other. Within each group, some samples indicated the weak genetic introgression from the other. All sampled individuals of both P. wulianensis and P. ningshanica were assigned to the group represented by P. adenopoda (Fig. 2B). However, approximately 10% of the genetic composition of P. ningshanica derived from the cluster represented by P. davidiana, while more than 90% was from P. adenopoda. Similar results were obtained based on PCA analyses of genetic polymorphisms and that two groups were identified to be, respectively, represented by P. davidiana and the other three (Fig. 2C).  We have combined sequences of four cpDNAs for each individual into one cp-DNA sequence. We aligned the cpDNA sequences of all individuals and identified 2, 1, 1 and 2 sequences for P. davidiana, P. adenopoda, P. wulianensis and P. ningshanica, respectively. The total length of the aligned cpDNA sequence was 1866 bp with 9 vari-able sites amongst different sequences from four species (Fig. 3B). NJ clustering of all different cpDNA sequences from four species similarly identified two tentative groups: one comprised of P. davidiana and P. ningshanica, while the other included those from P. adenopoda, P. wulianensis and P. ningshanica. We identified 1, 2, 1 and 2 different ITS sequences for the sampled individuals for P. davidiana, P. adenopoda, P. wulianensis and P. ningshanica. We aligned these ITS sequences from four species, which were 552 bps long with 1 variable site amongst all the different sequences from four species (Suppl. material 1: Table S3; Fig. 3A). NJ analyses of the ITS dataset identified two tentative groups: one comprised 4 sequences from P. adenopoda and P. ningshanica while the other, all four species.

Discussion
Statistical analyses based on geometric morphometric measurements are highly successful at separating similar species (Villemant et al. 2007;Francuski et al. 2009), even when the individual character shows the overlapped variations between them (Lumley and Sperling 2010;Buck et al. 2012). Especially, geometric morphometrics could differentiate the overall changes in the gross morphology (Rohlf and Marcus 1993). Poplar leaves are ideal for geometric morphometric analyses, as they are two-dimensional, easily imaged and the venation provides many points that are clearly homologous and straightforward to landmark accurately. In addition, flower traits are highly static across the genus without variations and leaf characters are therefore used to classify different species (Dickmann and Stuart 1983;Eckenwalder 1996;Fang et al. 1999;Dickmann and Kuzovkina 2008). We tried to classify four popular species based on geometric morphometric analyses of leaf traits. Our results obviously suggested that P. davidiana and P. adenopoda differed distinctly from each other. P. wulianensis and P. ningshanica could not be distinguished from each other and they together clustered into one subgroup, which obviously belonged to the P. adenopoda group (Fig. 1B). Therefore, this statistical clustering indicated that both P. wulianensis and P. ningshanica may belong to the polymorphic P. adenopoda.
Genetic evidence, based on nuclear SSR loci, similarly recognised the distinct species boundary between P. davidiana and P. adenopoda ( Fig. 2A, B). However, all sampled individuals of P. wulianensis belong to the P. adenopoda group without distinct introgression from P. davidiana. All sampled individuals of P. ningshanica shared similar genetic compositions, together belonging to the P. adenopoda group but with obvious genetic introgressions from the P. davidiana group. These individuals comprise the obvious backcrosses from P. adenopoda. Similarly, sequence variations from five DNAs (ITS, matK, trnH-psbA, trnG-psbK and psbK-psbI) seem to support these inferences. The connected sequences of four cpDNAs distinguished P. davidiana and P. adenopoda while all P. wulianensis individuals shared the same cpDNA sequences with P. adenopoda. We found two types of cpDNA sequences in P. ningshanica (Fig. 3B), clustering respectively with those from P. davidiana and P. adenopoda, which further suggested the hybrid origin of P. ningshanica. However, the initial hybrids must have repeatedly backcrossed with P. adenopoda, which resulted in the high genetic similarity of the sampled individuals of P. ningshanica to P. adenopoda but with introgression with P. davidiana (Fig. 3B). The interspecific hybrids in the genus Populus could be F1, F2 to multiple generation backcrossing hybrids (Braatne et al. 1992;Bradshaw et al. 2000;Feng et al. 2013;Jiang et al. 2016). We failed to find stable ITS differences between P. davidiana and P. adenopoda. It is highly probably that the gene flow, mediated by interspecific hybrids, had caused the concerted evolutions and indistinct differences in the ITS sequence variations (Feng et al. 2013;Jiang et al. 2016).
Overall, multiple lines (Figs 1B, 2B, C;) of evidence suggested that P. wulianensis was described based on the intraspecific variations of the polymorphic P. adenopoda and individuals ascribed to P. ningshanica are, in fact, hybrids between P. adenopoda and P. davidiana with the repeated backcrosses to the former. Both taxa should be treated accordingly in the taxonomic revision of the genus Populus. Additional