Evaluating the monophyly of Mammillaria series Supertextae (Cactaceae)

Abstract Mammillaria (Cactaceae) taxonomy has been historically problematic due to the morphological variability and sympatry of the species. This has led to several proposals for infrageneric classification, including subgeneric, section and series categories. Mammillaria ser. Supertextae is one of 15 series and is made up of a variable set of species that are mainly distributed in southern Mexico and Central America. However, the phylogenetic relationships within M. ser. Supertextae and its relationship to other Mammillaria taxa are far from fully understood. Here we attempt to elucidate these relationships using complete terminal sampling and newly obtained chloroplast marker sequences and comparing them to Mammillaria species sequences from GenBank. Our phylogenetic analyses showed that M. ser. Supertextae comprises a well-supported monophyletic group that diverged approximately 2.1 Mya and has M. ser. Polyacanthae as its sister group; however, relationships within M. ser. Supertextae remain unresolved. The topology obtained within M. ser. Supertextae must also be interpreted under the distribution shared by these taxa, but it is difficult to differentiate ancestral polymorphisms from possible introgression, given the short time elapsed and the markers used. Our results show that the infrageneric units of M. haageana and M. albilanata can be considered independent evolutionary units. We also suggest that the relationship between M. haageana and M. albilanata is convoluted because their distribution overlaps (mainly towards southern Mexico), with genetic differences that possibly indicate they represent more than two taxonomic entities. One possible explanation is that there could still be gene flow between these taxa, and we might be witnessing an ongoing speciation process.


Introduction
Mammillaria Haw. (Cactaceae, Cactoideae, Cacteae) is the most diverse genus within the cactus family, with a broad range of recognized species, ranging from 163 (Hunt et al. 2006) to 181 (Pilbeam 1999) up to 320 species (Reppenhagen 1992). Mammillaria is characterized by tuberculate stems, definite dimorphic areoles not connected by a groove, flowers that arise from the base of the tubercles and not apically, and seeds with testa cell walls that are par-concave and undulated (Bravo-Hollis and Sánchez-Mejorada 1991;Lüthy 1995;Anderson 2001;Hunt et al. 2006). This genus, together with Coryphantha (Engelm.) Lem., Escobaria Britton & Rose, Neolloydia Britton & Rose, and Ortegocactus Alexander, is integrated into the Mammilloid clade (Butterworth et al. 2002). However, it has been proposed that Mammillaria is a polyphyletic group. Breslin et al. (2021), using plastid genomes, confirmed previous studies showing Mammillaria is nonmonophyletic, as currently circumscribed, so they proposed that the Mammilloid clade be circumscribed in three monophyletic genera, Mammillaria s.s., Coryphantha and Cochemiea s.l., as previously suggested by Vázquez- Sánchez et al. (2013). Furthermore, the taxonomy of Mammillaria has historically been difficult due to large morphological variability, phenotypic plasticity, sympatric distribution of species, and suspected hybridization events.
Within Mammillaria, there are 15 recognized series (Hunt et al. 2006), one of which is M. ser. Supertextae D.R. Hunt, distributed from western Mexico to Central America and even the Caribbean islands (Pilbeam 1999). Mammillaria ser. Supertextae is a clear example of the species delimitation problem within the genus, illustrated by the number of accepted species ranging from 8 to 27, although the most recent taxonomic proposal recognizes only 9 taxa (Table 1). Morphometric and molecular studies have attempted to assess the proposed infrageneric relationships for Mammillaria, but none have been specifically directed at M. ser. Supertextae. Lüthy (1995) performed a detailed morphological analysis of the genus with a phenetic approach, in which he included four species of M. ser. Supertextae (M. albilanata Backeb., M. dixanthocentron Backeb. ex Mottram, M. haageana Pfeiff., and M. huitzilopochtli D.R. Hunt), showing that M. ser. Supertextae is characterized by the presence of extracellular crystals; however, the trait is not exclusive to M. ser. Supertextae, as the M. ser. Leucocephalae Lem. ex Schumann also shows this characteristic. Butterworth and Wallace (2004) conducted a molecular study using two chloroplast markers (rpl16 and psbA-trnH), including five species of M. ser. Supertextae (same as Lüthy but including M. supertexta Mart. ex Pfeiff.); although the species were grouped together, the support values were low (BS = 63, PP = 0.99), and a phylogenetic relationship could not be established with the remainder of Mammillaria. More recently, the complete sequencing of the chloroplast genome of eight Mammillaria species confirms that four M. ser. Supertextae taxa (M. crucigera Mart., M. supertexta, M. huitzilopochtli, M. haageana subsp. san-angelensis (Sánchez-Mej.) D.R. Hunt) represent a clade (Solórzano et al. 2019;Hinojosa-Alvarez et al. 2020).
To disentangle the evolution of Mammillaria, we decided to focus on elucidating the phylogenetic relationships of M. ser. Supertextae. We included all taxa proposed by Hunt et al. (2006), except for M. halbingeri Boed., as, according to Reppenhagen (1992), the species was not reported again. We also included 12 localities of M. haageana and seven of M. albilanata. All these species together constitute one taxonomic complex within M. ser. Supertextae (Arias et al. 2012). We chose two chloroplast markers, the rpl16 intron and the intergenic spacer psbA-trnH. In Cactaceae, both markers have been used to resolve phylogenetic relationships (Korotkova et al. 2010;Cruz et al. 2016;Barrios et al. 2020), including Mammillaria (Butterworth and Wallace 2004;Vázquez-Sánchez et al. 2013;Hernández-Hernández et al. 2014); therefore, there are many sequences available in GenBank that can be used to expand the sampling of terminals, including sister groups and outgroups essential for testing the monophyly (Korotkova et al. 2017) of M. ser. Supertextae. The main objective of this study was to test the monophyly of M. ser. Supertextae and estimate its divergence time by broadening the sample of terminals within the series.

