The Philippines is an archipelagic country of around 7000 islands that supports immense biological diversity (Myers et al. 2000), including more than 9430 native species of vascular plants of which 50.1% are endemic (Pelser et al. 2011). In the past few centuries, the Philippines has lost 80.6% of its forest cover (Laurance 2007), resulting in the reduction of the size and genetic connectivity of populations of forest inhabiting species (Young et al. 1996; Ricketts 2001; Harris and Reed 2002). Such changes are expected to result in a loss of genetic diversity and, consequently, of the evolutionary potential of species, compromising their long-term persistence (Ellstrand and Elam 1993; Jamieson et al. 2008). Estimates of genetic diversity and genetic connectivity can inform conservation management aimed at reducing genetic diversity loss in isolated forest habitats (Marshall et al. 1999; Schwartz et al. 2007; Jamieson et al. 2008; Frankham et al. 2010; Luque et al. 2012). However, to our knowledge, detailed information about patterns of genetic diversity and connectivity of plant populations in fragmented forest habitats in the Philippines (e.g. Pelser et al. 2017, 2018) is scarce.

Tetrastigma loheri Gagnep. is a species of woody vines in the grape family (Vitaceae) that is native to Borneo and the Philippines (Pelser et al. 2011). It is potentially a suitable species for studying patterns of genetic diversity and connectivity among fragmented forest ecosystems in various parts of the Philippines because T. loheri is a widespread and commonly encountered forest species (Pelser et al. 2016) and can therefore provide data from a sufficiently high number of individuals to yield statistically well-supported genetic patterns. However, the species delimitation of T. loheri is presently unclear and this taxon might be part of a species complex in the Philippines. The T. loheri s. l. complex (Pelser et al. 2016) is composed of T. loheri and four previously and currently recognized Tetrastigma species in the Philippines that display morphological similarities: Tetrastigma diepenhorstii (Miq.) Latiff, Tetrastigma philippinense Merr., Tetrastigma stenophyllum Merr., and Tetrastigma trifoliolatum Merr. (Miquel 1861; Gagnepain 1910; Merrill 1912, 1914, 1916). These five species all have glabrous and coriaceous leaves that almost always have three leaflets (and never more than that) in combination with petiolules that are longer on the terminal leaflet than on the lateral leaflets, a leaflet margin that is generally toothed, and a rather compact inflorescence when compared to other Tetrastigma species (Miquel 1861; Gagnepain 1910; Merrill 1912, 1914, 1916).

Tetrastigma loheri only differs from the other four species in the complex in relatively subtle morphological differences (Table 1). Both T. diepenhorstii and T. trifoliolatum are recorded as having pubescent inflorescences, whereas those of T. loheri are reported to be glabrous (Miquel 1861; Gagnepain 1910; Merrill 1914). In addition, T. diepenhorstii has elliptic to ovate (vs. lanceolate to oblong) leaflets and longer petioles (6.5–28 cm vs. up to 4 cm) than T. loheri (Miquel 1861; Gagnepain 1910; Wan Zakaria et al. 2016, 2017). As described in its protologue (Merrill 1914), T. trifoliolatum has larger leaflets (14–20 × 7–10 cm vs. 4.5–12 × 2.5–5 cm), longer petioles (ca. 12 cm vs. up to 4 cm), and longer petiolules of the lateral leaflets (3–4 cm vs. up to 0.5–1 cm) than T. loheri (Gagnepain 1910). Merrill (1912) described T. philippinense as having pubescent petals whereas those of T. loheri are glabrous according to Gagnepain (1910) and Merrill (1912). However, Merrill (1916) later synonymized T. philippinense under T. loheri. Although he did this without providing a detailed argument, this suggests that he considered their morphological differences too minor to merit taxonomic distinction at the level of species. Tetrastigma stenophyllum reportedly has narrower leaves than T. loheri (Merrill 1916). Merrill (1916) also mentioned that T. stenophyllum differs from T. loheri in having pubescent ovaries and stellate stigma lobes but did not explicitly mention the presence or type of indumentum of the ovaries of T. loheri nor the details of the morphology of its stigma lobes.

Table 1.

Overview of morphological differences between Philippine members of the Tetrastigma loheri s. l. complex and their distribution as reported in the literature (Miquel 1861; Gagnepain 1910; Merrill 1912, 1914, 1916; Latiff 2001; Pelser et al. 2011; Wen et al. 2013; Wan Zakaria et al. 2016, 2017). Character states not reported in the literature are marked —.

As a whole, the T. loheri s. l. complex displays substantial morphological variation, particularly in leaf morphology, such as in the size, shape, and length/width ratio of the leaves and leaflets, the shape of the leaflet margin, venation patterns of leaflets, and in the length of petioles and petiolules (Pelser et al. 2016). To further complicate species delimitation within this complex, some plants cannot be easily accommodated in any of its five species, because they show a conflicting combination of species-level diagnostic character states. In order to be able to use T. loheri for conservation genetic studies in the Philippines, it is important that the species delimitation of the T. loheri s. l. complex is resolved, because the population genetics methods that are required for these studies can only provide meaningful results if the individuals included in the datasets are conspecific.

The reproductive morphology of plants can provide useful characters for species delimitation and identification (e.g. Brown and Gilmartin 1984; Sytsma 1988; Smith et al. 2008; Gahagen et al. 2016). However, the reproductive characters of T. loheri s. l. are poorly known. This is in part because most herbarium specimens of Tetrastigma that are available for study lack flowers or fruits (Pelser et al. 2016). This is perhaps because plants might flower and fruit infrequently and are therefore rarely encountered in the field, or because these reproductive parts might be difficult to see if they are primarily produced on parts of the vines that grow in the canopy of dense forest vegetation and are therefore seldom collected (Pelser et al. 2016). Furthermore, Tetrastigma species are dioecious and staminate and pistillate flowers are therefore present on different individuals. This complicates species delimitation studies that use morphological approaches, because it can be difficult to determine if staminate and pistillate plants belong to the same species if differences in vegetative morphology between species are wanting. Because of the challenges related to using reproductive characters to aid species delimitation in Tetrastigma, vegetative characters are a more accessible source of morphological information for taxonomic studies in T. loheri s. l. Among vegetative characters, leaf shape can be taxonomically informative (e.g. Meade and Parnell 2003; Plotze et al. 2005; Du et al. 2007; Cope et al. 2012), but morphometric studies of leaf shape have not previously been performed for T. loheri s. l.

