{'fraction_non_alphanumeric': 0.03449995524468371, 'fraction_numerical': 0.024244594197153562, 'mean_word_length': 5.055753445872696, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Time-Calibrated Species Tree Phylogeny of the New World Cycad Genus Zamia L. (Zamiaceae, Cycadales)

Premise of research. The genus Zamia L. (Zamiaceae), consisting of 79 species, is the most species-rich and widely distributed cycad genus in the New World and arguably the most morphologically and ecologically diverse genus in the Cycadales. However, a strong phylogenetic framework for this genus is still lacking. Methodology. We used a multilocus sequence data set of 10 independent loci (nine single-copy nuclear genes [SCNGs] and one plastid) and extensive taxon sampling (ca. 90% of species) to infer phylogenetic relationships within Zamia. We implemented a concatenated matrix analysis with maximum parsimony, a partitioned maximum likelihood (ML) analysis, and a time-calibrated Bayesian species tree-estimation approach. Diversification, historical biogeography, and ancestral character state reconstruction analyses were conducted using the species tree topology that was the most morphologically and geographically congruent. Pivotal results. We inferred a robust phylogenetic tree for the genus with a strong geographic delimitation of clades and found that four morphological characters typically used for diagnostic purposes in the genus exhibit a high degree of homoplasy. The stem group of Zamia was estimated to have originated at 68.28 Ma (95% highest posterior density [HPD] 51.0–84.5) and the crown group at 9.54 Ma (95% HPD 9.0–10.62). The majority of species richness in the genus appeared during the Pliocene and Pleistocene, with the highest diversification rates occurring in clades comprising Caribbean and South American species. Biogeographic analysis suggests a Caribbean or Mesoamerican origin for Zamia with subsequent dispersal to the Central American Isthmus and South America, where the genus reaches its maximum species and morphological diversity. Conclusions. The high degree of convergence found in multiple morphological characters has historically confounded efforts to elucidate species relationships using nonphylogenetic methods. This study presents the most species-comprehensive, well-resolved hypothesis of phylogenetic relationships within Zamia and provides a strong phylogenetic framework for further studies.


Introduction
The genus Zamia L. (Zamiaceae, Cycadales) is widely considered to be the most ecologically and morphologically diverse genus of the extant cycads (Norstog and Nichols 1997). With its 79 accepted species (Calonje et al. 2019), it is the most species-rich and broadly distributed genus among the New World genera of the Zamiaceae. It is restricted primarily to the Neotropical region (sensu Sclater 1858; Morrone 2014) with only the northernmost Zamia integrifolia L.f. populations in Florida and southeast Georgia extending into the Nearctic region (sensu Escalante et al. 2013;Morrone 2014). The distribution of the genus can be spatially divided into three separate areas of endemism: (1) a Caribbean group restricted to islands on the Bahama Archipelago and Greater Antilles, as well as to the mainland states of Florida and Georgia in the southeastern United States, (2) a Mesoamerican group extending from Tamaulipas, Mexico, to northern El Salvador, and (3) a Central and South American group occurring from southern Nicaragua to Bolivia, the southernmost extent of the genus ( fig. 1).
Despite the remarkable species diversity and broad geographic distribution of Zamia, the only formal infrageneric classification of the genus was published decades ago in Schuster's (1932) monograph of cycads. Schuster included 26 species in his treatment of Zamia and divided the genus into three geographically defined sections (Caribaeae, Mexicano-Meridionales, and Centrali-Meridionales) that are broadly congruent with the Caribbean, Mesoamerican, and Central and South American species groups outlined above. However, Schuster did not designate types for these sectional names, making their taxonomic applica-tion uncertain. Furthermore, several of Schuster's proposed sectional names are illegitimate because they included the type of the name of the genus and did not repeat the generic name unaltered as its epithet, as required by nomenclatural rules (art. 22.1, Turland et al. 2018). Such is the case with Schuster's "Caribbaeae," which is illegitimate because it includes the type of the genus (Zamia pumila L.) and had to be called section Zamia to be considered a legitimate name. Schuster's work is also notorious for contravening the rules of priority and other nomenclatural rules, and it includes very elaborate and often nonsensical infraspecific hierarchies (see Johnson 1959;Hill 1996). Consequently, Schuster's classification has been disregarded by most authors, although some of Schuster's sectional names have been typified and adopted in formal classification work, most notably in Cycas (see Hill 1995). However, no authors since Schuster have attempted formal classification of supraspecific taxa within Zamia. Subsequent authors have instead mostly limited themselves to informally classifying and surmising interspecific relationships within subsets of the genus based on similarities in geographic distribution (Norstog and Nichols 1997), morphology (e.g., Schutzman et al. 1988;Taylor et al. 2008;Calonje et al. 2010Calonje et al. , 2011, anatomy (Acuña-Castillo and Marín-Méndez 2013), karyology (Caputo et al. 1996), genome sizes (Zonneveld and Lindström 2016), and other characters. Nevertheless, the phylogenetic relationships within the genus remain unclear, as only a limited number of studies have involved a species-level phylogenetic component.
Most notably, Caputo et al. (2004) conducted a phylogenetic analysis of Zamia, including 23 of the 79 currently accepted species in a parsimony-based phylogenetic analysis of sequences of the internal transcribed spacer 2 of nuclear ribosomal DNA in combination with a morphological data set. The analysis found several clades to be more congruent with geographic distribution than with morphological similarities, leading the authors to suggest that convergent evolution of morphological characters is pervasive in the genus. Geographically congruent clades identified included those composed of Central American, North American, and Caribbean species, as well as a combined clade including South American and Central American species ( fig. 2A). Clugston et al. (2016) examined the phenological phases of 11 species of Zamia in a phylogenetic context, presenting a tree based on a maximum parsimony (MP) analysis using nucleotide sequence data from two plastid genes (matK, rbcL), two SCNGs (CAB, NEEDLY), and one high-copy nuclear gene (26S). This analysis showed strong support for a Caribbean clade and a monophyletic South American group, and moderate support for a Central American lineage ( fig. 2B).
To examine the diversification ages of extant cycads, Nagalingum et al. (2011) used the nuclear gene phytochrome P (PHYP) with a few additional sequences from the chloroplast genes matK and rbcL to produce a time-calibrated phylogeny of 199 taxa, including 31 species of Zamia. The results indicated that most living cycad species are the result of quite recent diversification of a very ancient lineage and are therefore much younger than previously thought. This was suggested previously by Treutlein and Wink (2002) in an earlier phylogenetic study using rbcL plastid sequences and also found to apply generally to extant gymnosperms by Crisp and Cook (2011). Although infrageneric relationships were not discussed by Nagalingum et al. (2011), their phylogeny supported the monophyly of Zamia, albeit with low node support throughout most of the genus ( fig. 2C). Similarly, in an article using the same molecular markers but with the inclusion of additional species and fossil calibration constraints, Condamine et al. (2015) also found poor node support within genera. Although both divergence time estimation studies included a large number of cycad species in their analyses, they similarly obtained generally low node support within genera, which was likely due to the sparse and often nonoverlapping locus sampling between species and the limited informative sites present in the markers selected.
As illustrated above, the phylogenies published to date that focus on Zamia or contain a considerable sampling of the ge-nus are of limited utility due to low taxon sampling (e.g., Clugston et al. 2016) and/or weak phylogenetic information resulting from a small number of molecular markers (e.g., Caputo et al. 2004;Nagalingum et al. 2011). In this article, we use a larger multilocus sequence data set of 10 independent loci (nine SCNGs and one plastid) and extensive taxon sampling (190% of species) to infer phylogenetic relationships within the genus Zamia using MP, ML, and Bayesian timecalibrated species-tree inference methods. We discuss our results in the context of clade age estimation, diversification, ancestral character state reconstruction, and biogeographical history.

Taxonomic Sampling
DNA was isolated from 113 individual leaflet samples of Zamia belonging to 70 of the 79 currently accepted species in the genus ) and from single samples of Microcycas calocoma (Miq.) A.DC. and Stangeria eriopus (Kunze) Baill. that were included as outgroups for phylogenetic analyses. Multiple samples were collected for some species of Zamia, particularly from those with broader geographic distributions. The leaflet samples were obtained primarily from plants of known provenance or pedigree cultivated in botanic gardens and other botanical collections (app. A).

Outgroup Selection
Microcycas calocoma and Stangeria eriopus, both belonging to monotypic genera, were selected as outgroups for phylogenetic analyses, with the latter as the functional outgroup. Microcycas is strongly supported as the sister genus to Zamia and Stangeria as sister to Zamia 1 Microcycas in a recent genuslevel phylogeny of the Cycadales that used multiple phylogenetic methods (Salas-Leiva et al. 2013) as well as in other phylogenetic studies (e.g., Rai et al. 2003;Bogler and Francisco-Ortega 2004;Chaw et al. 2005;Zgurski et al. 2008).

Marker Selection
Sequences of nine independent, cycad-specific SCNGs and one chloroplast gene (table 1) were used to infer the phylogeny of the genus Zamia using various phylogenetic methods. Of the markers used, seven (40S, ATG2, GroES, HTS, LiSH, PEX4, PMP22, WRKY4) were identified and developed in our laboratories (Salas-Leiva et al. 2014); the remaining ones (CyAG, psbK/I) were derived from previous publications (table 2). The loci 40S, ATG2, CyAG, and GroES were previously used to infer the generic relationships of the Cycadales (Salas-Leiva et al. 2013), and 40S, ATG2, CyAG, GroES, LiSH, PEX4, and WRKY4 were applied to assess the genetic diversity and genetic structure of Bahamian zamias (Salas-Leiva et al. 2017). The single chloroplast gene locus included in the study (psbK/I) was selected because it had displayed the best performance (amplification success and high no. diagnostic sites) in a previous work that evaluated seven different chloroplast genes for barcoding Mexican species of Zamia (Nicolalde-Morejón et al. 2011). 288 Fig. 2 Phylogenetic tree topologies for Zamia resolved by other authors. Some taxon names are revised based on current taxonomy and nomenclature; original names are shown in brackets. A, Caputo et al. (2004): Maximum parsimony consensus tree using combined morphological data and nrDNA ITS2; numbers above branches indicate branch support (up to five-step-longer trees). B, Clugston et al. (2016): Maximum parsimony strict consensus tree found using combined nuclear (26S, CAB, NEEDLY) and plastid (matK, rbcL) genes; numbers above branches are jackknife support values. C, Nagalingum et al. (2011): Maximum likelihood tree generated from concatenated matrix combining one nuclear (PHYP) and two plastid (matK, rbcL) genes; numbers above branches are bootstrap support values above 70%.

DNA Extraction, Amplification, and Sequencing
DNA was extracted from 20 to 100 mg of fresh or silica-geldesiccated leaflet tissue using the FastDNA spin kit (MP Biomedicals, Santa Ana, CA) and the FastPrep-24 instrument (Qbiogene, CA). The extracted DNA was quantified using a GeneQuant Pro RNA/DNA calculator-spectrophotometer (GE Healthcare Life Sciences) or NanoDrop 2000 spectrophotometer (Thermo Scientific) and diluted to 10 ng/mL. Amplifications contained 1# amplification Buffer with 2 mM MgSO 4 , 10 mM of dNTPs, 0.2 mg/mL bovine serum albumin, 10 mM of forward and reverse primer, 0.05 U/mL Taq DNA polymerase (New England BioLabs, Ipswich, MA), and 10 ng/mL of template DNA brought to a total volume of 15 mL with nuclease-free H 2 0. DNA samples were amplified using C1000 and S1000 thermal cyclers (Bio-Rad Laboratories, Hercules, CA). Most genes were amplified in single fragments, with the exception of CyAG and GroES, which were each amplified in two separate fragments using separate primer pairs. For all loci except GroES, the following thermal profile was used: 957C for 2 min, 35 cycles of 957C for 30 s, annealing temperature (507-607C) for 1 min, 727C for 1 min, and a final extension at 727C for 7 min. Both fragments of the locus GroES were amplified using touchdown PCR (957C for 2 min; three cycles of 957C for 30 s, 557C for 1 min, 727C for 1 min; three cycles of 957C for 30 s, 547C for 1 min, 727C for 1 min; three cycles of 957C for 30 s, 537C for 1 min, 727C for 1 min; three cycles of 957C for 30 s, 527C for 1 min, 727C for 1 min; three cycles of 957C for 30 s, 527C for 1 min, 727C for 1 min; three cycles of 957C for 30 s, 517C for 1 min, 727C for 1 min; 30 cycles of 957C for 30 s, 507C for 1 min, 727C for 1 min; and final extension of 727C for 10 min).
The PCR amplifications were evaluated by electrophoresis using 1.2% agarose gel stained with GelRed (Biotium, Hayward, CA) and a size standard ladder (100 bp; New England Biolabs). PCR products were purified using exonuclease I (New England Biolabs) and shrimp alkaline phosphatase (USB Products-Affymetrix, Santa Clara, CA), incubating at 377C for 1 h, followed by 807C for 20 min. Single sequencing reactions used 1-2 mL of purified PCR product and were performed using the ABI Big Dye Terminator 3.1 cycle sequencing kit (Applied Biosystems, Carlsbad, CA), followed by ethanol cleanup. Labeled fragments were visualized on an ABI 3730 automatic DNA sequencer (Applied Biosystems), and the nucleotide sequences were manually edited with Sequencher 4.9 (Gene Codes Ann Arbor, MI). Sequence data were obtained for all samples across all loci except for S. eriopus, which failed to amplify with the PCR markers used for the 40S and LiSH loci. All DNA sequence data (1148 sequences) in this study were deposited in GenBank (app. A).

Sequence Alignment
Sequences for each locus were aligned using MAFFT 7 (Katoh and Standley 2013) and/or aligned manually with Sequencher 4.9. The aligned lengths ranged from 498 bp (PEX4) to 1991 bp (GroES). Sequence data for different loci were joined using Mesquite 3.2 for concatenated analyses (Maddison and Maddison 2017 Taylor & Holzman), necessitating the manual realignment of the inverted regions. The inversions could be dealt with by either reverse complementing the inverted sequence region or inserting a 180 base pair gap in the alignment for all sequences that lacked the inversion (Morrison 2009). We chose the latter, as it was not clear whether the single-nucleotide polymorphisms (SNPs) observed resulted from mutations prior to or after the inversion event.

Phylogenetic Analyses
Maximum parsimony and maximum likelihood analyses. Aligned, concatenated sequences for the 10 loci were analyzed using MP performed in PAUP 4.10b (Swofford 2003) for each gene separately and for the concatenated matrix. Heuristic searches were conducted using 1000 random stepwise addition replicates, with the MultTrees option in effect, and saving up to 10 minimum length trees per search for swapping using tree bisection reconnection (TBR) branch swapping. Jackknife analysis was conducted for the concatenated matrix with 37% deletion probability for each character (JK; Farris et al. 1996; 1000   Lahaye et al. 2008. replicates with simple stepwise-addition, TBR branch swapping, saving no more than 10 trees per replicate, and retaining only groups with frequencies 150%). A partitioned ML analysis was conducted using IQTREE (Chernomor et al. 2016) as implemented on the CIPRES Science Gateway (Miller et al. 2011) using the automated method selection with ModelFinder (Kalyaanamoorthy et al. 2017), and support was calculated by running 1000 bootstrap replicates. Time-calibrated species tree. A time-calibrated species tree analysis was performed using the multispecies coalescent model of StarBeast (Heled and Drummond 2010) as implemented in BEAST (ver. 2.4.4;Bouckaert et al. 2014). Samples were assigned to species according to individual taxonomic assessments, resulting in a species tree topology with 77 terminals.
Separate, uncorrelated, log-normal, relaxed-clock models (UCLD) and random starting trees were assigned to each partition, and a birth-death process was applied as branching process prior (or tree prior), as it has shown to be a better fit for cycads than the traditionally used Yule process (Condamine et al. 2015). The following priors were set following Condamine et al. (2015): a uniform prior between 0 and 10 with a starting value of 0.1 for the mean growth rate, a uniform prior between 0 and 1 with a starting value of 0.5 for the relative death rate, an exponential prior with a mean of 0.33 on the standard deviation of the UCLD model, and a uniform prior between 0 and 1 on the mean of the UCLD model. Site models were estimated independently for each partition (table 3) using a Bayesian approach as implemented by the bModelTest package (Bouckaert and Drummond 2017) available in BEAST.
The species tree was calibrated using 95% HPD age estimates obtained from a previously published time divergence analysis of the Cycadales based on six fossil calibrations and the birth-death tree prior (Condamine et al. 2015), as well as from a recently discovered Zamia fossil (Erdei et al. 2018). Prior distributions for all calibrated nodes were conservatively set to uniform using the minimum and maximum age bounds outlined below. The age intervals for the tree root node (74.3-147.7 Ma) and the crown node of Zamia (9-22.1 Ma), as well as the maximum age bound (84.5 Ma) for the stem node of Zamia, were set according to the 95% HPD estimates provided by Condamine et al. (2015). The minimum age bound for the stem group of Zamia was constrained to the minimum age estimate (33 Ma) of a fossil cycad leaflet of unknown terrestrial origin recently found preserved in marine sediments of the Gatuncillo Formation in Central Panama, which is as-signable to Zamia based on cuticular micromorphology (Erdei et al. 2018). The fossil age is considered to be 35-33 Ma based on nannoplankton and foraminiferal biostratigraphy, and the minimum age estimate (33 Ma) was used to constrain the minimum age bound for the stem group of Zamia. This age is comparable to the lower bound of the 95% HPD estimate (34.2 Ma) obtained by Condamine et al. (2015).
Two independent runs of 750 million Markov chain Monte Carlo (MCMC) iterations were conducted, sampling every 2500 iterations and preceded by a 100 million iteration burnin period. Log and tree files were combined using Logcombiner (ver. 2.4.4, included in BEAST), resulting in 600,000 trees. The log output was evaluated using Tracer (ver. 1.6; Rambaut et al. 2014), and a maximum clade credibility (MCC) tree was created from these trees using TreeAnnotator (ver. 2.4.4, included in BEAST) and visualized using Figtree (ver. 1.4.3;Rambaut 2012). The MCC tree topology recovered was used for the diversification, biogeographic, and ancestral character state reconstruction analyses presented here.

Diversification Analyses
Diversification rates were estimated per-million-year interval spanning 10-0 Ma for the crown node of Zamia as well as several major clades within the genus using the dates derived from the time-calibrated species tree. The rates were obtained using Foote's (Foote 2000) equation for per-capita origination rates as modified to apply to molecular tree lineages by Nagalingum et al. (2011). We examined changes in diversification rates for each of the major clades using the APE package for R (Paradis et al. 2004) to calculate the g statistic (Pybus and Harvey 2000) from the branching times in the time-calibrated species tree. The g statistic is a metric of the distribution of speciation times in molecular phylogenies, with negative values indicating early diversification and positive values indicating late branching. A Monte Carlo constant rates (MCCR) test was conducted to determine whether the distribution of diversification rates significantly differed from a constant rate (P < 0:05). The MCCR test simulates a distribution for the g statistic by simulating trees that take into account the incomplete lineage sampling of the observed tree and providing a critical value for rejecting a constant diversification rate. The MCCR test was conducted using the LASER package in R (Rabosky 2006) by simulating 10,000 trees under a pure birth model. Lineage through time (LTT) plots were prepared using the time-calibrated species tree using the APE package in R. We examined the stem and crown nodes of Zamia as well as several major clades within the genus.
Net diversification rates were calculated according to Magallon and Sanderson (2001) using the package Geiger (Harmon et al. 2008) as implemented in R. The method requires the input of the extinction rate as a fraction of the speciation rate (e), but since the extinction rate cannot be reliably estimated, we followed previous recommendations (Magallon and Sanderson 2001;Crisp and Cook 2011) to provide net diversification rate estimates based on several values ranging from low (e p 0) to high (e p 0:9).
The Comet model ) from the R package TESS ) was used to test for shifts in speciation or extinction rates as well as for the signature of mass extinction in the crown node of Zamia. We selected a single prior for  Hodell et al. 2001;Herbert et al. 2016). The survival rate prior was set to 0.05, representing a major extinction event.
We used BAMM 2.5.0 (Rabosky 2014) to test evolutionary rate heterogeneity in the topology of the crown node of Zamia using the default priors and settings recommended by the developers.

