First molecular phylogenetic insights into the evolution of Eriocaulon (Eriocaulaceae, Poales)

Eriocaulon is a genus of c. 470 aquatic and wetland species of the monocot plant family Eriocaulaceae. It is widely distributed in Africa, Asia and America, with centres of species richness in the tropics. Most species of Eriocaulon grow in wetlands although some inhabit shallow rivers and streams with an apparent adaptive morphology of elongated submerged stems. In a previous molecular phylogenetic hypothesis, Eriocaulon was recovered as sister of the African endemic genus Mesanthemum. Several regional infrageneric classifications have been proposed for Eriocaulon. This study aims to critically assess the existing infrageneric classifications through phylogenetic reconstruction of infrageneric relationships, based on DNA sequence data of four chloroplast markers and one nuclear marker. There is little congruence between our molecular results and previous morphology-based infrageneric classifications. However, some similarities can be found, including Fyson’s sect. Leucantherae and Zhang’s sect. Apoda. Further phylogenetic studies, particularly focusing on less well sampled regions such as the Neotropics, will help provide a more global overview of the relationships in Eriocaulon and may enable suggesting the first global infrageneric classification. Electronic supplementary material The online version of this article (10.1007/s10265-019-01129-3) contains supplementary material, which is available to authorized users.


Introduction
Eriocaulon L., commonly known as pipeworts, is a cosmopolitan genus of ephemeral and perennial aquatic and wetland plants of the Eriocaulaceae family (Poales). The genus includes c. 470 species (WCSP 2019) and is most speciesrich in Asia (c. 220 species), the Americas (c. 122 species) and Africa (c. 111 species), with its centres of diversity in the tropics (Stützel 1998). Species of Eriocaulon primarily grow in seasonal or permanent wetlands while some inhabit shallow rivers and streams with an apparent adaptive morphology of elongated submerged stems. Two subfamilies are recognised in Eriocaulaceae, i.e. Eriocauloideae with diplostemonous flowers and glandular petals, and Paepalanthoideae with isostemonous flowers and eglandular petals (Giulietti et al. 2012;Ruhland 1903). Together with the African genus Mesanthemum Körn., which was recently revised by Liang et al. (2019), Eriocaulon is placed in subfamily Eriocauloideae. Subfamily Paepalanthoideae is largely restricted to the Americas.
Despite the ecological importance of Eriocaulon as a species-rich genus of wetland plants, no attempts have been made to reconstruct a molecular phylogeny for the genus. Only a few Eriocaulon species have been included in the sampling of family level studies (e.g. de Andrade et al. 2010;Giulietti et al. 2012). A molecular phylogenetic study including a broad sampling covering much of the taxonomic, morphological and geographic variation within the genus is needed to assess whether the infrageneric taxa suggested in the existing regional infrageneric classifications of Eriocaulon circumscribe monophyletic groups. It is a first step in providing insights into the evolution of the genus and to enable establishing a new infrageneric classification for the whole genus in the future.
None of the published regional infrageneric classifications have yet been scrutinised using molecular phylogenetic data. A molecular study by de Andrade et al. (2010) on Eriocaulaceae included just five species of Eriocaulon while the study of Giulietti et al. (2012) included just four species. Of the species sequenced in these studies, E. cinereum R.Br. is the only one that has been included in the published infrageneric classifications. The objectives of this study are to: (1) construct a molecular phylogeny of Eriocaulon, and (2) critically assess the existing regional infrageneric classifications of Eriocaulon.
The PHYC gene (a distinct member of the phytochrome gene family) was selected as a nuclear marker, based on its phylogenetic utility as a single or low copy nuclear locus (Mathews and Donoghue 1999;Samuel et al. 2005). Fragments of a part of exon 1 of PHYC were amplified by PCR using Comm_PHYC_P1F (Hertweck et al. 2015) and the newly designed AlisPHYC-1R (5′-GCA TCC ATTTCMACA TCY TCCCA). The PCR cycling conditions were 94 °C for 90 s; then 35 cycles of 94 °C for 45 s, 60 °C for 30 s, 72 °C for 90 s; and finally, 72 °C for 10 min. The fragments obtained were digested with ExoSAP-IT and directly sequenced.
The PCR products were cleaned using ExoSAP-IT (GE Healthcare, Piscataway, NJ, USA) purification, and then amplified using ABI PRISM Big Dye Terminator v.3.1 (Applied Biosystems, Foster City, CA, USA) using the same primers as those used for the PCR amplifications. DNA sequencing was performed with an ABI PRISM 377 DNA sequencer (Applied Biosystems). Automatic base-calling was checked by eye in Genetyx-Win v.3 (Software Development Co., Tokyo, Japan). All sequences generated in the present study have been submitted to the DNA Data Bank of Japan (DDBJ), which is linked to GenBank, and their accession numbers and voucher specimen information are presented in Table S2.