In addition to morphological characters, patterns of evolutionary relationships can be used to determine the number of species in species complexes and their delimitation (e.g. Costea et al. 2008; Jargeat et al. 2010; Salicini et al. 2011; Gutiérrez-Gutiérrez et al. 2012; Li and Yan 2013). However, the phylogenetic relationships within the T. loheri s. l. complex are still unclear. To date, T. philippinense and T. stenophyllum have not been included in phylogenetic analyses and previous studies resulted in conflicting patterns of relationships among T. diepenhorstii, T. loheri, and T. trifoliolatum (Chen et al. 2011; Wen et al. 2013; Pelser et al. 2016; Habib et al. 2017).

This study aimed to determine if the T. loheri s. l. complex in the Philippines is composed of more than one species by using an integrative taxonomic approach in which different lines of evidence are used to support species recognition (Sukumaran and Knowles 2017; Luo et al. 2018). First, we expanded the existing molecular phylogenetic dataset of Pelser et al. (2016) by adding the first nuclear and plastid DNA sequence data from specimens of T. philippinense and T. stenophyllum, more specimens from the other taxa in the T. loheri s. l. complex, and by sequencing an additional DNA region (the external transcribed spacer; ETS) that, despite its proven ability to resolve species-level relationships in other families (e.g. Baldwin and Markos 1998; Bena et al. 1998; Linder et al. 2000; Markos and Baldwin 2001; Sallares and Brown 2004; Pelser et al. 2010), has not previously been used for phylogenetic studies in Vitaceae. These data were subsequently used for species delimitation using generalized mixed Yule coalescent (GMYC; Pons et al. 2006; Monaghan et al. 2009; Fujisawa and Barraclough 2013) and Poisson tree processes (PTP; Zhang et al. 2013) analyses. These methods aim to identify putative species by discriminating between-species coalescence from within-species coalescence using information from branching rates. Because it is not uncommon that different species delimitation analyses result in incongruent results (Carstens et al. 2013), it is important that putative species delimited by GMYC and PTP are validated with other lines of evidence (Talavera et al. 2013; Zhang et al. 2013). Therefore, we used a conservative, multifaceted approach to species delimitation in which putative GMYC- and PTP-delimited species are only accepted as species if they are also diagnosably morphologically distinct. Thus, we adopted a unified species concept (De Queiroz 2007) for our study, in which species are defined as (segments of) separately evolving metapopulation lineages (primary defining property) that are both phylogenetically and morphologically distinct (secondary defining properties).

Materials and Methods

Phylogenetic Sampling Strategy—This study aimed to resolve the species delimitation of the T. loheri s. l. complex in the Philippines by expanding the dataset of Pelser et al. (2016), which is composed of DNA sequence data of the internal transcribed spacer (ITS) and five plastid regions (atpB-rbcL, psbA-trnH, rps16, trnL, and trnL-F). We primarily focused on adding DNA sequence data of the ITS and ETS regions for additional specimens of the T. loheri s. l. complex, because previous results (Pelser et al. 2016) and preliminary data of the present study showed that these nuclear regions are substantially more informative for providing phylogenetic resolution within the T. loheri s. l. complex than the plastid DNA regions used by Pelser et al. (2016). A smaller plastid dataset consisting of a subset of specimens from the clades that were identified in preliminary phylogenetic ITS and ETS analyses was compiled with the aim of providing better resolution and support for the relationships between these clades.

In total, DNA sequence data of 104 specimens of which 64 are members of the Philippine T. loheri s. l. complex (Appendix 1) from 33 populations throughout the Philippines were used in this study. In addition, DNA sequences of 33 specimens of other species of Tetrastigma (including representatives of seven Philippine species; Appendix 1) were selected as phylogenetic context for determining the relationships among specimens representing T. diepenhorstii, T. loheri, T. philippinense, T. stenophyllum, and T. trifoliolatum. Sequences of representatives of a few other genera of Vitaceae were included in the datasets to root the Tetrastigma clade. To increase sampling coverage for T. loheri s. l. in the Philippines, 21 of the 104 specimens were newly collected for this study in regions that were not previously sampled. Type specimens of T. loheri, T. philippinense, T. trifoliolatum, and T. stenophyllum were not available for this study and specimens from their type localities that morphologically conform to these species as circumscribed in their protologues (Gagnepain 1910; Merrill 1912, 1914, 1916) were therefore collected to inform taxonomic conclusions instead.

DNA Extraction and Sequencing—During the fieldwork for this study, fresh leaf samples from each collected specimen were preserved in silica gel. Voucher specimens were deposited in CAHUP, CANU, CEBU, PNH, and PUH (acronyms following Thiers 2020) (Appendix 1). Approximately 6 mg of dried leaf tissue was pulverized with two metal beads for 2 minutes in an Oscillating Mill MM400 (Retsch GmbH, Haan, Germany). Total genomic DNA was subsequently extracted using a QIAGEN DNeasy Plant mini kit (Germantown, Maryland) following the manufacturer's protocol.

The ITS and the plastid regions were PCR amplified and sequenced with the primers that were used by Pelser et al. (2016) (Table 2). The ETS region was sequenced using primers newly developed for the present study (Table 2). The ETS region was amplified in a 15 µl reaction containing 1 µl of DNA, 2 × Kappa plant PCR buffer, 4.5 pmol each of the forward and reverse primers, and 0.3 U Kappa3G Taq polymerase. The ITS region was amplified in a 15 µl reaction containing 1 µl of DNA, 5 × Promega Go Taq Flexi Green buffer, 6 pmol each of the forward and reverse primers, 3 nmol of each dNTPs, 37.5 nmol of MgCl₂, and 0.12 µl of Go Taq polymerase. The atpB-rbcL region was amplified in a 15 µl reaction containing 1 µl of DNA, 5 × Promega GoTaq Flexi Green, 6 pmol each of the forward and reverse primers, 1.875 nmol of each dNTPs, 56.25 nmol of MgCl₂, and 0.18 µl of GoTaq polymerase. The psbA-trnH region was amplified in a 10 µl reaction containing 1 µl of DNA, 2 × Kappa Ready mix, and 4 pmol each of the forward and reverse primers. The rps16 and trnL-F regions were amplified in a 25 µl reaction containing 1 µl of DNA, 2 × Kappa plant PCR buffer, 7.5 pmol each of the forward and reverse primers, and 0.5 U of Kappa 3G Plant DNA polymerase. The trnL region was amplified in a 15 µl reaction containing 1 µl of DNA, 2 × Kappa Ready mix, and 6 pmol each of the forward and reverse primers.