Ancestral Character State Reconstruction
We used Mesquite 3.02 (Maddison and Maddison 2017) to study the evolution of four morphological characters on the calibrated species tree by performing ancestral character state reconstruction using parsimony. The morphological characters used were presence of prickles on the petiole, prominence of teeth along leaflet margin, prominence of leaflet veins, and arborescence (table 4). These characters were selected because they are considered informative and are commonly used to diagnose different species within the genus.
Zonneveld and Lindström (2016) examined genome sizes of most Zamia species and inferred a Mexican center of origin for Zamia with a migration south to South America based on increasing DNA content in the southernmost species. Furthermore, they proposed Mega-Mexico, Caribbean Island, and South American biogeographic groupings based on genome size, with several Central American species falling into different groupings based on purported ancient hybridizations. We examined these inferences by examining the evolution of genome size using the species tree topology. For this purpose, we used the function contMap from the phytools R package (Revell 2012) to map the haploid genome sizes (1C) compiled by Zonneveld and Lindström (2016) onto a pruned phylogeny for the crown group of Zamia. We used the function phylosig to compute the phylogenetic signal using Pagel's l. This statistic reflects the similarity of the observed trait distribution to the result of a passive Brownian motion process, with a perfect fit indicated by a lambda of 1. We conducted this analysis to test whether genome size passively reflects the relationships between species and clades instead of being related to a north to south increase.