Materials and methods
The present study included a total of 123 taxa, 111 species of Mammillaria, 5 closely related genera ( Sánchez et al. (2013). For the remaining 28 taxa of M. ser. Supertextae, we generated new sequence data. DNA was extracted from 40 mg of silica-dried (24 h) stems. The samples were stored at -80 °C, and 12 hours later, they were triturated in a TissueLyser II (Qiagen, Venlo, Netherlands) at 29 rpm for 25 s twice. Extraction was performed with the DNeasy Plant Mini Kit (Qiagen, Hilden, Germany) following the manufacturer's instructions, and the elution volume was 35 µl twice in Milli-Q water. The rpl16 intron and psbA-trnH intergenic spacer were amplified using standard PCR protocols. The rpl16 region was amplified using the primers from Hernández-Hernández et al. (2011), rpl161F (5ʹ-GCTATGCTTAGTGTGTGACTCGTT-3ʹ) and rpl163R (5ʹ-CTTCTATTTGTCTAGGCGTGATCC-3ʹ), by initially denaturing the DNA for 5 min at 94 °C, followed by 28 cycles of 1 min at 94 °C, 50 s at 55 °C, and 2 min at 72 °C, and a final extension of 4 min at 72 °C. The psbA-trnH intergenic spacer was amplified with primers from Korotkova et al. (2010), CApsbA (5ʹ-CCGTGCTAACCTTGGTATGG-3ʹ) and CAtrnH (5ʹ-CCGCGAATGGTGGATTCACAAT-3ʹ). PCR conditions were 2 min at 94 °C; followed by 29 cycles of 30 s at 94 °C, 30 s at 52 °C; and 1 min at 72 °C; and a final extension of 7 min at 72 °C. Amplifications were performed using 0.6 U of Platinum Taq polymerase according to the manufacturer's protocol (Invitrogen, Carsbad, California, USA), 4 mM mixed dNTPs (Invitrogen, Thermo Fisher Scientific, Waltham, USA), 1.5 mM MgCl2, 16 mg/mL BSA, 0.25 µM each primer, and 30-50 ng of genomic DNA in a reaction volume of 25 µL. The PCR products were sequenced in an Applied Biosystems Sequencer Model 3730xL at the Laboratorio de Biología Molecular de la Biodiversidad y de la Salud, Instituto de Biología, UNAM.
The sequences were aligned in GUIDANCE2 (v. 2.02, Sela et al. 2015) using MAFFT (v. 7.407, Katoh and Standley 2013). The algorithm was implemented with 100 iterations (-msaProgram MAFFT -MSA_Param "\-globalpair \--maxiterate 100" -bootstraps 100). The resulting matrix was imported into PHYDE (v. 0.9971, Müller et al. 2019) to manually edit ambiguously aligned sites. Both genes were concatenated, and data partitions were determined with the program PARTITIONFINDER (v. 2, Lanfear et al. 2017) using the Bayesian information criterion and greedy search. The model with the best fit for both markers was TVM + I + G. Indels and inversion were manually coded with MESQUITE (v 3.6, Maddison and Maddison 2019) using the simple coding method of Simmons and Ochoterena (2000).
Bayesian inference (BI) analysis was performed using MRBAYES (v. 3.2.1, Ronquist et al. 2012). The analysis was run with four Markov chain Monte Carlos (MCMC) (nchains = 4) and ten million generations (ngen = 10000000), sampling trees every 100 generations (samplefreq = 100) and discarding the first 25% as burnin. All parameters were monitored with TRACER (v. 1.7.1, Rambaut 2018) until they had effective sample sizes (ESS) of greater than 200. Maximum likelihood (ML) analysis was conducted using RAXML (v. 8.2.12, Stamatakis 2014) with molecular and binary partitioning, calculating the conditional likelihood of no invariant data and considering CAT for the heterogeneity of rate. The correction was made with the Lewis (2001)  To estimate divergence times, we used the credibility interval around the estimated age of the Mammilloid clade (5.83-12.56 Mya; Hernández-Hernández et al. 2014). We inferred a time-calibrated phylogenetic tree using a BI approach implemented in BEAST (v. 2.6.1, Bouckaert et al. 2019). Analysis of the concatenated matrix used the uncorrelated lognormal relaxed clock (Drummond et al. 2006) for a total of 20 million generations of MCMC, sampling once every 10000 trees and discarding 15% as burn-in using TREEANNOTATOR v. 2.6.0.

Results
The overall sequence matrix for the two genes included 2257 bp and 8 encoded indels. We excluded 1045 bp in rpl16 due to uncertain homology. The final length of the aligned matrix for rpl16 was 897 bp and 315 bp for psbA-trnH, with 168 and 69 potentially informative sites, respectively. The BI and ML analyses produced trees with similar topologies (Fig. 1). Mammillaria ser. Supertextae was recovered as a monophyletic group (PP = 1, BS = 84), supported by a transversion in psbA-trnH. Mammillaria ser. Polyacanthae was also recovered as a monophyletic group (PP = 1, BS = 91), supported by 2 transitions in psbA-trnH. The sister relationship between M. ser. Supertextae and M. ser. Polyacanthae (PP = 0.99, BS = 80) is supported by a deletion in rpl16.
A polytomy formed within M. ser. Supertextae, where four clades were formed (S1, S2, S3, and S4), three of which are defined by specific polymorphisms (i.e., clade S1 a transversion in rpl16, clade S3 an inversion in rpl16). In three clades, at least one terminal of M. albilanata subsp. oaxacana was confirmed (S1: CC044 and CC046; S2: CC040; S3: CC036); in addition, part of its geographic distribution was common to M. haageana ( Fig. 2A, B). The 12 terminals of M. haageana are distributed into two clades. S1 is formed by six terminals of M. haageana (CC024, CDMX, CC023, CC025, CC027, and CC045), two terminals of M. albilanata subsp. oaxacana referred  oaxacana (CC036) shows an inversion in rpl16 of 46 bp, and its sister group was M. lanata. Clade S4 is made up of two species, one of which corresponds to M. columbiana, which is distributed from Yucatan, Mexico to Colombia and Venezuela ( Fig. 2A); the second species is M. eriacantha, which is distributed in Veracruz, Mexico ( Fig. 2A, B).

Discussion
The concatenation of two matrices (rpl16 and psbA-trnH) and extensive sampling (eight of nine species, according to Hunt et al. 2006; Table 1) helped to recover M. ser. Supertextae as a monophyletic group, consistent with previous molecular phylogenetic studies that included only five (Butterworth and Wallace 2004)   uncertain, and it was placed within M. ser. Polyacanthae due to the size of its flower (Bravo-Hollis and Sánchez-Mejorada 1991;Hunt et al. 2006). Remarkably, we show that M. eriacantha is nested within M. ser. Supertextae, as previously proposed by Reppenhagen (1992) and Lüthy (1995) based on the presence of extracellular crystals. Our results also show that M. ser. Polyacanthae is the sister series of M. ser. Supertextae and that they are part of the Mammillaria sect. Subhydrochylus (Hunt 1987). Within M. ser. Supertextae, phylogenetic relationships were not resolved; however, both M. haageana and M. albilanata appeared in more than one clade.
Intron rpl16 was demonstrated to be the most variable and informative marker compared to the intergenic spacer psbA-trnH, which is consistent with previous studies in Cactaceae (Korotkova et al. 2010;Vázquez-Sánchez et al. 2013;Cruz et al. 2016). The molecular characteristics that define the monophyly of M. ser. Supertextae and M. ser. Polyacanthae were found in psbA-trnH, while a deletion in rpl16 supports the sister relationship. Although the deletion in rpl16 is partial in M. ser. Supertextae and M. ser. Polyacanthae, a total deletion has been reported in members of M. ser. Stylothelae (Butterworth et al. 2007), showing that deletions in rpl16 in Mammillaria may be a strong characteristic for the identification of infrageneric groups. It remains to be defined whether other polymorphisms, such as the inversion of another chloroplast gene, trnF-GAA, are also diagnostic of these clades within and between series, as recently reported for M. crucigera, M. huitzilopochtli, and M. supertexta (Solórzano et al. 2019).
The Mammilloid clade originated approximately 8.62 Mya (95% HPD = 5. 83-12.56;Hernández-Hernández et al. 2014), and within Mammillaria, it is likely that M. ser. Supertextae is a recently divergent group that originated in the Neogene-Quaternary transition approximately 2.1 Mya (95% HPD = 0.91-3.47). In some geographic regions, the M. ser. Supertextae species have undergone tectonic, erosive, alluvial and volcanic changes for millions of years; during the Pleistocene, these processes continued, giving rise to the current geomorphology (Siebert and Carrasco-Núñez 2002;Medina-Sánchez et al. 2020). Paleontological and molecular evidence suggests that glacial climate cycles that occurred during the last 2.5 Mya affected the distribution, diversity, and genetic structure of plant and animal populations (Gámez et al. 2014;Scheinvar et al. 2016;Cornejo-Romero et al. 2017). The hypotheses suggest that during the Pleistocene, these species sought refuge during adverse environmental conditions and expanded again when conditions improved (Scheinvar et al. 2016). Our results show that within M. ser. Supertextae, four clades are formed, two of which have distinctive climatic and topographic characteristics: Clade S1 (M. haageana, M. albialanta, M. dixanthocentron and M. supertexta), with species that are distributed in warm zones mainly at altitudes ranging from 447 to 2318 meters in thorn and tropical deciduous forests; and Clade S2 (M. haageana, M. albilanata and M. flavicentra), with species that are distributed in temperate zones at altitudes that range mainly from 1285 to 2518 meters in pine-oak forests. The environmental, geological and topographic differences between closely related species produced during climatic changes suggest differential selection pressures and local adaptation, which could have driven the speciation process (Mastretta-Yanes et al. 2015;Aquino et al. 2021), as has been suggested for Mammillaria pectinifera (Cornejo-Romero et al. 2014), Cephalocereus columna-trajani (Cornejo-Romero et al. 2017) and the genus Epithelantha (Aquino et al. 2021). Mammillaria haageana and M. albilanata represent a complex that extends widely in southern Mexico. Our results show that the infrageneric units of M. haageana and M. albilanata can be considered independent evolutionary units. It is possible that the variation in these inhabited environments promotes divergence in these taxa, although more in-depth studies are needed to understand and corroborate the hypotheses raised here.
The chloroplast marker sequences that we used (rpl16 and psbA-trnH) were not sufficient to establish the relationships among the taxa within M. ser. Supertextae. This was not surprising, as chloroplast markers have been used to resolve relationships at the species level; however, they have limitations when the species are closely related (Yan et al. 2018). This is because recently diversified groups may generate complicated genetic patterns, such as incomplete lineage sorting and hybridizations and/or introgressions (Li et al. 2016;Goetze et al. 2017), which may be true for M. ser. Supertextae. In other taxa (e.g., Petalidium Nees, Acanthaceae; Tripp et al. 2017), these problems have been addressed using multiple-locus methods to infer genetic trees, although they require nuclear markers that are not linked with levels of sequence variation according to phy-logenetic questions (Eaton and Ree 2013). To date, no effective nuclear markers have been developed for Cactaceae, and existing markers provide fewer informative sites than chloroplast markers (Cruz et al. 2016). Recently, proposals for nuclear markers have been generated through mining strategies to test hybridization in Opuntia species (Granados-Aguilar et al. 2020). Currently, several methodologies have been designed that allow biological questions to be answered using a reduced representation of the genome (Anderson et al. 2017;Choquet et al. 2019;David et al. 2019). This confers advantages when working with nonmodel species such as M. ser. Supertextae, since genomic markers can be genotyped in many individuals at low cost, and in most cases, it is not necessary to have a priori information such as a reference genome (da Fonseca et al. 2016).
The taxonomic proposals of M. ser. Supertextae species have been mainly based on interpretations according to the author's experience (Table 1), and their relationships have not been specifically tested under phylogenetic methods. Our methods are intended to be systematic (explicit criteria) and reproducible. Under this scheme, the results obtained show that within M. haageana and M. albilanata, there are genetic differences possibly indicating that these species comprise more than one taxonomic entity. Nevertheless, when distinguishing between M. haageana and M. albilanata, the task becomes difficult because both share similar distributions and habitats (mainly in southern Mexico; Fig. 2), and the morphological differences have not been well defined (Arias et al. 2012). A possible hypothesis is that there could still be gene flow between these taxa, and we might be witnessing an ongoing speciation process.

Conclusion
By including most of the species recognized by Hunt et al. (2006), our results show that M. ser. Supertextae is monophyletic, and we corroborate that M. eriacantha is part of the series as previously proposed. We find that M. ser. Polyacanthae is the sister series, as proposed by (Hunt 2011). The results also showed that M. ser. Supertextae is a recently diverged group. This is a first approximation to understand the evolutionary processes within M. ser. Supertextae. Future work should test sequencing techniques that allow genomic markers to be genotyped in many individuals since it is possible that conflicts in the phylogeny were the result of reticulate evolution. Furthermore, disentangling this problem will require a comprehensive pool of approaches regarding morphology and ecology, opening an avenue to develop M. ser. Supertextae as a model for studying complex evolutionary processes in Mammillaria.