For ETS, the PCR conditions were: initial denaturation at 94°C for 2 minutes, followed by 35 cycles of denaturation at 94°C for 30 sec, annealing at 51°C for 30 sec, and extension at 72°C for 1 minute. A final extension of 72°C for 10 minutes was used. For ITS, the PCR conditions were: 94°C initial denaturation for 3 minutes followed by 94°C for 1 minute, 52°C for 1 minute, and 72°C for 1 minute for 35 cycles, and a final extension of 72°C for 3 minutes. For atpB-rbcL, psbA-trnH, and trnL, the PCR conditions were: initial denaturation at 97°C for 3 minutes, followed by 37 cycles of: 20 sec at 94°C, 30 sec at 50°C, and 40 sec at 72°C; final extension at 72°C for 5 minutes. The PCR conditions for rps16 and trnL-F were: initial denaturation at 95°C for 10 minutes, followed by 40 cycles of 20 sec at 95°C, 15 sec at 50°C, and 30 sec at 72°C; final extension for 30 sec at 72°C.

The PCR products were purified using the Promega Wizard® SV Gel and PCR Clean-Up System (Madison, Wisconsin) following the manufacturer's protocol. Cycle sequencing was carried out with BigDye® Terminator v3.1 (ThermoFisher, Auckland, New Zealand) using the same primers as for the PCR amplification. The sequenced samples were run on an ABI 3130xL Genetic Analyzer at the University of the Canterbury. Geneious v. 6.1.8 (Biomatters, Auckland, New Zealand) was used for editing the sequencing trace files.

DNA Sequence Alignment and Phylogeny Reconstruction—A total of 372 DNA sequence accessions (Appendix 1) were used to compile five datasets for phylogenetic analysis: ITS (84 specimens), ETS (63 specimens), the combined ITS-ETS data (i.e. the nuclear dataset; 85 specimens), the combined plastid data (atpB-rbcL, psbA-trnH, rps16, trnL, and trnL-F; 85 specimens), and the combined nuclear and plastid dataset (104 specimens). Of these 372 accessions, 108 were newly generated for this study and 264 were obtained from GenBank. Sequences of Ampelocissus, Causonis, and Cayratia (Vitaceae) were used to root the Tetrastigma phylogenies.

DNA sequences were aligned using the Geneious Alignment method in Geneious. Insertions and deletions (indels) were subsequently coded as binary characters with Gapcode.py v. 2.1 (distributed by Richard Ree, Field Museum, Chicago, Illinois, USA,  http://www.bioinformatics.org/∼rick/software.html), which uses the simple indel coding method of Simmons and Ochoterena (2000).

jModelTest2 2.1.6 (Guindon and Gascuel 2003; Darriba et al. 2012) was used on the Cyberinfrastructure for Phylogenetic Research Science Gateway v. 3.3 portal (CIPRES;  https://www.phylo.org; Miller et al. 2010) to select a model of DNA substitution for each dataset. The number of substitution schemes to be considered was set at three. Using the Akaike information criterion (AIC), GTR + G was selected as the model for all DNA regions except psbA-trnH, for which the HKY + G model was used. The Markov k model (Lewis 2001) was used for indel characters.

Bayesian inference (BI) analysis as implemented in MrBayes 3.2.6 (Ronquist et al. 2012) was used for phylogeny reconstruction on the CIPRES cluster. The Markov Chain Monte Carlo (MCMC) analyses were set for 100,000,000 generations and were performed using four chains with a temperature setting of 0.001 and one tree saved per 100 generations. The analyses were set to stop when the average deviation of the split frequencies between both simultaneous analyses reached a value less than 0.01, suggesting convergence. The first 25% of the sampled trees were discarded as burn-in. The consensus trees resulting from the BI analyses were visualised using FigTree v. 1.4.2 ( http://tree.bio.ed.ac.uk/software/figtree/). Clades with posterior probabilities (pp) ≥ 0.95 are considered strongly supported.

Ultrametric Tree Generation—The combined nuclear and plastid dataset resulted in the Bayesian consensus tree with the highest resolution (see Results) and was therefore used in BEAST v. 1.10.4 (Drummond and Rambaut 2007; Suchard and Rambaut 2009) to generate an ultrametric tree, which is required for GMYC analyses. For this, the sequence alignment was divided into two partitions in BEAUti v. 2.5.1 (Drummond et al. 2012): partition 1 contained the ITS, ETS, atpB-rbcL, rps16, trnL and trnL-F regions, while partition 2 only contained the psbA-trnH region. This was done to specify a different model of substitution for each partition following the results of the previously mentioned jModeltest analyses: the GTR model for partition 1 and the HKY substitution model for partition 2. The site model was set to Gamma with a category count of four and a substitution rate of 1.0. Insertion/deletion characters were not used for this analysis. All monophyletic groups that are present in the phylogenetic tree that resulted from the analysis of the combined nuclear and plastid dataset were constrained in the BEAST analysis to ensure that they were retained in the phylogeny after ultrametric conversion.

Table 2.

Primers and their respective sequences for amplification of the DNA regions used in this study.

Four xml files were created in BEAUti to correspond to four combinations of clock model and tree prior: 1) strict clock with Yule model (SY), 2) strict clock with coalescent constant population model (SC), 3) relaxed clock log normal with Yule model (RY), and 4) relaxed clock log normal with coalescent constant population model (RC). The marginal likelihood (ML) score of each combination was estimated using a Nested Sampling approach (Russel et al. 2018) to identify the combination that has the highest ML score for the dataset provided. A script for the Nested Sampling analysis was incorporated in each of the xml files using the MCMC-toNS-editor application within BEAST and with the following settings: a particle count of 10 with a subchain length of 10,000 and an Epsilon value of 1.0e-6. The modified xml files were run in BEAST on CIPRES. The highest likelihood score was obtained by the RC combination of clock model and tree prior (-15,074.38). Hence, an xml file with the RC combination was run in BEAST on CIPRES to generate the ultrametric tree for the GMYC analyses. The MCMC simulation chain length was set to 50M generations. Sampling was done at every 1000th generation. Tracer v. 1.7.1 (Rambaut et al. 2018) was used to inspect the effective sample sizes, which were at least 304 for each parameter, suggesting convergence. Trees were summarised using Tree Annotator v. 2.5.1 with a burn-in percentage of 10%, using the maximum clade credibility tree as the target tree type and common ancestor heights as the node heights.