Biogeographic Analyses
The historical biogeography of Zamia was inferred using a statistical dispersal-vicariance analysis (S-DIVA; Yu et al. 2010) as implemented in RASP 3.2 (Yu et al. 2015). The MCC tree resulting from the time-calibrated species tree analysis was selected as the condensed tree, and the 600,000 individual postburn-in trees obtained from the analysis were resampled to a more manageable 14,000 trees using Logcombiner 2.4.4. Of these, 100 were randomly selected in RASP for the analysis. Areas coded for the analysis include the following: Ap Caribbean, B p Mesoamerica, C p Isthmus, D p South America, and E p South Africa. The maximum number of ancestral areas at each node was constrained to four, and the estimated probabilities of the ancestral areas were visualized on the condensed tree.
We also reconstructed the historical biogeography of Zamia using BioGeoBEARS (Matzke 2013). Analyses were conducted using the dispersal-extinction-cladogenesis (DEC) and DIVAlike models, as well as variants of these models with the jump parameter (DEC1J and DIVA1J). All analyses were nonstratified and with equal probability of transition between areas. Ancestral area number occupancy was limited to two areas at the time.

Data Availability
The following data used in or resulting from the analyses in this article have been uploaded and are available on the FigShare public file repository at https://figshare.com/s/34b73962882c8 adb0c73. They include (1) the partitioned sequence alignment used for all phylogenetic analyses, (2) the results of the Bio-GeoBEARS biogeographical analyses, and (3) the following phylogenetic trees: time-calibrated MCC species trees as well as individual locus trees resulting from the StarBEAST analysis, trees resulting on concatenated matrix and partitioned ML analyses, and trees based on MP analyses of concatenated matrix and individual loci.

Single Locus and Concatenated Matrix Sequence Characteristics
Single locus MP analyses generally yielded poorly resolved consensus trees (see data availability), indicating that the incongruence between them is most likely soft incongruence because of the weak phylogenetic information of most individual gene trees, rather than hard incongruence, which is due to significantly conflicting topologies (Wendel and Doyle 1998). The range of parsimony informative characters (table 5) ranged from 29 (PEX4) to 127 (GroES), or by percentage, from 5.1% (LiSH) to 10.2% (ATG2). Consistency indexes (CI) were above 0.82, and retention indexes (RI) were above 0.89 for all examined to-

Major Clades Recovered
The MP concatenated analysis, ML partitioned analysis, and Bayesian inference (BI) species tree approach recovered the genus Zamia as monophyletic and sister to Microcycas. All analyses recovered the same broad topology consisting of the following geographically delimited major clades ( fig. 3): (1) Caribbean clade, consisting mostly of Caribbean Island species (a single species, Zamia integrifolia, also reaches Florida), which is sister to the rest of the genus consisting of primarily of mainland American species; (2) Fischeri clade, consisting of three Mexican endemic species and is itself sister to the rest of the genus excluding the Caribbean clade; (3) Mesoamerica clade, including all other species occurring in Mesoamerica to the exclusion of the Fischeri clade and Zamia soconuscensis; (4) Isthmus clade, consisting primarily of Panamanian and Costa Rican species; and (5) South America clade, consisting primarily of South American endemic species. The latter two clades are sister to each other and together form the most species-rich and morphologically diverse clade in the genus, referred to here as the Central-meridional clade. All the above clades were strongly supported by all analyses except for the Mesoamerica clade, which was recovered by all analyses, albeit with poor support (MP p 57.8%, ML p 50%, BI p 0.63).

Maximum Parsimony, Maximum Likelihood, and Species Tree Inference Analyses
Although all phylogenetic methods resolved the same broad topologies and recovered the same major clades ( fig. 3), the species tree analysis recovered several relationships more congruent with morphological similarities and/or geographic proximity than the other two analyses and will therefore be used in this article to discuss phylogenetic relationships, node age estimates, and biogeographic patterns within the genus Zamia. The MP