Molecular phylogenetic analyses
Sequences were aligned using MAFFT v.7.058 (Katoh and Standley 2013) and then inspected manually. Analyses were independently performed for ptDNA (matK, rbcL, rpoB, rpoC1) and PHYC datasets respectively to identify possible incongruences between different genomic regions. All 199 ingroup and the 13 outgroup accessions were represented in the ptDNA dataset, while 55 ingroup and five outgroup accessions were represented in the PHYC dataset. The ptDNA dataset consisted of concatenated gene alignments with 145 or 68% of accessions represented for matK, 197 or 93% for rbcL, 122 or 58% for rpoB and 41 or 19% for rpoC1.
Phylogenies were reconstructed using maximum parsimony (MP), maximum likelihood (ML), and Bayesian inference (BI; Yang and Rannala 1997). In the MP analysis in PAUP* v.4.0b10 (Swofford 2002), a heuristic search was performed with 100 random addition replicates with treebisection-reconnection (TBR) branch swapping, with the MulTrees option in effect. The MaxTrees option was set at 100,000. Bootstrap analyses (Felsenstein 1985) were performed using 1,000 replicates with TBR branch swapping and simple addition sequences. The MaxTrees option was set at 1,000 to avoid entrapment in local optima. For the ML analysis, the RAxML BlackBox online server (https ://www.raxml -ng.vital -it.ch/) was used, which supports GTR-based models of nucleotide substitution (Stamatakis 2006). The maximum likelihood search option was used to find the best-scoring tree after bootstrapping. The gamma model of rate heterogeneity was selected. Statistical support for branches was calculated by rapid bootstrap analyses of 100 replicates (Stamatakis et al. 2008).
BI analyses were conducted with MrBayes v.3.2.2 (Ronquist and Huelsenbeck 2003; Ronquist et al. 2012) run on the CIPRES portal (Miller et al. 2010) after the best models had been determined in MrModeltest v.3.7 (Nylander 2002); these models were GTR + I + G and GTR + G for ptDNA and PHYC datasets, respectively. Analyses were run for 6,335,000 and 1,500,000 generations for ptDNA and PHYC datasets, respectively, until the average standard deviation of split frequencies dropped below 0.01, sampling every 1,000 generations and discarding the first 25% as burn-in.

Molecular dating
A species tree was used to conduct a molecular dating analysis. A multispecies coalescent method (Heled and Drummond 2009) implemented in BEAST v.1.7.2 (Drummond et al. 2006;Drummond and Rambaut 2007) was performed to reconstruct a species tree. *BEAST was run using a multilocus dataset (ptDNA and PHYC) utilising all 212 ingroup and outgroup samples assigned to the 84 operational taxonomic units (OTUs) that were retrieved as clades in the phylogenetic analyses above. For the purposes of this analysis, species resolved as non-monophyletic or that contained multiple lineages are represented multiple times in the resulting tree (i.e. E. cinereum R.Br., E. latifolium Sm., E. nepalense J.D.Prescott ex Bong., E. plumale N.E.Br. and E. setaceum L.).
A relaxed molecular clock as implemented in BEAST v.2.4.4 (Drummond et al. 2006) was used and run on the CIPRES portal (Miller et al. 2010). Uncorrelated lognormal distributed substitution rates for each branch were used. The tree was rooted by constraining Eriocaulaceae. Previous divergence time estimates between Eriocaulaceae and Xyridaceae of 105 mya (million years ago) provided by Janssen and Bremer (2004), Bouchenak-Khelladi et al. (2014) and Hertweck et al. (2015) were used as a calibration point. These dates were set as a mean age with stdev = 0.1 and a normal distribution. A Yule speciation process was used as tree prior. The default settings of BEAUti v.2.4.4 were used for the other parameters. Two runs of ten million generations of the MCMC chains were run, sampling every 1,000 generations. Convergence of the stationary distribution was checked by visual inspection of plotted posterior estimates using Tracer v.1.6 (Rambaut et al. 2014). After discarding the first 1,000 trees as burnin, the samples were summarised in the maximum clade credibility tree using TreeAnnotator v.1.6.1 (Drummond and Rambaut 2007) with a posterior probability (PP) limit of 0.5 and summarizing mean node heights. The results were visualised using FigTree v.1.3.1 (Rambaut 2009).

Molecular phylogeny
The ptDNA dataset for four genes included 4,445 aligned characters, of which 889 were parsimony informative. Analysis of this dataset yielded the imposed limit of 100,000 MP trees (tree length = 2,590 steps; consistency index = 0.64; retention index = 0.88). The strict-consensus MP tree, the RAxML tree, and the MrBayes BI 50% consensus tree showed no incongruent phylogenetic relationships; thus only the BI tree is presented here (Fig. 2a). Eriocaulon is broadly divided into two lineages: clade I and clades II-XII. The clade II is resolved as sister to clades III-XI. Clade III is resolved as sister to clades IV-XII. The relationships among clades IV-XII are less well resolved, except the weakly supported clades VIII-IX, yet each clade is highly supported. Singleton X is differentiated from clades XI-XII.
The PHYC dataset included 979 aligned characters, of which 324 were parsimony informative. Analysis of this dataset yielded the imposed limit of 100,000 MP trees (tree length = 1,181 steps; consistency index = 0.53; retention index = 0.81). The strict-consensus MP tree, the RAxML tree and the MrBayes BI 50% consensus tree showed no incongruent phylogenetic relationships; thus, only the BI tree is presented here (Fig. 2b). The labelling of PHYC tree follows the ptDNA tree. Eriocaulon is broadly divided into two lineages: Clade I and clades/singletons II-III, V-IX and XI-XII. Singleton III is resolved as sister to clades/singletons II-XII. Singleton VI and clade VII are retrieved as sister lineages, as are clades/singletons II, V, VIII-IX and XI-XII. The relationships in the latter group are less resolved. Clade II is strongly supported. Clade XII plus Eriocaulon_schimperi_K3065 which belongs to clade XI in the ptDNA analysis are strongly supported as a natural group. Members of clade XI except Eriocau-lon_schimperi_K3065 are retrieved as a clade.

Molecular dating
The divergence time for each clade was estimated using the calibration point between Eriocaulaceae and Xyridaceae of 105 mya provided by Janssen and Bremer (2004), Bouchenak-Khelladi et al. (2014) and Hertweck et al. (2015). The most recent common ancestor (MRCA) of the Eriocaulaceae family was estimated as early Paleogene with the Eriocauloideae MRCA as mid-Paleogene. The approximate divergence time for the MRCA of Eriocaulon was estimated as late Paleogene to early Neogene (21.66 mya; 95% HDP = 15.88-28.36 mya) (Fig. 3). Most of the species diversity of Eriocaulon appears to have originated in the the last 10 mya (Fig. 3).

Phylogeny and systematics of Eriocaulon
We reconstructed the phylogenetic history of Eriocaulon using both ptDNA and PHYC datasets with the aim of assessing the existing regional infrageneric classifications (Table S1). Although our taxon sampling is not sufficiently comprehensive to cover Mueller's (1859) sectional classification for Australian species of Eriocaulon, selected species listed in the infrageneric classifications proposed by Fyson (1919Fyson ( , 1921Fyson ( , 1922, Ma (1991Ma ( , 1997, Ansari andBalakrishnan (1994, 2009) and Zhang (1999) were sampled (Table S1, S2). Here, using the ptDNA tree (Fig. 2a), we discuss whether and how the results support these infrageneric classifications of as well as the previous molecular phylogeny of de Andrade et al. (2010). There is little congruence between our molecular results and previous morphology-based infrageneric classifications. However, some similarities can be found, as detailed below.
In de Andrade et al. (2010) ptDNA tree, Eriocaulon cinereum was retrieved as sister to the other four species including E. decangulare L. Our ptDNA tree recovered a similar topology in which E. cinereum of clade II branches off before E. decangulare of clade IV (Fig. 2a). Eriocaulon cinereum belongs to sect. Leucantherae Fyson characterised by pale anthers and a smooth seed coat (Fyson 1919(Fyson , 1921(Fyson , 1922Zhang 1999), and recognised by Ansari andBalakrishnan (1994, 2009) as their sect. XII. This group is represented by clade II of the ptDNA tree (Fig. 2a), and hence upheld by both morphological and molecular evidence. Ma (1991Ma ( , 1997 classified Chinese Eriocaulon into subgen. Trimeranthus of 27 species and subgen. Eriocaulon accommodating E. decemflorum. Although neither subgen. Trimeranthus nor most of its sections or series are supported in our results, it is noteworthy that E. decemflorum is retrieved as a single species lineage ( Fig. 2a clade X). It was also placed as the only member of sect. Nasmythia by Zhang (1999) based on its dimerous flowers and seed coat structure. Ansari andBalakrishnan (1994, 2009) proposed an infrageneric classification of the Indian species of Eriocaulon into twelve sections, primarily based on seed surface structure. These are mostly not supported by our molecular analysis. For instance, E. nepalense, E. parviflorum (Fyson) R. Ansari & N.P. Balakr. and E. xeranthemum Mart. are grouped in their sect. III, but are scattered in clades I, XI and XII of the ptDNA tree (Fig. 2a). Similarly, E. truncatum and E. hamiltonianum and are grouped in their sect. VII but are here placed in clades XI and XII, respectively (Fig. 2a). Zhang (1999) carried out a morphology-based cladistic analysis of the 71 East Asian species studied. The resulting cladograms divided the species into four groups. None of these groups are supported in our ptDNA and PHYC trees (Fig. 2). However, noteworthy in our results is the sister relationship between clade I and the rest of Eriocaulon (Fig. 2). Clade I accommodates relatively robust and large species, i.e. E. australe R.Br., E. cuspidatum Dalzell and E. sexangulare L. Zhang (1999)

Species distribution and taxonomy
Some species of Eriocaulon are known to have a wide distribution in the Old World tropics, such as E. cinereum and E. setaceum (Cook 1996(Cook , 2004. In the present study, E. cinereum is divided into two lineages, one from Africa and the other from Asia, although both fall within clade II (Fig. 2a). Similarly, samples of E. setaceum from Africa and Asia showed genetic variation ( Fig. 2a clade XII). On the other hand, no significant differentiations are observed among samples of E. truncatum from Africa and Asia ( Fig. 2a clade XI). Species such as E. cinereum and E. truncatum are common in rice fields (Cook 1996), probably contributing to their widespread distribution around the world.
Clade I includes a subclade of 12 accessions of Eriocaulon australe and E. sexangulare. Prajaksood et al. (2012) reduced E. australe to a variety of E. sexangulare (E. sexangulare var. australe (R.Br.) Praj. & J.Parn.). These taxa differ in E. australe having hairy leaves, sheaths, involucral bracts and receptable (Prajaksood et al. 2012;Zhang 1999). From our results it appears that this character may not be phylogenetically informative, and therefore, the varietal status of E. australe is supported.
Clade VII comprises Eriocaulon alpestre Hook.f. & Thomson ex Körn., E. buergerianum Körn., E. sikokianum Maxim., E. hondoense, E. miquelianum and E. taquetii with limited genetic variation among samples (Fig. 2 clade VII). This clade corresponds to subgen. Spathopeplus sect. Apoda of Zhang (1999), a mainly Asian group of species with female sepals connate into a spathe and seeds with T-shaped projections. Our results reflect the taxonomic complexity of this group in Japan, e.g. E. sikokianum Maxim. (accepted name: E. miquelianum Körn.) and E. hondoense (accepted name: E. taquetii Lecomte), while E. buergerianum_ TCMK_493_K needs critical re-identification because this species is clearly diagnosed by floral morphology (Ma et al. 2000;Miyamoto 2015 Phillips 1998Phillips , 2010Phillips , 2011. Further phylogenetic work is required to investigate the morphology-based hypotheses for species delimitation in Asian and African species of Eriocaulon. Several of the 20 samples of Eriocaulon from Cambodia are phylogenetically unique and distinguished from other known species in both ptDNA and PHYC trees (Fig. 2). An in-depth morphological analysis may reveal whether the collections belong to undescribed species.

Evolutionary history of Eriocaulon
Our results show that Eriocaulon originated in the late Paleogene to early Neogene (Fig. 3), and most species diversity originated in the last 10 mya. With Eriocaulon occurring in (permanent or ephemeral) wetlands across the tropics, the increased speciation during this time may be due to drift arising from fragmentation of suitable habitats associated with aridification since mid-Miocene. During the same period, Poales lineages like Cyperaceae and Poaceae that evolved adaptation to aridification (e.g. C4 photosynthesis, growing in non-wetland habitats) exhibited increased diversification (Bouchenak-Khelladi et al. 2014).

Future perspectives
For the time being, we refrain from suggesting a new infrageneric classification of this polymorphic and widespread genus until the morphological characters used in previous studies (e.g. seed surface structure, anther colour, floral structure; Ansari andBalakrishnan 1994, 2009;Zhang 1999) can be thoroughly investigated for their phylogenetic informativeness. Further phylogenetic studies, particularly focusing on less well sampled regions such as the Neotropics, will help provide a more global overview of the relationships in Eriocaulon and may enable suggesting the first global infrageneric classification. Further research may also aid towards understanding closely related species groups in Africa and Asia. This study should be seen as a step towards achieving the aim of a natural classification.