GMYC Species Delimitation Analysis—The ultrametric tree obtained from the combined nuclear and plastid dataset was used as input for a single threshold GMYC analysis using the package “splits” (Ezard et al. 2009) in the R environment 4.0.1 (R Development Core Team 2020). GMYC aims to find nodes in ultrametric phylogenies that mark the point of transition between speciation and coalescence, whereby nodes before this transition reflect speciation events and more recent nodes indicate coalescence events (Pons et al. 2006; Fontaneto et al. 2011; Fujisawa and Barraclough 2013). Lineages that form putative species are marked by these points of transition. Before the analysis, outgroups were removed because this improves the performance of the species delimitation analysis (e.g. García-Melo et al. 2019). The function gmyc using a single threshold method was used to identify the nodes that mark the transition between speciation and coalescence events (Pons et al. 2006; Fujisawa and Barraclough 2013). The function gmyc.support was used to provide AIC-based support values (similar to probability values; Fujisawa and Barraclough 2013). A likelihood ratio test within the gmyc function was used for a subset of the dataset that is composed of only T. loheri s. l. specimens (65 samples) to test if the hypothesis that there is more than one species in the T. loheri s. l. complex is significantly better supported than the hypothesis that this complex is only composed of a single species.

PTP Species Delimitation Analysis—In contrast to the GMYC method, the PTP method does not require ultrametric trees. It models speciation and coalescence events directly using the branching rates (Zhang et al. 2013) and uses heuristic algorithms to classify the branches into those indicating species level processes and those indicating population level processes (Zhang et al. 2013; Tang et al. 2014). PTP assumes that the number of substitutions between species is significantly higher than the number of substitutions within species, resulting in differences between branch lengths associated with speciation events and those that are associated with coalescence events (Zhang et al. 2013). The non-ultrametric BI phylogenetic tree generated from the combined nuclear and plastid dataset was used for the PTP analysis. It was converted into a Newick file using FigTree and was uploaded to the bPTP server ( https://species.h-its.org/ptp/) which provides a Bayesian implementation of the PTP model for species delimitation (Zhang et al. 2013). The outgroups were removed from the analysis because this has been shown to provide more accurate results (Zhang et al. 2013) while identical sequences were retained in the dataset. The number of MCMC generations was set to 500,000, thinning was set to 100, and the burn-in was 10%. The likelihood scores were visually inspected for convergence which was indicated by consistent high likelihood scores in the trace file (Zhang et al. 2013). The accuracy of species delimitation at transition nodes was indicated by the value of their posterior probability (Zhang et al. 2013).

Identifying Morphological Support for GMYC and PTP Groups—To determine if the putative species delimited by the GMYC and PTP methods are diagnosably distinct in their morphology, a leaf shape geometric morphometric study of a subsample of 52 T. loheri s. l. herbarium specimens (Fig. 1) from the DNA sequence dataset was conducted using terminal leaflets and right-hand lateral leaflets (in abaxial view). For this, we used an outline analysis method (Adams et al. 2004), which is a commonly used geometric morphometric approach (e.g. Jensen et al. 2002; Viscosi et al. 2009; Klein et al. 2017) that captures leaf shape using margin geometry (Rohlf and Marcus 1993; Adams et al. 2004; Viscosi and Cardini 2011). At least three samples from each putative species were included (Fig. 1). One representative mature leaf from each herbarium specimen with at least two leaflets and without or with only minor damage was photographed for these analyses. For specimens with a missing or damaged right lateral leaflet, the left lateral leaflet was used. This was done by using the mirror image of the left lateral leaflet. For the few specimens with a missing leaflet apex, the shape of the leaf apex was redrawn in Adobe Photoshop CC v. 19.1.6 by estimating the position of the leaflet apex and using it as a guide to complete the leaflet shape. Furthermore, a few of the specimens had some level of herbivory resulting in holes in the lamina. These holes usually extended to a portion of the margin of the leaflet. This was addressed by reconstructing the missing areas of the lamina by redrawing the leaflet margin in Photoshop.

Outline analysis involves digitizing points along an outline of a structure and turning these point coordinates into quantitative variables using a mathematical transformation such as the elliptical Fourier transform (Giardina and Kuhl 1977; Kuhl and Giardina 1982; Bonhomme et al. 2014). The so-called Fourier coefficients obtained were then analysed in multivariate analyses (Adams et al. 2004). The terminal and lateral lamina of each specimen were cropped out from the leaf photographs separately and were converted to black images with white background using Photoshop. These photos were imported in “Momocs” (Bonhomme et al. 2014) in R using the function import_jpg which also extracted the outlines of the leaflet laminae from the photos. These outlines were converted to coordinates using the function Out. Two landmarks were placed on the leaflet outlines to guide the alignment of the lamina shapes: one at the apex of the petiolule, and the other one at the apex of the lamina of the associated leaflet. These landmarks allowed us to normalize outlines for position, size, rotation, and position of the first point. After alignment and scaling, the shape outlines were analyzed with an elliptical Fourier analysis using the function efourier. Normalization was set to ‘false’ since the standardization was previously done. Eleven harmonics were used for the lateral laminae and 12 for the terminal laminae, generating 44 and 48 harmonic coefficients respectively. Each set of harmonics had a cumulative power of 99% indicating that the number of harmonics generated was sufficient to describe the shape (Bonhomme et al. 2014). The 92 harmonic coefficients of the terminal and lateral laminae were combined to form a leaf shape dataset. A principal component analysis (PCA; method outlined below) of this dataset was performed and the scores of the first seven principal components, which captured 95% of the total variation, were used as synthetic shape variables. These seven variables were combined with a non-leaf shape dataset (10 vegetative characters other than leaf shape; Table 3) into a new matrix to form a combined morphological dataset comprising a total of 17 characters. Combining geometric morphometric data with other morphometric data has been shown to improve the correct classification of specimens into species (e.g. Ginter et al. 2012). Random forest analyses (Breiman 2001; Liaw and Wiener 2002) of this combined dataset were used to determine if the putative species delimited by GMYC and PTP are morphologically diagnosable. Because random forest analysis requires groups to be composed of at least two specimens, and four species delimited in the GMYC and PTP analyses were only represented by a single individual in the dataset (Fig. 1), these putative species were not included in the random forest analyses.

Fig. 1.

Bayesian inference phylogeny of the concatenated nuclear and plastid dataset. Clades with posterior probabilities of at least 0.95 are marked with red circles. Accessions illustrated with leaflet outlines (left: lateral leaflet; right: terminal leaflet) indicate plants included in the morphometric study. Plants collected from type localities are marked with an asterisk (^*). Highlighted groups are discussed in the text. Species delimited by the generalized mixed Yule coalescent (GMYC) and the Poisson tree processes (PTP) methods are indicated by blue (< 0.95 AIC-based support/posterior probability) and red bars (≥ 0.95 AIC-based support/posterior probability). Narrow delimitations (sensu stricto) of Groups A, D, F, and G are designated as A s. s., D s. s., F s. s., and G s. s., respectively. The islands or the regions within an island in the Philippines where the GMYC or PTP groups are found are indicated in the tree and in the map.

A random forest analysis produces multiple classification trees where each tree is generated from a bootstrapped training set typically containing two-thirds of the samples of the original dataset (Breiman 2001). A sample is run down each tree and is classified to the group to which most of the trees assign the sample. A k-folding cross-validation test is performed using the same classification process on the out-of-the-bag data (OOB; i.e. the testing set which contains the remaining one-third of the samples). The result of the cross-validation test is then compared to the group assigned to the sample in the original dataset and the OOB error rate is determined as the average proportion of the OOB samples that were incorrectly classified (Breiman 2001).

The random forest analyses were run in R using the “randomForest” package (Liaw and Wiener 2002). The analyses were executed using ntrees = 100,000 and mtry = 4. The optimal mtry value was determined by obtaining the square root of the total number of variables. Since the dataset has a different number of specimens in each group, down-sampling was used to balance group sizes (Chen et al. 2004). The number of specimens to be drawn from each group was set to the total number of specimens of the smallest group: three for the GMYC and two for the PTP-defined groups.

Table 3.

Non-leaf shape characters and character states used for morphological analyses. °These veins are referred to as second-order lateral veins by Pole (1991).

After inspecting the results of the phylogenetic analyses and the GMYC and PTP analyses, two additional species delimitation hypotheses were explored by performing random forest analyses using all 52 specimens for which morphometric data was generated. These analyses were aimed at testing morphological support for larger, more encompassing groups than those identified by the GMYC and PTP analyses. This approach was used because there are indications that these methods can result in false positives as a result of over splitting (Luo et al. 2018). Hypothesis 1 is a delimitation in which three species are recognized within the T. loheri s. l. complex (Groups A s. s., the core T. loheri s. l. clade, and H; Fig. 1). Hypothesis 2 recognizes two species within the complex (Groups H and A s. s. + the core T. loheri s. l. clade). Following the same down-sampling approach as previously mentioned (Chen et al. 2004), three specimens from each putative species were drawn for the random forest analysis for Hypothesis 1, and four for Hypothesis 2.

Principal Component and Coordinates Analyses—Since random forest analysis is a supervised clustering method, the conclusions obtained from this analysis directly rely on the assumption that the topology of the phylogenetic tree that was used for the GMYC and PTP analyses is correct. However, some nodes of this tree are poorly resolved or not strongly supported (Fig. 1). We therefore combined this approach with complementary multivariate approaches using the separate and combined morphological datasets to determine if groups of individuals that can be identified are morphologically distinct and therefore putative species under a unified species concept that uses morphological distinction as evidence of lineage separation (De Queiroz 2007). PCA and principal coordinates analysis (PCoA) are multivariate analysis approaches that can be used to reduce the dimensionality of a dataset through new projections of the variables that maximize the variance in a lower dimensional space (Abdi and Williams 2010). The leaf shape dataset was analyzed in R using the PCA function of “Momocs”, which wraps the base R prcomp and also uses the variance-covariance matrix. PCoA is a multi-dimensional reductional method similar to PCA but is capable of representing a mixture of qualitative and quantitative data through the use of an intermediate distance matrix. Distance matrices of the combined morphological dataset (composed of the non-leaf shape dataset and the leaf shape PCA scores) and the separate non-leaf shape data set were generated using Gower's distance (Gower 1971) in R studio with the function daisy with the metric type set to “gower” using the statistical package “cluster” (Kaufman and Rousseeuw 2009). The PCoAs of these datasets were performed in R studio using the function cmdscale.

All datasets and supplemental figures in this study are archived in the Dryad Digital Repository:  https://doi.org/10.5061/dryad.pg4f4qrks (Obico et al. 2021).

Results

Phylogenetic Analyses—Five DNA sequence datasets were compiled for this study (Table 4). The BI consensus trees obtained from the separate ITS and ETS datasets (Supplemental Figs. S1–S2) do not contain clades that are in strongly supported (i.e. posterior probabilities ≥ 0.95) incongruent positions, and both datasets were therefore combined. A BI analysis of this combined ITS-ETS dataset resulted in a tree that is better resolved and contains a larger number of strongly supported clades (Fig. S3). The BI consensus tree that was reconstructed from the concatenated dataset of five plastid regions (atpB-rbcL, psbA-trnH, rps16, trnL and trnL-F; Fig. S4) does not contain strongly supported clades that are incongruent with those in the ITS-ETS tree. We therefore also combined the plastid dataset with the ITS and ETS dataset, with the aim of using all available DNA sequence data for resolving the relationships among the Tetrastigma specimens that we included in our study. A BI analysis of this combined nuclear-plastid dataset resulted in a further increase in phylogenetic resolution and node support (Fig. 1) and this tree was therefore used for the species delimitation analyses.

Within the T. loheri s. l. complex, nine groups (labelled A–I) of T. loheri s. l. specimens can be identified that are each composed of accessions that either form clades in phylogenies obtained from all or most datasets, or that are at least consistently placed in each other's phylogenetic vicinity (Figs. 1, S1–S4). Group H is strongly supported by all datasets (Figs. 1, S1–S4). Groups B–G and I form a clade (core T. loheri s. l. clade) that is strongly supported by all datasets except the plastid dataset, but the relationships among these groups are poorly resolved and supported (Figs. 1, S4).

All specimens collected at the type localities of T. loheri, T. philippinense, T. stenophyllum, and T. trifoliolatum are found nested within the core T. loheri s. l. clade (Fig. 1). Those of the former three species are resolved as members of Group D and specimens from the type locality of T. trifoliolatum are members of Group G. However, two specimens identified as T. trifoliolatum from New Guinea and Peninsular Malaysia are placed outside of the T. loheri s. l. complex (Fig. 1). The combined nuclear and plastid phylogeny (Fig. 1) also shows some geographic structure at the level of islands or regions within the island of Luzon.

GMYC and PTP Species Delimitation Results—The number of species identified in the GMYC analysis is 33 with a 95% confidence interval between 11–41. Within the T. loheri s. l. complex, GMYC identified 12 species. Nine of these are composed of more than one specimen (Fig. 1). Groups B, C, E, G, H, and I were all delimited as species by the GMYC modelling. The delimitations of three additional putative species identified in the analysis each closely correspond to Groups A, D, and F (here referred to as Groups A s. s., D s. s., and F s. s.), but do not include one of their specimens, which are each individually segregated as species. Of the other Tetrastigma species that were included in the combined nuclear and plastid dataset and represented by more than one specimen, only T. harmandii, T. papillosum, and T. scariosum were recovered in the GMYC analyses. The delimited species all have AIC-based support values less than 0.95. A likelihood ratio test that was conducted with a dataset that only contains T. loheri s. l. specimens to specifically test the hypothesis that this complex is composed of more than one species resulted in a likelihood value of 597.1079 for the hypothesis that there are no distinct species groups within T. loheri s. l. and a value of 601.3039 for the hypothesis that more than one species group exists. The latter hypothesis is better supported (p = 0.015).

Table 4.

Details of the DNA sequence datasets used in this study.

The results of the PTP delimitation estimated a mean number of species of 32.55 with a 95% confidence interval between 17 and 61. PTP delimited nine species within T. loheri s. l. which consist of six groups of specimens and three individual specimens (Fig. 1). The six groups are BCD, FI, A s. s., E, H, and G s. s. Group BCD is the combination of Groups B, C, and D s. s. Likewise, Group FI is the combination of F and I. Of the remaining species included in the dataset that are represented by more than one specimen, only T. harmandii and T. lawsonii (King) Burkill were identified by PTP. Only some of the PTP-delimited species outside of the T. loheri s. l. group are strongly supported (i.e. pp ≥ 0.95) and all of these consist of individual specimens (Fig. 1).

Random Forest Analyses—Random forest analyses were used to determine if the putative species delimited with the GMYC and PTP methods are diagnosably distinct when using the combined morphological dataset. These analyses resulted in a 66% out-of-the-bag (OOB) error of classifying specimens to GMYC groups within the T. loheri s. l. complex. At least 20% of the specimens of each GMYC group was misclassified by the model, except for specimens of Group H (0% classification error; Table 5). The OOB estimate of error of classifying samples to T. loheri s. l. groups as delimited by PTP was 60% and each PTP group had a classification error of at least 40%, except for specimens of Group H (0% classification error; Table 6).

Table 5.

Random forest confusion matrix for putative species delimited with the generalized mixed Yule coalescent method within Tetrastigma loheri s. l. Classification error is the percentage of specimens that were misclassified.

Alternative species delimitation Hypothesis 1 (i.e. three species: A s. s., core T. loheri s. l. clade, H) had a 19% OOB error rate and Hypothesis 2 (i.e. two species: H, A + core T. loheri s. l. clade) had a 10% OOB error rate. Each group had a classification error of at least 18% for Hypothesis 1 and 10% for Hypothesis 2, except for Group H which had 0% classification error for both hypotheses (Tables 7, 8).

Principal Component and Coordinates Analyses—The first principal component axis of the PCA of the leaf shape data (Fig. 2) explains 52.9% of the total variation. This axis shows variation in leaflet aspect ratio, i.e. roundness (Fig. 2). The second principal component axis, which explains 23.7% of the total variation, shows variation in the shape of the apex of the lateral leaflets and their bilateral symmetry (Fig. 2). In the PCoA plot of the combined morphological data, both the first and second axes, which explain 12.2% and 10.4% of the variation respectively, seem to show variation in leaflet aspect ratio along their axes if the individual leaflet shapes are plotted along with the position of the specimens in the PCoA plot (Fig. 3). Neither the PCA ordination plot of the leaf shape data (Fig. 2) nor the PCoAs of the non-leaf shape dataset (not shown) and the combined morphological dataset (Fig. 3) show morphologically distinct clusters of specimens nor a clear clustering of specimens from the same GMYC or PTP groups or from the same island or region in the Philippines. Although most members of H and GMYC group B are closer to each other than individuals from other groups in the PCoA plot of the combined morphological data (Fig. 3), subsequent linear discriminant analysis based on the GMYC groups showed poor discrimination of these two groups (results not shown).

Discussion

GMYC and PTP Species Delimitation—The results of a likelihood ratio test as part of a GMYC analysis of a combined nuclear and plastid DNA sequence phylogeny in which only specimens of T. loheri s. l. are included suggests that this complex is composed of more than one species (null hypothesis of no distinct species groups within T. loheri s. l. rejected: p = 0.015). In fact, the GMYC analysis of the complete dataset that also includes representatives of other Tetrastigma species suggests that there may be 12 species within T. loheri s. l. (Fig. 1). Also, the results of the PTP analysis indicate the presence of more than one species within the complex (Fig. 1). The delimitation of the nine T. loheri s. l. species that resulted from the PTP analysis is largely congruent with those delimited by the GMYC analysis, but the former is mostly more broadly delineated and therefore fewer in number. Whereas the GMYC method identified specimen Groups B, C, and D s. s. as separate species, the PTP model grouped these together into Group BCD. Similarly, the GMYC model identified Groups I, F s. s., and a single accession of a specimen from Central Luzon (Barcelona 3868) as distinct, but the results of the PTP analysis instead suggest that these specimens are conspecific (i.e. Group FI). Only Group G was more narrowly delimited by the PTP method than by the GMYC method (Fig. 1).

Table 6.

Random forest confusion matrix for putative species delimited with the Poisson tree process method within Tetrastigma loheri s. l. Classification error is the percentage of specimens that were misclassified.

Table 7.

Random forest confusion matrix for putative species within Tetrastigma loheri s. l. for alternative species delimitation Hypothesis 1 (see text). Classification error is the percentage of specimens that were misclassified.

Table 8.

Random forest confusion matrix for putative species within Tetrastigma loheri s. l. for alternative species delimitation Hypothesis 2 (see text). Classification error is the percentage of specimens that were misclassified.

Fig. 2.

Ordination plot from a principal component analysis of the leaf shape data of 52 individuals of Tetrastigma loheri s. l. using the first two principal components. Specimens collected from type localities of species are indicated with their species names. Letters indicate the GMYC-PTP groups in which the specimens were placed (Fig. 1). Barcelona 3667, 3726, and 3868 are referred to in the Discussion. The paired leaflet outlines refer to the leaflet lamina shapes in which the top shape corresponds to the terminal leaflet and the bottom shape to the lateral leaflet. The color of the symbols indicates the Philippine regions or islands where the specimens were collected.

Although both species delimitation methods suggest that the T. loheri s. l. group is composed of multiple species and similar delimitation hypotheses were obtained, the differences between these delimitations indicate that care should be taken when interpreting these results. They might suggest that the performance of the GMYC and PTP modelling was not optimal in the present study. This is further indicated by the large confidence intervals of the GMYC and PTP species estimations (i.e. GMYC 11–41 species; PTP 17–61 species). In addition, statistical support for the putative species within the T. loheri s. l. complex was poor in both delimitations (i.e. all AIC support values and posterior probabilities < 0.95). Furthermore, only a few species of the other Tetrastigma species included in the analyses represented with more than one individual were correctly identified by the GMYC (i.e. T. harmandii, T. papillosum, T. scariosum) and PTP methods (i.e. T. harmandii, T. lawsonii). Non-optimal performance of the GMYC and PTP modelling could be a result of errors in phylogenetic reconstruction or insufficient phylogenetic resolution due to a lack of informative characters. However, although indeed only a few clades obtained high posterior probabilities in the combined plastid and nuclear phylogeny (e.g. Groups A, H, core T. loheri s. l. clade; Fig. 1), none of these was identified as a potential species by the GMYC or PTP analyses with high statistical support.

Fig. 3.

Ordination plot from a principal coordinate analysis of the combined morphological dataset (leaf shape PCA scores and non-leaf shape data) of 52 individuals of Tetrastigma loheri s. l. using the first two principal components. Specimens collected from type localities of species are indicated with their species names. Letters indicate the GMYC-PTP groups in which the specimens were placed (Fig. 1). Barcelona 3667, 3726, and 3868 are referred to in the Discussion. The paired leaflet outlines refer to the leaflet lamina shapes in which the left shape corresponds to the lateral leaflet and the right shape to the terminal leaflet. The color of the symbols indicates the Philippine islands or regions where the specimens were collected.

Finding Morphological Support for GMYC and PTP Groups—The present study aimed to determine if the putative species delimited by GMYC and PTP are diagnosably distinct in vegetative morphological characters. Such groups could be recognized as species under the unified species concept (De Queiroz 2007), assuming phylogenetic as well as morphological distinction as evidence of lineage separation.

To determine if the GMYC and PTP-delimited groups are diagnosable using vegetative morphology, random forest analyses were carried out using a combined vegetative morphological dataset (i.e. leaf shape data and several other vegetative characters). Because the random forest method requires groups to be composed of at least two specimens, putative species to which only one specimen was assigned (i.e. singletons) could not be included in these analyses. The four singletons identified within T. loheri s. l. (Barcelona 3667, Barcelona 3726, Barcelona 3868, and Chen and Lu 158) might represent species that are rare and additional representatives of these lineages need to be studied to clarify their taxonomic status. However, although Chen and Lu 158 was not available for the morphometric study, the results of the PCoAs and PCAs do not indicate that the other three specimens are morphologically distinct from the other T. loheri s. l. specimens (Figs. 2, 3).

The overall random forest classification errors for the GMYC and PTP delimitation hypotheses (excluding singletons) of the T. loheri s. l. complex are high (66% and 60% respectively). However, this error rate decreased when a three-species hypothesis was explored (19%; Table 7). The two-species hypothesis had the lowest error rate (10%; Table 8). The specimens of Group H obtained a 0% classification error in all confusion matrices (Tables 5–8). Although this result may suggest that Group H shows some level of morphological distinction, some specimens from other groups were incorrectly classified to Group H (Tables 5–8). Furthermore, Group H specimens did not form a morphologically distinct cluster in the PCA and PCoA ordination plots (Figs. 2, 3) and random forest analyses using the separate leaf shape and non-leaf shape datasets (not shown) resulted in high classification errors. Therefore, in conclusion, the vegetative characters used in this study do not appear to indicate that the putative species identified by the GMYC or PTP methods are diagnosably distinct (Tables 5, 6), even when broader species delimitation hypotheses are considered (Tables 7, 8).

Principal Component and Coordinates Analyses—PCAs and PCoAs of the separate (e.g. Fig. 2) and combined (Fig. 3) morphological datasets were carried out to determine if there are morphologically distinct groups of specimens that were not identified by the supervised clustering (i.e. random forest) analyses. Such groups might be found if the T. loheri s. l. complex is composed of more than one species, but these species were not recovered by the GMYC and PTP analyses, for example, as a result of incorrect phylogenetic relationships in the topology used for these species delimitation analyses. However, none of the PCAs and PCoAs revealed morphologically distinct groups of specimens and therefore also do not indicate that more than one species should be recognized within the T. loheri s. l. complex in the Philippines.

Species Delimitation in the T. loheri s. l. complex—Our results do not recognize more than one species within the T. loheri s. l. complex in the Philippines with the species delimitation criteria used, because they do not indicate the presence of groups that are both phylogenetically and morphologically distinct. This is corroborated by the results of the PCA and PCoA analyses of T. loheri s. l. specimens and indicates that this lineage is most likely a species that is very variable in its vegetative morphology. This is well illustrated by the two Philippine specimens that were collected at the type locality of T. trifoliolatum (Obico 983 and Obico 984). Both were collected in each other's vicinity and are resolved as each other's closest relatives (Figs. 1, S1–S4). However, the two T. trifoliolatum specimens are substantially different in their vegetative morphology as indicated by their different positions in the various morphometric ordination plots and their leaflet shape outlines (Figs. 2, 3). Morphological differences that are associated with habitat, elevation, or other ecological factors were not observed (results not shown). Overall, it therefore seems most likely that the large morphological disparity in the vegetative characters that were studied is a result of phenotypic plasticity, although heteroblasty/heterophylly (e.g. Baumgartner et al. 2020) and local adaptation cannot be entirely excluded.

Although the available evidence suggests that the T. loheri s. l. complex might only be composed of a single species in the Philippines, it is possible that future studies of reproductive or other characters that were not included in this study will provide evidence in support of recognizing more than one species. This should be taken into account when taxonomic conclusions are drawn from this study. In addition, despite that the GMYC and PTP analyses did not result in strongly supported species delimitation hypotheses and that none of the GMYC and PTP groups were found to be morphologically distinct, some of the inferred groups might be biologically meaningful, because they are composed of specimens that are mostly from the same island (e.g. Groups A s. s., F, and I; Fig. 1) or from the same region within an island (e.g. Groups B, C, D; Fig. 1). Because these islands and some of the regions are separated by geographic features that may prevent or significantly reduce gene flow among them (e.g. bodies of water), it is possible that the correspondence between the species delimitation groups and biogeographic patterns in the data is a result of their reproductive isolation. This could therefore mean that these groups are morphologically cryptic or incipient species. Alternatively, however, the biogeographic patterns in the phylogenies might not be an indication of the existence of multiple species or incipient species within the T. loheri s. l. complex, but rather provide an explanation for obtaining species delimitation hypotheses that are too narrowly delineated as a result of over splitting (Luo et al. 2018). That is because geographic structuring of intraspecific genetic variation might result in identifying a partially genetically isolated population as a separate species (Papadopoulou et al. 2008; Fujisawa and Barraclough 2013; Luo et al. 2018). Although reduced gene flow between some of the GMYC and PTP groups as a result of geographic barriers is an intuitively appealing explanation especially in an archipelago like the Philippines, geographic structuring can also be explained by spatial autocorrelation: individuals of a species that are living closer together are expected to be more genetically similar than individuals that live further apart (Meirmans 2012).

Of all possible putative species considered, most support exists for recognizing Group H as taxonomically distinct at the species level. Resolved as more distantly related to the other T. loheri s. l. groups than they are to each other, Group H is phylogenetically distinct and forms a strongly supported clade (Figs. 1, S1–S4). In addition, Group H is biogeographically unique among the GMYC and PTP groups, because it is composed of specimens from islands that are relatively distant from each other: Leyte, Mindanao, and Panay (Fig. 1). In addition, although Group H individuals have been found to be sympatric with T. loheri s. l. specimens from other groups on each of these islands, they are genetically distinct. For example, Obico 986, which is a member of Group H, was collected in the same forest on Leyte where two of the Group G individuals (Obico 983 and Obico 984) were found. Although such patterns of sympatry could be a result of recent long-distance dispersal, they might also indicate reproductive isolation in sympatry and, as such, be considered as evidence for recognizing Group H as a distinct species (e.g. Ma et al. 2019).

Taxonomic Implications—If the T. loheri s. l. complex indeed only consists of one species in the Philippines and considering that Philippine specimens identified as T. diepenhorstii, T. loheri, T. philippinense, T. stenophyllum, and T. trifoliolatum are resolved as its members, taxonomic changes to the species-level classification of Tetrastigma might be required. Such changes require certainty regarding the correct application of taxonomic names. This is problematic because type specimens were not included in this study. In an attempt to mitigate this, specimens from the type localities of four of the five species that have been reported from the Philippines and that were identified as putative members of the T. loheri s. l. complex because of morphological similarities were included in the analyses presented here. These specimens show the diagnostic morphological features as outlined in the species protologues (Gagnepain 1910; Merrill 1912, 1914, 1916). Unfortunately, specimens from the type locality of T. diepenhorstii in Sumatra (Miquel 1861) were unavailable for study.

Four specimens identified as T. trifoliolatum were included in this study. Two Philippine specimens (Obico 983 and Obico 984) are nested within Group G of the T. loheri s. l. complex, whereas two specimens from Peninsular Malaysia (Wen 8350) and New Guinea (Wen 10758) form a clade with specimens of T. papillosum and T. crenatum (Fig. 1). This raises doubt about the correct identification of these specimens. Because the Philippine specimens were collected from the type locality of T. trifoliolatum in Leyte, it is most likely that they represent this species. Unfortunately, specimens of Wen 8350 and Wen 10758 at US were not accessible for this study. Hence, it was not possible to verify their identification. If Obico 983 and Obico 984 indeed represent T. trifoliolatum and if future studies would provide evidence in support of recognizing Group G as a distinct species, then the name T. trifoliolatum is available for this species.

Tetrastigma trifoliolatum was synonymized with T. diepenhorstii by Latiff (2001), who compared the types of both names and considered them to be conspecific. The two specimens identified as the latter species that were included in our analyses are resolved as members of the T. loheri s. l. complex, although they are placed in different groups (i.e. Groups A and E; Fig. 1). One specimen (Chen and Lu 158), which was identified as a singleton in the GMYC and PTP analyses, is from Indonesia and the other (Wen 8261) was collected in the Philippines (Mt. Makiling, Southern Luzon). Wen 8261 could not be located at US, but the vegetative morphology of Chen and Lu 158 (PE photo!) is in agreement with the description in the protologue of T. diepenhorstii (Miquel 1861). It is therefore possible that the type of T. diepenhorstii is indeed a member of the T. loheri s. l. complex, although this remains to be confirmed.

Specimens from the type localities of T. loheri, T. philippinense, and T. stenophyllum are resolved as closely related to each other. They are placed in Group D s. s. (GMYC) or Group BCD (PTP; Fig. 1). Under the assumption that these specimens indeed represent these three species, T. loheri, T. philippinense, and T. stenophyllum are therefore synonyms even if the narrowest species delimitation hypothesis provided by the GMYC and PTP analyses is adopted. Further, if T. diepenhorstii, T. loheri, T. philippinense, T. stenophyllum, and T. trifoliolatum are conspecific, T. diepenhorstii is the correct name for T. loheri s. l., because this name has priority (Turland et al. 2018). However, because neither the type specimen nor specimens from the type locality of T. diepenhorstii were included and because this name has not been used for Philippine plants, other than for Wen 8261 (Pelser et al. 2011), which was not available for study, we suggest the name T. loheri to refer to the Philippine representatives of this taxon in the interim, because this name has priority over T. philippinense, T. stenophyllum, and T. trifoliolatum.