Time-Calibrated Species Tree
Single-locus MCC trees produced by StarBeast analysis generally yielded poorly resolved trees (see data availability) compared with the species tree topology ( fig. 4). The effective sample size (ESS) scores from the time-calibrated species tree analysis were all above 200, indicating good levels of convergence of the MCMC. The mean stem age of Zamia was calculated at 68.28 Ma (95% HPD 51.0-84.5) and the crown age at 9.54 Ma (95% HPD 9.0-10.62; fig. 5; table 6). The crown node that includes the Mainland Zamia clade is the oldest within the genus, with a mean age of 5.97 Ma (95% HPD 3.9-8.2). The Caribbean clade, despite being sister to the rest of the genus, has the youngest crown age-1.19 Ma (95% HPD 0.61-1.86)-of the major clades. The Fischeri clade, despite consisting of only three species, is relatively old, with a crown age of 2.59 Ma (95% HPD 0.75-4.62). The Mesoamerica clade, with a crown age of 3.9 Ma (95% HPD 2.54-5.34) is of comparable age to the Central-meridional clade, which is 3.8 Ma (95% HPD 2.49-5.17). Within the Central-meridional clade, the South America clade, with a crown age clade of 2.62 Ma (95% HPD 1.71-3.56), is slightly older than the Isthmus clade, which is 2.35 Ma (95% HPD 1.43-3.32).

Diversification Analyses
Net diversification rates. Net diversification rates were highest when lower extinction rate (e) values were used in the calculations (table 8). The net diversification rate for the crown of Zamia ranged from 0.22 to 0.39. Among the major clades examined, using e values of 0 and 0.7, the Caribbean clade had the highest diversification rates (1.26, 1.26, 0.88), followed by the South America clade (1.04, 1.03, 0.81) and the Isthmus clade (0.88, 0.65, 0.37). The lowest diversification rates were found in the species-poor Fischeri clade (0.16, 0.10, 0.05), followed by the Mesoamerica clade (0.60, 0.46, 0.27). In the analysis  3 Broadly congruent tree topology obtained from maximum parsimony (MP; PAUP), maximum likelihood (ML; RAXML), and Bayesian species tree (BEAST) analyses (BA) of nine single-copy nuclear gene loci and one plastid locus. All methods resolved the same major clades and relationships between them. Underlined numbers to the right of nodes are posterior probabilities from the BA, numbers below branches in italics are bootstrap support values for ML, and numbers above branches in boldface are jackknife support values for the MP analysis. Numbers in parentheses next to clade names represent the number of taxa included in these clades. using ε p 0:9, the diversification rate was slightly higher for the South America clade than the Caribbean clade. However, ε values of 0.9 are exceptionally high (Magallon and Sanderson 2001) and may be unrealistic for such a short time frame, so the results using ε p 0:9 should be interpreted with caution.
Per-million-year diversification and diversification rate heterogeneity. Most diversification above the crown node of Zamia occurred in the Pliocene and Pleistocene epochs (figs. 6, D1), with only the split between the Caribbean and Mainland clade and between the Fischeri clade and the rest of the Mainland Fig. 4 Maximum clade credibility species tree cladogram obtained from analysis in StarBeast using nine single-copy nuclear genes and one plastid locus. Posterior probabilities displayed above branches. Circles at nodes are labeled with letter codes identifying clade names with major clades color coded. Inset map shows the geographic provenance of the samples used in this study, color coded by major clade membership. clade occurring in the Miocene. The highest per-million-year diversification rates were found in the Caribbean clade followed by the South America clade ( fig. 6; table 9), which were the two clades that also had the highest net diversification rates. Among the major Zamia clades, the MCCR tests only found diversification rates differing significantly from constant (P < 0:05) in the South America clade. In this case, the negative g statistic indicated rapid, early diversification of the South America clade followed by a decrease in diversification rate over time. Similarly, the BAMM analysis did not uncover shifts in diversification rates across the tree. Finally, the TESS analysis showed no evidence of shifts in extinction or diversification rates and only marginal evidence of mass extinctions at around 6 and 2.5 Ma ( fig. D2).

Ancestral Character State Reconstruction
Morphological characters. All four morphological traits examined exhibited homoplasy, and the ancestral states recovered include subterranean stems, smooth petioles, entire leaflet margins, and leaflets without prominent veins (figs. D1-D4). Arborescence ( fig. D1) is absent in the Caribbean clade and rare in the Mesoamerica clade, where it is restricted to only three species. The ancestral state for the Central-meridional clade is ambiguous, leading to an ancestral state of arborescence for the Isthmus clade and subterranean stems for the South Amer-ica clade. Arborescence is most prevalent in the Isthmus clade, with only two species possessing subterranean stems. Within the South America clade, arborescence seems to have evolved separately four times. Regarding petiole armature, unarmed petioles are universal within the Caribbean clade but extremely rare among mainland species (fig. D2). Both leaflet margin dentation ( fig. D3) and leaf venation prominence ( fig. D4) are extremely homoplastic characters that appear to have evolved independently multiple times within the genus and may be of limited use as diagnostic characters.
Genome sizes in Zamia. We recovered a strong phylogenetic signal on the evolution of genome size within the crown group of Zamia (Pagel's l p 0.9721436), with both genome expansion and contraction occurring within different lineages ( fig. D5).

Biogeographic Analyses
Historical biogeography of Zamia. Optimization of ancestral areas on the species tree using S-DIVA ( fig. 7) required six dispersal and five vicariance events (table 7), with most speciation events occurring within geographical areas. Dispersal events occurred from the Mesoamerica (B) region to the Isthmus (C) and South America (D) regions, from the Isthmus to the South America region, and from South America to the Isthmus region.
The crown node of Zamia 1 Microcycas (node 152, probability of ancestral range at node, P p 0:52) is optimized in two areas, one being the Caribbean (A) region, the other being a composite area including the Caribbean (A) and Mesoamerican (B) regions, with a dispersal event from the composite area to the Caribbean region. The crown node of Zamia (151, P p 1) is optimized in a composite area consisting of the Caribbean (A) and Mesoamerica (B) regions with a vicariance between these two regions separating the Caribbean clade from the rest of the genus on the American mainland. The Caribbean clade is ancestrally optimized in the Caribbean (A) region with no events hypothesized (85, P p 1), and stasis within all inner nodes (78-84, all with P p 1). The crown node of the clade including all mainland zamias (150, P p 0:63) is optimized in the Mesoamerican region (B) with no events hypothesized. The crown node of the Fischeri clade is optimized to the Mesoamerica region (87, P p 1) with no events hypothesized. The clade including all Zamia except the Caribbean and Fischeri clade (149, P p 0:21) has a poorly supported optimization among four areas, with the highest support value found for the Mesoamerican region (B, 62.6%). The other three areas were all composite areas that included the Mesoamerica region. Two dispersal events are hypothesized from the Mesoamerica region to the Isthmus and South America regions. The Mesoamerica clade (104, P p 1) is optimized in the Mesoamerican region with high probability (P p 1), no events hypothesized and stasis within all inner nodes (88-103, all with P p 1). The node that includes Zamia soconuscensis as sister to the Central-meridional clade has a poorly supported optimization (148, P p 0:34) consisting of three composite areas, all of which include the Mesoamerican region and vicariance between Mesoamerica and a composite area consisting of the Isthmus and South America regions.
The split of the Central-meridional clade (147, P p 0:98) into the Isthmus and South America clades was attributed to a vicariance event occurring 3.8 Ma (95% HPD 2.49-5.17) Lindstr., Zamia cunaria, and Zamia ipetiensis D.W.Stev., which is optimized to a composite area including the Isthmus (C) and South America (D) regions with one vicariance event between these two regions.
The ancestral areas recovered for nodes using the DIVA and DEC models in BioGeoBEARS and the S-DIVA model in RASP were identical and differed only slightly in their area probabilities. Implementing the jump parameter, a jump between Mexico and South America (in the DEC1J) or Mexico and Central America (in the DIVA1J) at the base of the Central-meridional clade is favored. However, given the issue that recently emerged with the jump parameter estimation (Ree and Sanmartín 2018), the results of the DEC1J and DIVA1J must be viewed with caution. As the BioGeoBEARS analyses recovered the same ancestral areas as the S-DIVA analysis, albeit with slightly different area probabilities, we will limit our discussion of the historical biogeography of Zamia to the results of the RASP S-DIVA analysis. The BioGeoBEARS analyses are available in the FigShare file repository (see data availability).

Utility of psbK/I for Barcoding and Phylogenetic Analyses
Despite the high rate of amplification success across all cycad genera and a high number of diagnostic sites for DNA barcoding demonstrated for psbK/I (Nicolalde-Morejón et al. 2011), its utility for barcoding and phylogenetic inference may be limited by the fact that sizable inversions are present in several taxa. If the inversions are not detected or dealt with during multiple sequence alignment, they could prove problematic because they would result in multiple false-positive SNPs. For example, the unusually high number of plant barcoding diagnostic sites reported for Zamia pseudoparasitica by (Nicolalde-Morejón et al. 2010) appears to be due to an inversion present in this species that was neither reverse complemented nor staggered during the alignment process. A comparable situation has also been found in the chloroplast barcoding candidate locus trnH-psbA, where frequent inversions in multiple angiosperm lineages present a challenge to the gene's utility as a plant barcoding region (Whitlock et al. 2010).

Phylogenetic Relationships and Major Clades within Zamia
Below we consider the phylogenetic relationships within Zamia, providing brief descriptions of the major clades recovered in our analyses as well as for some notable subclades that we define for the first time below. Informal names are provided for these monophyletic groups for the purposes of discussion clarity ( fig. 4).
Monophyly of Zamia. The monophyly of Zamia was strongly supported (fig. 4, node A; BI/MP/ML p 1/100/100) Fig. 6 Per-million-year diversification rates and g values for major Zamia clades. Only the g value for the South America clade (asterisk) was significant, indicating that diversification was initially high and that it decreased over time. The g value for all other clades was not significant, indicating a constant diversification rate. and its sister relationship to Microcycas was recovered in all analyses. Zamia consists of two sister clades: the Caribbean clade, consisting primarily of Caribbean Island species, and the Mainland clade, which largely consists of mainland American species.
Caribbean clade. The Caribbean clade ( fig. 4, node A; BI/ MP/ML p 1/100/100) consists of nine species with subterranean stems, unarmed petioles, and the same 2np16 chromosome number (Olson and Gorelick 2011). The group is primar- Fig. 7 Historical biogeography of Zamia inferred using a statistical dispersal-vicariance analysis as implemented in RASP using the maximum clade credibility chronogram tree shown in figure 4. The proportions of colors in circles represent possible ancestral geographic ranges at each node, and the rings around the circles represent dispersal and/or vicariance events. The two basalmost nodes are out of scale; their mean ages are indicated next to the asterisks. Inset map shows the geographic distribution of the Zamia samples color coded by their area assignations. 300 ily distributed on several islands in the Greater Antilles and Bahama Archipelago, with a single mainland species restricted mostly to Florida, albeit with a few historical collections in extreme southeast Georgia (Duncan 1979). Caribbean zamias are universally accepted as a distinct lineage (e.g., Eckenwalder 1980a; Stevenson 1987;Sabato 1990), yet their current classification (e.g., Stevenson 1987;González-Géigel 2003;Osborne et al. 2012) remains controversial. This classification is primarily based on a historic overreliance on vegetative morphological characters, especially leaflet macromorphology, for species delimitations. Consequently, populations exhibiting similar leaflet morphologies are currently considered conspecific (González-Géigel 2003) despite sometimes occurring on geographically distant islands that in many cases were never connected by land. This scenario would require multiple occurrences of longdistance dispersal over the Caribbean Sea. An alternative explanation is that some of these morphotypes may have evolved independently multiple times within the clade. Recent genetic studies have begun to address this question and provided further insight into the genetic relationships within this group (Meerow 2007;Meerow et al. 2012;Salas-Leiva et al. 2017;Meerow et al. 2018). Ongoing systematics research should help clarify the currently controversial taxonomy of Caribbean zamias in the near future.
Mainland clade. The Mainland clade ( fig. 4, node B, BI/ MP/ML p 1/100/100) is sister to the Caribbean clade and consists chiefly of species occurring on the American mainland, with only a few species also occurring on nearshore continental islands (e.g., Zamia hamannii A.S.Taylor, J.L.Haynes & Holzman). It is the most species-rich group in Zamia (90% of species) and has a broad geographic range that extends from Mexico to Bolivia. It consists of multiple geographically defined clades discussed below.
Fischeri clade. The Fischeri clade ( fig. 4, node C, BI/MP/ ML p 1/100/100) is sister to all remaining mainland species and consists of three species occurring in northeastern Mexico. Zamia fischeri Miq. and Zamia vazquezii D.W.Stev., Sabato & De Luca are morphologically very similar, with the latter segregated from the former based on differences in leaf size, leaflet shape and texture, and differing chromosome counts (Stevenson et al. 1998). Our analyses confirm that the two species form a monophyletic group, and we report for the first time this group's sister relationship to Zamia inermis Vovides, J.D.Rees & Vázq.Torres. All three species occur in relative geographic proximity and have unarmed or sparsely armed petioles and stems that are subterranean or very short. In addition, Z. inermis shares the same chromosome count (2np16) as Z. fischeri (Vovides 1983;Stevenson et al. 1998) and has a reproductive phenology similar to that of Z. vazquezii (Griffith et al. 2012).
Mesoamerica clade. This group ( fig. 4, node D) of 18 species is restricted to the Mesoamerican dominion (sensu Morrone 2014), which it shares with Zamia soconuscensis and the three Fischeri clade species. The clade was recovered in all analyses, albeit with low support (below 70%). It is divided into the sister Tuerckheimii and Mexico subclades. The Tuerckheimii subclade ( fig. 4, node E), recovered only in the species tree analysis with high support (BI p 0.96), is a group consisting of six species endemic to Guatemala, Honduras, and Belize. The Mexico subclade ( fig. 4, node B; BI/MP/ML p 1/98.4/100) consists of mostly Mexican endemic species with subterranean stems, armed petioles, and dentate leaflet margins. This group is further divided into sister Purpurea and Furfuracea subclades. Although the Purpurea subclade (fig. 4, node G) was only poorly supported in the species tree analysis (BI p 0.64), its constituent species have been recognized as morphologically distinct by several authors (e.g., Schutzman 1984;Schutzman and Vovides 1998;Pérez-Farrera et al. 2012), albeit always with the inclusion of Zamia standleyi Schutzman, which in our results was resolved as deeply embedded within the Tuerckheimii subclade ( fig. 4,   Phylogenetic placement of Zamia soconuscensis. Surprisingly, the Mexican species Zamia soconuscensis was recovered as sister to all species from the Central American Isthmus region and South America and does not appear to be closely related to other Mesoamerican species. This inclusive clade was strongly supported in the MP and ML analyses but poorly supported in the species tree analysis (BI/MP/ML p 0.4/98/100). The Soconusco region in Chiapas where Z. soconuscensis occurs is believed to have been a primary Pleistocene floristic refuge (Toledo 1982), so this species may be the sole survivor of a once larger group of species eliminated from surrounding areas during Pleistocene glaciations. Alternatively, the species may have other close relatives in neighboring Guatemala, as on-site observations point to several as yet undescribed species of Zamia in the region. The potential discovery of related species and/or the inclusion of additional loci may aid in phylogenetic reconstructions to explain or reject the unusual placement recovered for Z. soconuscensis in our analyses.
Central-meridional clade. The Central-meridional clade ( fig. 4, node I; BI/MP/ML p 0.98/99.6/98) includes more than 60% of all species in the genus, making it by far the most speciesrich clade within the mainland zamias. It includes all species occurring from southern Nicaragua south to the southern limit of the genus's range in Brazil and Bolivia. It is separated by a large distribution gap from Mesoamerican Zamia populations, which only extend to northern Honduras and El Salvador. This group consists of two sister clades: one that is mainly species endemic to the Central American Isthmus (Isthmus clade) and another that is mainly South American species (South America clade).
Isthmus clade. The Isthmus clade ( fig. 4, node J; BI/MP/ ML p 1/99.9/100) consists of 15 species mostly restricted to the Central American Isthmus region, with only Zamia obliqua extending into Colombia in adjacent South America. Most species are arborescent, with only two (Z. nana and Z. dressleri) having subterranean stems. Notable groups within the Isthmus clade include the Acuminata subclade ( fig. 4, node K; BI p 0.92) with four species restricted to the Pacific side of the Cordillera Central mountain range in Costa Rica and Panama, the Obliqua subclade ( fig. 4, node L; BI/ML p 0.88/87), consisting of arborescent species with very sparsely armed petioles, and the Skinneri subclade ( fig. 4, node M; BI/MP/ML p 0.83/100/ 100), composed of five arborescent species with large, broad, prominently veined leaflets and occurring primarily on the Atlantic side of the Cordillera Central from central Panama through southern Nicaragua. Zamia dressleri, a species with similar leaflet morphology but with a subterranean stem, has previously been considered part of the Z. skinneri Warsz. ex A.Dietr. species complex (Taylor et al. 2008), but this relationship was not supported by our analyses which recovered the species as sister to the Obliqua subclade.
South America clade. The South America clade ( fig. 4, node N; BI/MP/ML p 1/100/100) consists of 28 species endemic to South America, with the exception of Zamia manicata Linden ex Regel, which extends into Panama, and sister species Zamia cunaria and Zamia ipetiensis (Cunaria clade; fig. 4, node R; BI p 0.89), which are Panamanian endemics. The phylogenetic relationship of the latter two species to other species occurring in South America remains unclear, as the phylogenetic analyses resulted in conflicting and poorly supported topologies. However, these species share extensive morphological similarities with Zamia pyrophylla from the Colombian Chocó, including subterranean stems holding only one to two leaves, tomentum present on strobili axes, and microsporangia on the adaxial side of the microsporophylls (Calonje et al. 2010).
The Manicata subclade ( fig. 4, node O; BI=MP=ML p 1=99:3= 100) is a group of four species occurring primarily in northern Colombia that have a remarkably variable leaflet morphology including the distinctly channeled petiolule with a gland-like collar of Zamia manicata, the membranaceous and prominently veined leaflets of Zamia disodon D.W.Stev. & Sabato, and the unusual midrib found in Zamia restrepoi (D.W.Stev.) A.Lindstr. The midrib found in Z. restrepoi leaflets is unique in the genus and so distinctive that the species was initially described within its own genus (Chigua D.W.Stev.) until it was eventually subsumed into Zamia (Lindstrom 2009). The strikingly diverse leaflet morphology found in this clade has historically obscured the close phylogenetic relationships of its constituent species, but they do share other characters such as acaulescent stems, toothed leaflets, relatively small seeds, and microsporophylls with a very short fertile section of lamina; these characters support the close relationship recovered in our analyses.
The Pacific subclade ( fig. 4, node P; BI=MP=ML p 0:78=70= 89) consists of five species occurring from the western foothills of the Andes to the Pacific coast. The species are variable in terms of leaflet vein prominence, stem habit, and leaflet margin dentation (figs. D1-D4) but are all extensively armed with prickles on the petioles, making it the most aggressively armed clade in the genus. Panamanian highland species Zamia lindleyi Warsz. Ex A.Dietr. was previously considered conspecific with Zamia chigua Seem. (Stevenson 1993) from the lowlands of the Colombian Chocó biogeographic region primarily due to having fernlike leaves with numerous narrow leaflets. However, several morphological and ecological differences in addition to the large geographic disjunction have since brought recognition of Z. lindleyi as a distinct and separate species from Z. chigua (Calonje et al. 2012b), a recognition supported by our analyses, which place the former within the Isthmus clade and the latter in the South America clade. Norstog (1980) believed that Pacific subclade species Zamia roezlii Linden (as Z. chigua; Norstog 1986) was the most primitive species in the genus because its stable rainforest habitat allowed it to retain what he considered primitive characters such as an arborescent habit, large spermatozoids, and a large asymmetrical karyotype. However, three of the four traits we examined for this species, including the arborescent habit, were not ancestral in our character state reconstructions (figs. D1-D4). Additionally, karyotype evolution in Zamia appears to be moving toward increased asymmetry with higher numbers of smaller chromosomes (Vovides and Olivares 1996; Olson and Gorelick 2011), meaning the large asymmetric karyotype of Z. roezlii is currently considered more derived than ancestral.
The Wallisii subclade ( fig. 4, node Q; BI/MP/ML p 0.78/93/ 99) consists of Zamia wallisii Braun and Zamia oligodonta E. Calderón & D.W.Stev., two subterranean-stemmed montane species with prominently veined, extremely broad, coriaceaous leaflets and ovoid female cones with few, relatively large seeds, 302 both occurring on the western flank of Colombia's Western Cordillera. Zamia montana Braun, an arborescent species that also has prominently veined leaflets and occurs in the same geographic region, was not sampled in this study but likely belongs to this subclade as well. In fact, Z. montana and Z. oligodonta were previously considered conspecific due to vegetative similarities (Lindstrom 2009), but recent field studies with the two species (Calonje et al. 2015) have clarified that the two species are significantly morphologically distinct (Calonje et al. 2015).
The The Eastern subclade ( fig. 4, node T; BI/MP/ML p 0.81/ 64/88) consists of species occurring in eastern Colombia and western Venezuela. All species occur to the east of the Andes mountain range with the exception of Zamia encephalartoides D.W.Stev., which occurs on the western side of the Cordillera Oriental mountain range. The taxon treated here as Zamia lecontei Ducke is from Amazonian Venezuela, following the species circumscription by Stevenson (2004). However, the type locality for this species is a great distance away in Brazil, and additional taxonomic work is required to determine whether Brazilian, Venezuelan, and Colombian populations of this species are truly conspecific.
The Amazonian subclade ( fig. 4, node U; BI/MP/ML p 0.64/ 70/95), recovered with high support only in the ML analysis, consists of seven species endemic to the Amazon basin. Zamia poeppigiana Mart. & Eichler, a large arborescent Amazonian species, shares many morphological similarities with Zamia lindenii Regel ex André, which occurs on the Pacific side of the Andes. The two were previously considered conspecific (Stevenson 2001) until additional taxonomic work helped clarify morphological differences between them (Lindstrom 2009;Calonje et al. 2011). Their placement in the separate Pacific and Amazonian subclades in our analyses justify their recognition as separate species.

Ancestral Character States and Homoplasy of Macromorphological Vegetative Characters
Homoplasy was prevalent within all the macromorphological characters examined (i.e. arborescence, petiole armature, leaflet margin dentation, and leaflet vein prominence), and no single character or combination thereof appeared to be diagnostic for the major clades within Zamia (figs. D1-D4). This is in accordance with previous findings by Caputo et al. (2004), who found all characters from an extensive morphological and micromolecular data set to be homoplastic, as well as with inferences from population-based studies in the Caribbean clade (Salas-Leiva et al. 2017;Meerow et al. 2018). The rampant homoplasy pervasive in the genus has resulted in several species that share morphological resemblance but that are not closely related. This is perhaps best exemplified by species with prominently veined leaflets, a character unique to the genus and so distinctive that it was the primary reason for the segregation of the genus Aulacophyllum Regel from Zamia (Regel 1876). This character state, found in approximately one-quarter of all Zamia species, has evolved multiple times throughout the genus and is not diagnostic of any clade ( fig. D4). Similarly, several large arborescent South American species with prominently toothed leaflet margins (Zamia tolimensis, Z. lindenii, and Zamia poeppigiana), which were previously considered to belong to the same species complex (Calonje et al. 2011), are actually members of three separate lineages within the South America clade. The pervasiveness of homoplasy within macromorphological vegetative characters typically considered informative complicates their diagnostic utility, and future systematic research will need to consider an array of characters, including morphological, anatomical, reproductive, and genetic characters.

Genome Sizes in Zamia
We recovered a strong phylogenetic signal for the evolution of genome size within the crown group of Zamia (Pagel's l p 0.97), indicating a strong statistical dependence between DNA content and our species tree topology. Zonneveld and Lindström (2016) posited a southward migration of Zamia based on increasing genome size to the south. However, genome size appears to have evolved stochastically within the genus, with both increases and decreases in DNA content occurring in different clades. Zonneveld and Lindström (2016) noted patterns between genome size and geographic distribution, identifying three major biogeographic groups matching the Mesoamerican, Caribbean, and South and Central American areas of endemism ( fig. 1) for the genus. However, the latter two groups are not diagnosable based on genome size alone, as they include overlapping values, a situation that Zonneveld and Lindström (2016) attributed to ancient hybridizations between the two groups. Our results instead suggest that these similarities are due to both groups having genome sizes that have remained stable since the diversification of the crown node of the genus ( fig. D5). Similarly, the Mesoamerican group, while diagnosable by its low genome size, in our phylogenetic analyses represents a paraphyletic assembly that includes members of the separate Mesoamerica and Fischeri clades as well as Zamia soconuscensis.

Divergence Age Estimation
The divergence between Zamia and Microcycas was estimated to be 68.28 Ma (95% HPD 50.99-84.5), an estimate older than those presented by Condamine et al. (57 Ma),Salas et al. (36.5 Ma),and Nagalingum et al. (31.1 Ma). Our estimate suggests that the divergence between the two genera happened near the Cretaceous-Paleogene boundary, which was marked by the catastrophic Chicxulub asteroid impact that may have wiped out up to 60%-70% of plant species in tropical North America (Nichols and Johnson 2009). The impact likely had a devastating effect on contemporary cycad flora in the Americas, but the aftermath likely presented a speciation opportunity for surviving cycad lineages. The estimated stem age for Microcycas is older than the present Antillean islands, which were formed less than 40 Ma (i.e., late Eocene; Iturralde-Vinent 2006). This suggests the ancestral range of Microcycas may have orig-CALONJE ET AL.-A TIME-CALIBRATED SPECIES TREE PHYLOGENY OF ZAMIA L.
inally included mainland populations that eventually became extinct. However, such interpretation should be taken cautiously as the fossil record for Zamia is scarce and is, as yet, inexistent for Microcycas. As the latter is a monotypic genus, there are no diversification events to infer that could clarify its biogeographic history.
The mean crown age of Zamia was estimated at 9.54 Ma (95% HPD 9.0-10.6 Ma). This age is younger than reported by Condamine et al. (2015;14.6 Ma, using new fossil data set and birth/death prior), slightly younger than reported by Nagalingum et al. (11.25 Ma, BI using nuclear markers), and slightly older than reported by Salas et al. (8.2 Ma). Our crown age estimate for Zamia marks the divergence of the Caribbean clade from the Mainland clade and is concordant with the timesince-divergence estimate (Sarich 1977) of 10.8 Ma between Zamia pumila L. (sensu Eckenwalder 1980a) and Zamia splendens reported by Walters and Decker-Walters (1991) based on allozyme differentiation. However, our use of a stem calibration for the genus represents a rather conservative approach that might push the age of the crown group toward the present. We thus suggest to interpret our estimates as minima that will be more properly resolved by employing more sophisticated ways of dealing with fossil information.

Biogeography of Zamia
Geographic delimitation of clades. Three separated major areas of endemism are readily identifiable in the geographic distribution of the genus Zamia: the Caribbean Islands and southeast United States, Mesoamerica, and Central and South America ( fig. 1). The latter two groups are separated by a considerable distribution gap extending from northern Salvador/Honduras to southern Nicaragua. The Zamia species occurring in the Caribbean Islands and southeast United States are all monophyletic (Caribbean clade), as are all species from Central and South America (Central-meridional clade). However, in the Mesoamerica region we found two separate clades with overlapping geographic distributions (Fischeri and Mesoamerica clades) as well as the unusual situation with Z. soconuscensis, which in our analyses appears more closely related to the Central-meridional clade than to any Mesoamerican species. With the notable exception of Z. soconuscensis, we found a strong pattern of phylogenetic congruence with geographic distribution in Zamia, with most of the major clades recovered in our analyses exclusively occurring in particular geographic regions with little or no range overlap. Indeed, our results bolster similar findings by Caputo et al. (2004) by demonstrating that this pattern is widespread throughout the entire genus.
The strong congruence between geographic distribution and phylogenetic relationships recovered in our analysis is a pattern consistent with the limited dispersal and pollen-transfer abilities found in cycads (Tang 1993;Yang and Meerow 1996), which most likely limit the ability of these clades to expand beyond geographically contiguous and ecologically suitable conditions following isolation via vicariance or dispersal events. Although rare, long-range dispersal events have occurred, as evidenced by the distribution of Caribbean clade species on landmasses that have never had land connections between them, and these events have played a role in shaping the evolutionary history of the ge-nus (Salas-Leiva et al. 2017). Eckenwalder (1980b posited that overwater dispersal of Caribbean zamias could occur via beached rafts of drift material carried by strong currents or during hurricanes. Additionally, several authors have reported birds as short-range dispersal agents of Zamia (Eckenwalder 1980b;Gómez 1993;Tang 1993;Chemnick 2007). Long-range dispersal could also be possible via larger birds that are able to swallow whole seeds and disperse them via defecation, such as the great Curassow (Crax rubra L.), a known dispersal agent of Ceratozamia whitelockiana Chemnick & T.J.Greg. (Chemnick 2007). The distribution range of Zamia, by far the broadest in the New World cycads and including several landmasses never connected by land, suggests that its dispersal abilities have played an important role in its evolutionary history.
Historical biogeography of Zamia. Inferring the historical biogeography of Zamia with any certainty is burdened by an extremely poor fossil record consisting of a single fossil that can be unequivocally assigned to Zamia based on morphological and epidermal characters (Erdei et al. 2018). The fossil, estimated to be of late Eocene to earliest Oligocene age (ca. 33-35 Ma), was found in marine sediments in central Panama.
Little is known about its terrestrial origin, but since parts of the Panama Arc were emergent since at least 30 Ma (O'Dea et al. 2016), it is possible that leaflets may have been carried by currents from nearby islands or more distant land. Interestingly, the cuticular micromorphology of the Panamanian fossil more closely resembles that of extant Caribbean clade species than that of Central-meridional clade species (Erdei et al. 2018). This may imply that the elongated adaxial cuticular cell shape and stomatal band width and ratio found in Caribbean species are ancestral character states. Two additional Paleogene (i.e., Oligocene 30 Ma) fossils from Collazo Shale (Puerto Rico) were originally described as separate species of Zamia (Hollick 1928) but are apparently conspecific (Erdei et al. 2018). Despite their Zamia-like macromorphology, their confirmation as Zamia is not possible due to poorly preserved cuticles (Erdei et al. 2018). If assignable to Zamia, the fossils would indicate the earlier presence of the genus in the Caribbean, albeit likely of an extinct lineage, as Puerto Rico was below sea level between the Late Oligocene and Early Pliocene (van Gestel et al. 1998). Hence, the current presence of Zamia on the island may be younger than ca. 5 Ma (Meerow et al. 2012), which is in agreement with our estimations suggesting that the diversification of extant Caribbean species occurred less than 2 Ma and within the last 340,000 yr for Puerto Rican and Dominican species (fig. 5).
The young age retrieved for the crown group of the Caribbean clade contrasts with its relatively old stem age, perhaps resulting from high turnover rates caused by glacio-eustatic sea level fluctuations, which can greatly influence island biotas (Weigelt et al. 2016). In fact, the effects of historic sea level fluctuations are readily observable in the genetic structure of extant species. Sea level fluctuations are thought to have reinforced isolation in Zamia erosa O.F.Cook & G.N.Collins limestone hill populations in northern Puerto Rico by eliminating intervening populations in the lowlands (Meerow et al. 2012). Similarly, populations of Zamia on five islands in the Bahamian archipelago cluster genetically into two groups (Salas-Leiva et al. 2017) that correspond to two separate historic landmasses (Little Bahamas Bank and Great Bahamas Bank), which encompassed these islands during the Last Glacial Maximum (33.0 and 26.5 Ka) 304 when sea levels were as much as 120-135 m below current levels (Clark and Mix 2002).
Based on the limited fossil record for Zamia, Erdei et al. (2018) suggested an Eocene or earliest Oligocene origin of the genus in the Central American-Caribbean region with subsequent dispersal south and northward. Our molecular biogeographic and divergence time estimation analyses instead suggest an older Late Cretaceous to early Eocene origin for the stem group of Zamia and a Caribbean or Mesoamerican origin for both the stem and crown nodes of Zamia with subsequent dispersal of the genus south into the Central American Isthmus region and South America. Despite the discovery in the Central American Isthmus region of the only anatomically verified Zamia fossil (albeit of unknown terrestrial provenance), a Caribbean or Mesoamerican origin for the genus appears more convincing due to the presence of the sister genus Microcycas as well as the two earliest diverging Zamia lineages (Caribbean and Fisheri clades) at the northernmost range of the genus's geographic range. Our biogeographic analysis suggests a northern origin for the crown group of Zamia with vicariance caused by a southward distribution shift in the Late Miocene, a subepoch characterized by increasing drought, enhanced seasonality, and a restructuring of terrestrial biological communities (Herbert et al. 2016). Discoveries of additional fossils of Zamia or its sister genus Microcycas may help to better recreate the still poorly understood early biogeographic history of the genus.
There is currently considerable disagreement regarding dating of the emergence of the Panama land bridge, with several authors advocating the traditional view of ca. 2.8 Ma (O'Dea et al. 2016) and others suggesting the isthmus has been in place for more than 10 Ma (Montes et al. 2012;Jaramillo et al. 2017). Either way, the estimated divergence of the Central-meridional clade into the Isthmus and South America clades (3.8 Ma, 95% HPD 2.49-5.17) suggests that the two groups diverged when the Panama land bridge was either nearly completed (O'Dea et al. 2016) or had been in place for millions of years.
Diversification of the crown nodes for both clades began nearly simultaneously about a million years later, with the South America clade's mean crown age estimated at 2.62 Ma (95% HPD 1.71-3.56) and the Isthmus clade's at 2.35 Ma (95% HPD 1.43-3.32). This period coincides with the establishment of repeated Pleistocene glaciations beginning at 2.6 Ma, which converted much of Central America into recurrently dry environments (Bacon et al. 2016) leading to the fragmentation of wet forest habitat favored by the Central-meridional clade and likely driving diversification of both the Isthmus and South America clades.
The South America clade is considerably more species rich and has a broader geographic range than the Isthmus clade, perhaps in part because the dispersal potential of the latter has been limited by the geographic bottleneck of the isthmus's narrow landmass. Despite this bottleneck, the Isthmus clade has been remarkably successful. It has the third highest net diversification rate after the Caribbean and South America clades (table 8), and the Isthmus region hosts the greatest concentration of species per unit area in the genus, with the main taxonomic diversity hot spots occurring in central Panama and along the Costa Rican and Panamanian border ( fig. 8).
The Isthmus and South America clades appear to have diversified and dispersed mostly within their geographic regions, despite the fact that this diversification occurred when the Panama land bridge was already present or near completed. The genetic isolation of these clades mirrors similar findings reported for birds (Smith and Klicka 2010) and mammals (Webb 1976) in which biotic interchange facilitated by uplift of the isthmus declined during the Pleistocene, perhaps driven by the unavailability of ecological space. Reduced migration by birds and mammals, both known dispersal agents of Zamia (Tang 1989), may have affected transisthmian dispersal of the genus, and the ecological niches available to Zamia may have also filled rapidly in both landmasses, therefore impeding biotic exchange. Finally, the coevolution of Zamia species with insect pollinators also may have played a role in keeping the clades separate. In fact, whereas erotylid beetle pollinators of the genus Pharaxonotha Reitter are found throughout most of the distribution range of Zamia, weevil pollinators of the subtribe Allocorynina Sharp are absent from South America clade species (O'Brien and Tang 2015).

Diversification
Most diversification above the crown node of Zamia occurred in the Pliocene and Pleistocene epochs, with only the split between the Caribbean and Mainland clade and between the Fischeri clade and the rest of the Mainland clade occurring in the Miocene (figs. 5, 6). The vast majority of extant Zamia species appear to be of Pleistocene origin, suggesting that climactic oscillation during this epoch may have played a role in diversification.
The net diversification rates (sensu Magallon and Sanderson 2001) recovered for the crown group of Zamia are among the highest reported for gymnosperm genera by Crisp and Cook (2011), yet comparable to rates obtained using clade age estimates from Nagalingum et al. (2011) and Condamine et al.  table 8). The diversification rates were highest in the Caribbean clade, followed by the South America and Isthmus clades, the two major clades constituting the Central-meridional clade (table 8). The high diversification rate recovered in the Caribbean clade may be partly due to the very short generation times found in this clade, as its constituent species are able to reach reproductive age as early as 2 yr from seed (Salas-Leiva et al. 2017), therefore allowing for a faster rate of accumulation of neutral mutations compared with most mainland species.
It has been argued that the Magallon and Sanderson (2001) diversification rate estimate represents a more accurate estimator to infer rate shifts than BAMM . Indeed, our BAMM analysis suggests there is no significant variation in rate between clades. However, the difference in sample size between the clades identified here might indicate that the extreme estimates obtained for the two clades with the highest and lowest rates (Caribbean clade, n p 9, and Fischeri clade, n p 3) are very poor estimators of the subtending variation (Rabosky 2018). Thus, we advise that the results obtained here using the Magallon and Sanderson (2001) methods should be interpreted with caution. A similar issue might affect our estimation of the gamma statistic for different subclades. Nagalingum et al. (2011) found patterns of early diversification followed by subsequent decreases in diversification rates for all cycad genera, a pattern typical of an adaptive radiation (Rabosky and Hurlbert 2015), which they attributed to newly available seasonal environments as the ecological opportunity that triggered a global radiation following the mid-Miocene transition. We found a signal for significant early rapid speciation only in the South America clade ( fig. 6), perhaps an indica-

306
INTERNATIONAL JOURNAL OF PLANT SCIENCES tion of adaptive radiation triggered by the ecological opportunity presented by newly colonized South America. All other clades had constant diversification rates following a "museum model" of evolution (Wallace 1878; Fischer 1960) with gradual species accumulation coupled with low extinction rates. This model is supported by the BAMM analysis, in which no significant diversification rate shifts were detected within Zamia or any of its constituent clades, and by the MCCR test, which indicated a constant diversification rate for the genus and most major clades ( fig. 6; table 9). In addition, no significant mass extinction signals were detected in the genus using the Comet model ( fig. D2). The stable diversification history apparent in Zamia differs from the findings of Yessoufou et al. (2014) in the genus Encephalartos Lehm., in which a mass extinction event has been postulated to explain a pattern of rate heterogeneity through time. This difference could indicate the ability of Zamia species to track their favorite habitat during periods of cooling-drying, or it could signify different patterns of cooling and aridification between Africa and the Neotropics (Couvreur 2015).

Conclusions
Our results suggest a crown age of ∼10 Ma for the genus and an origin within the Mesoamerican or Caribbean region with eventual southward migration into the Central American Isthmus and South America. The diversification history of the genus appears relatively stable and characterized by gradual species accumulation and low extinction rates. Most extant Zamia species appear to have originated during the Pleistocene, suggesting an effect of climatic oscillations during this epoch on the diversification of the genus. We recovered a strong congruence between geographic distribution and phylogenetic relationships, a pattern consistent with the limited dispersal and pollen-transfer abilities found in cycads. Homoplasy was pervasive within macromorphological characters normally considered informative for the genus, suggesting that future classification work within the genus must rely mostly on genetic characters or explore a wider variety of potentially diagnostic characters, including additional morphological, anatomical, and reproductive characters. This study presents the most comprehensive interpretation of the evolutionary history of Zamia to date and should serve as a strong phylogenetic framework for subsequent research, including studies that may help increase our understanding of the mechanistic bases for the diversification of the genus. Future phylogenetic work, particularly with new high-throughput sequencing technologies, should help resolve some of the ambiguities that still remain in obtaining a fully resolved, well-supported phylogeny of the genus.

Acknowledgments
Funding for field and laboratory research was provided by the National Science Foundation (DEB 1050340), National Geographic Society (8965-11), Mohamed Bin Zayed Species Conservation Fund (project no. 0925331), Association for Zoological Horticulture, Institute of Museum and Library Services grants (MA-05-12-0336-12: Mission Based Collections Planning and MA-30-14-0123-14: Mission Based Collections Stewardship), SOS-Save Our Species (grant 2012A-035), Tinker Foundation Field Research Grant, Montgomery Botanical Center, US Department of Agriculture Agricultural Research Service, Fairchild Tropical Botanic Garden, Charles P and Dorothy Sacher, Dr. Lin Lougheed, and Christiane Tyson Research Fellowship granted to D. Salas-Leiva. We are grateful to the many botanic gardens, herbaria, and horticulturists that generously provided samples for this study (see app. A). Kyoko Nakamura provided laboratory training and support to M. Calonje; Fabien Condamine provided helpful guidance regarding age divergence estimation; and Miguel Angel Pérez Farrera provided relevant information about Mexican zamias. We are grateful to the two anonymous reviewers and to the editor Hervé Sauquet for helpful comments that allowed us to greatly improve previous versions of the manuscript. Finally, Claudia Patricia Gutierrez provided unconditional family support to M. Calonje during his dissertation work, from which this article originated.

Appendix A Accessions Sampled for Molecular Analyses
The following information is provided for each sample: species name and authorship, sample ID, sample provenance or pedigree (country and state or province), representative herbarium specimen (collector, herbarium repository), DNA sample source (source and botanical garden accession numbers), Genbank accession numbers (from left to right, dashes indicate no sequence available: 40S, ATG2, CyAg, GroES, HTS, LiSH, PEX4, PMP22, psbK/I, WRKY4. Herbarium acronyms follow Index Herbariorum (http://sweetgum.nybg.org/science/ih/). Abbreviations used for DNA sample sources: FTBG: Fairchild Tropical