Materials and methods

Plant material and taxon sampling. Plant material of Aporocactus species was collected from wild locations across the states of Hidalgo, Querétaro, Oaxaca, Puebla, and Veracruz in the springs of 2015 and 2016. Sampling included the type localities for the published names (when included in the protologue). For each locality, a section of stem was collected, and a fragment was subsequently herborized and deposited in MEXU; the second fragment was cultivated in the tempered greenhouse in the Botanical Garden of the Institute of Biology at UNAM (JB-IBUNAM), where a tissue sample was obtained, dried and stored in silica gel at -20 °C for subsequent DNA extractions. We included 50 taxa from Hylocereeae as ingroups, 21 of which corresponded to different localities of Aporocactus (Appendix 1), and the remaining 35 taxa corresponded to the genera Acanthocereus (Engelm. ex A. Berger) Britton & Rose, Disocactus Lindl., Epiphyllum Haw., Pseudorhipsalis Britton & Rose, Selenicereus (A. Berger) Britton & Rose, and Weberocereus Britton & Roses from the same tribe. The outgroup consisted of seven species from seven genera pertaining to the sister tribes: Bergerocactus Britton & Rose, Cephalocereus Pfeiff., Stenocereus (A. Berger) Riccob., Echinocereus Engelm., Deamia Britton & Rose, Myrtillocactus Console, Marshallocereus Backeb., and Leptocereus quadricostatus Britton & Rose. Sampled taxa in each analysis are described below.

Isolation, amplification and sequencing of DNA. For the isolation of total genomic DNA, most of the water-storing tissue was removed from the stems before the remaining cortex tissue was dehydrated in silica gel. The dried plant material was homogenized using a mixer mill (Retsch MM200, Haan, Germany) and extracted using the EZ-10 mini-prep kit for plant genomic DNA (Bio Basic, Inc., Ontario, Canada) following the manufacturer's protocol. The incubation time in the lysis buffer was increased to 120 min at 65 °C due to the tissue type. The concentration and purity of DNA (A260/A260 and A260/A230 ratios) were measured using a spectrophotometer (NanoDrop, peqLab, Erlangen, Germany). The original genomic DNA was stored at -20 °C and working dilutions with a standard concentration of 10 ng/μl were prepared for subsequent analysis in PCR assays. PCR amplification was performed for the rpl16 intron (^{Hernández-Hernández et al. 2011}), trnL-trnF intron (^{Taberlet et al. 1991}), psbA-trnH intergenic spacer (^{Sang et al. 1997}, ^{Tate & Simpson 2003}) and trnQ-rps16 intergenic spacer (^{Korotkova et al. 2010}, ^{Shaw et al. 2007}). The total volume for the standard sample was 25 µl, which consisted of 2.5 µl of 10X buffer, 0.5 µl dNTPs at 200 µM concentration, 1 µl of BSA, 0.75 µl of MgCl2, 0.3 µl F primer, 0.3 µl R primer, 1.25 µl of DNA Platinum Taq Polymerase (Invitrogen™) at 5 U/µl, 0.6 µl of total genomic DNA and 19.025 µl of H₂O. The markers that employed internal primers for sequencing were adjusted to a total volume of 50 µl. The PCR programmes used for each marker were as follows: 1) trnQ-rps16, denaturation at 95 °C × 2’, denaturation at 95 °C × 1’, annealing at 55 °C × 1’, extension at 72 °C × 1’, and extension at 72 °C × 7', for 35 cycles. 2) rpl16/trnL-trnF, denaturation at 95 °C × 2’, 94°C × 1’, annealing at 54 °C × 1’, extension at 72 °C × 1’ 30’’, and extension at 72 °C × 7’, for 30 cycles. 3) psbA-trnH, denaturation at 95 °C × 2, denaturation at 95 °C × 30’’, annealing at 55 °C × 1’, extension at 72 °C × 1’, and extension at 72 °C × 10’, for 30 cycles. The sequencing of the molecular markers was performed in the Laboratory of Genomic Sequencing of Biodiversity and Health from the Biology Institute at the National Autonomous University of Mexico (UNAM).

Sequence alignment. The sequences from Aporocactus samples were quality-checked, assembled and edited using Sequencher^® v. 4.8 (Gene Codes, Ann Arbor Michigan USA). The sequences for the species of the genera Acanthocereus, Disocactus, Epiphyllum, Pseudorhipsalis, Strophocactus, Bergerocactus, Cephalocereus, Deamia, and Marshallocereus were obtained from the database of the Laboratory of Systematics of Cactaceae from the Botanical Garden/Institute of Biology, UNAM (^{Arias et al. 2005}, ^{Cruz et al. 2016}, ^{Sánchez et al. 2014}, ^{Hernández-Hernández et al. 2011}, ^{Tapia et al. 2017}) (Appendix 1). Additionally, we included the rps3-rpl16 and trnK-matK sequences from ^{Korotkova et al. (2017)} to complete the matrix (Appendix 1). Individual sequences were cross-checked for possible assembly failures and subsequently stacked and subjected to primary alignment using the software BioEdit (^{Hall 1999}) and the integrated application ClustalW v.1.74 (^{Thompson et al. 1994}). Furthermore, individual marker matrices were realigned and corrected by eye using Mesquite^® software v. 3.03 (^{Maddison & Maddison 2016}).

Phylogenetic analyses. A phylogenetic analysis for delimiting species was performed by using four cpDNA markers (psbA-trnH, trnQ-rps16, rpl16, and trnL-F), including 21 samples of Aporocactus and 16 species from eleven genera of Hylocereeae. On the other hand, a phylogenetic analysis for recovered genus relationships used six cpDNA markers and included 35 species from 15 genera. For both analyses, the cpDNA matrix consisted of six markers: psbA-trnH, trnQ-rps16, rpl16, trnL-F, trnk-matk, and rps3-rpl16. The parameters of the Bayesian analyses were identical for both analyses and were performed in MrBayes v. 3.2.1 (^{Huelsenbeck & Roquist 2001}, ^{Ronquist et al. 2012}). The General Time Reversible model (GTR+I+G) was selected as the best substitution model using the Bayesian Information Criterion (BIC), as implemented in jModeltest v. 2.0 (^{Darriba et al. 2014}). The analyses consisted of 10 million generations, sampling of parameters and trees every 1,000 generations, and a burning of 25 % of the resulting trees. The convergence of the chains was evaluated visually from the resulting parameter archive of MrBayes using Tracer v. 1.6 (^{Rambaut et al. 2018}).

Ecological niche modeling. We constructed ecological niche models (ENMs) to predict the current distribution of suitable habitat of the recognized species of Aporocactus. Geographic coordinates of occurrence of each species were obtained from field collection, MEXU herbarium specimens, and unambiguous records from Naturalista (www.naturalista.mx). We discarded duplicate records, records with doubtful identity or geographic location and records from cultivated plants. The accessible area (M area, ^{Soberón & Peterson 2005}) was defined by the genus range based on the biogeographical provinces proposed by ^{Morrone et al. (2017)} and the distribution of pine-oak vegetation and cloud forest associate to those provinces (^{Rzedowski 1990}). Bioclimatic variables were used at an ~1 km² spatial resolution compiled by ^{Cuervo-Robayo et al. (2014)}. We masked those climate layers to the extent of the M area. To avoid collinearity, we discarded one of the bioclimatic variables that was highly correlated with another (Spearman correlation values > 0.79) for the study area. Nine variables were used in the final analysis (BIO2, BIO4, BIO10, BIO11, BIO13, BIO14, BIO15, BIO18, and BIO19). For each species, we constructed an ENM using MAXENT v. 3.4.1 (^{Phillips et al. 2017}) through package “dismo” in R v. 4.0.4 (^{R Core Team 2020}). We thinned occurrence points to 1 km² to avoid spatial autocorrelation. We built different models with 10,000 random background points and evaluated them with spatial-cross validation. We used no campling and different parametrization for Maxent, combining regularization multipliers in intervals of 0.5 ranging from 0.5 to 5, and feature class combinations of Linear, Quadratic, Hinge and Product: L, H, LQ, LH, LQH, and LQHP. We performed the evaluation process with the spatial cross validation procedure “random k-fold” (number of folds = 4) using the R package ENMeval v. 2.0.1 (^{Kass et al. 2021}) with R. Model selection was made based on the Akaike information criteria corrected for small sample sizes (ΔAICc), that reflects a comparison of the goodness-of-fit and parsimonious model (^{Muscarella et al. 2014}). We projected the models using the Maxent “cloglog” transformation. Finally, we evaluated variable importance with Maxent´s variable jackknife test (^{Phillips et al. 2006}). Final models were constructed with ten cross-validation replicates without extrapolation.

Niche identity and similarity. The differences between the niches of the species recognized in Aporocactus were evaluated by using niche overlap, niche identity, and niche similarity analyses in ENMtools (^{Warren et al. 2010}). Niche overlap was calculated through Schoener’s index (D) and Hellinger’s-based I index, which measures the similarity between predictions of habitat suitability (ENM) of one or more pairs of species (^{Warren et al. 2008}, ²⁰¹⁰). The niche identity test indicates whether the ENMs produced by two species are identical. The test pools the georeferenced data points for a pair of species, randomizes the taxon identities of these data points, and extracts two new samples with the same sizes as the two original samples. This process is replicated and generates a null distribution of overlap scores, which is compared with the empirical niche overlap scores (^{Warren et al. 2010}). The background similarity test compares the ENM of taxon “A” to an ENM created from n random points drawn from the geographic range of taxon “B”, which generates a null distribution of overlap scores (^{Warren et al. 2008}, ²⁰¹⁰). This method is subsequently repeated in the other direction for both taxa in the comparison (B vs. A background). Finally, the test compares the empirical niche overlap of two taxa to a null distribution of overlap scores generated. A total of 100 replicates were run for the niche identity test and background similarity test to assess the differences between the habitat suitability scores defined in the ENMs for both species.

Results

Species delimitation analysis. Four molecular markers were amplified for the ingroup and the outgroup species (Appendix 1). The matrix for the species delimitation analysis was 3,354 bp in length from four concatenated molecular markers (psbA-trnH, rpl16, trnL-F, and trnQ-rps16). Phylogenetic analysis for species delimitation recovered the genus Aporocactus as a monophyletic group (posterior probability (pp) = 1, Figure 1). Two main clades were observed for Aporocactus. One clade included 13 samples from the states of Queretaro, Hidalgo, and Veracruz, which represented the putative taxa A. flagelliformis, A. flagriformis, and A. leptophis. None of those taxa was recovered as a monophyletic group. This clade was well supported (pp = 1) by 11 substitutions: three in psbA-trnH, two in rpl16, and six in trnL-F. The second clade was composed of eight terminals from Oaxaca and Veracruz and included the putative taxa A. martianus and A. conzattii. This second clade was well supported (pp = 1) by four molecular sites: one in rpl16 and three in trnL-F (positions 1,769, 2,000, 2,456). Additionally, the three samples of A. conzattii were recovered in a monophyletic group (pp = 1).

Figure 1 Species delimitation in Aporocactus. Cladogram of the majority rule consensus tree from the Bayesian analysis of the concatenated trnQ-rps16, trnL-trnF, psbA-trnH, and rpl16 markers. Numbers above branches are the Bayesian posterior probability values. 

Phylogenetic relationships analysis. The alignment to infer the phylogenetic relationships of Aporocactus was 6,920 bp in length from six concatenated DNA markers (psbA-trnH, rpl16, trnL-F, trnQ-rps16, trnk-matk, and rps3-rpl16). The analysis to infer the phylogenetic relationships of Aporocactus recovered three principal clades with good support: the hylocereoid clade (pp = 0.9), the phyllocactoid clade (pp = 1), and the Acanthocereus clade (pp = 1) (Figure 2). Aporocactus was resolved as a well-supported monophyletic group (pp = 1) in the hylocereoid clade and was positioned in an early divergent group sister to Selenicereus and Weberocereus (pp = 0.9). In this analysis, the genera Disocactus, Epiphyllum, and Pseudorhipsalis were nested in the phyllocactoid clade, while Acanthocereus was recovered as the earliest diversified lineage in Hylocereeae (Figure 2). The relationship between hylocereoid and phyllocactoid clades in this analysis had low support (pp = 0.7).

Figure 2 Phylogenetic relationships in Aporocactus. Cladogram of the majority rule consensus tree from the Bayesian analysis of the concatenated trnk-matk, rps3-rpl16, trnQ-rps16, trnL-trnF, psbA-trnH, and rpl16 markers. Numbers above branches are the Bayesian posterior probability values. 

Distribution, ecological niche modeling, and niche comparison. Based on Figure 1, the A. flagelliformis clade was determined to be primarily distributed in the Sierra Madre Oriental (^{Morrone et al. 2017}) through Querétaro, Guanajuato, Hidalgo, northern Puebla, and central Veracruz; while the A. martianus clade occupies primarily Sierra Madre del Sur (^{Morrone et al. 2017}) from central Veracruz to southern Puebla and Oaxaca (Figure 3A). The distribution limits of both clades of Aporocactus were observed to converge in central Veracruz state, where Sierra Madre Oriental and Sierra Madre del Sur intersect with the Mexican Transvolcanic Belt (^{Morrone et al. 2017}). Both species were determined to be clearly distributed in pine-oak forests and cloud forests in those biogeographical regions. Accordingly, these clades were recognized as different species: A. flagelliformis and A. martianus. Those clades were determined to be congruent with the current taxonomy of the genus (see discussion).

Figure 3 Actual and potential distribution of Aporocactus. A) Actual distribution of the genus Aporocactus, PdO: Pico de Orizaba, CP: Cofre de Perote, Xa: Xalapa volcanic field, CHPs: Chiconquiaco-Palma Sola. B) ENM of Aporocactus flagelliformis. C) ENM of Aporocactus martianus. 

Selected ecological niche model (ENM) for A. flagelliformis presented LQ features and regularization multiplier of 0.5 (ΔAIC ≈ 0, Table S1 and Figure S1). The ENM showed the AUC value = 0.947 (S2). Projected ENM of A. flagelliformis added as suitable areas a number of pine-oak and cloud forests in Nuevo León, Tamaulipas, southern Veracruz, and Oaxaca (Figure 3B). The variable with the highest percent contribution in the A. flagelliformis ENM was BIO18 (precipitation of warmest quarter) (24.3 %), followed by BIO14 (precipitation of driest month) (19.3 %), and BIO4 (temperature seasonality) (16.5 %). Variables with the highest permutation importance were BIO4 (36.8 %) and BIO18 (18.6 %). In the case of A. martianus, selected ENM presented LQH features and regularization multiplier of 2 (ΔAIC ≈ 0, Table S2 and Figure S2). This ENM showed an AUC value = 0.928 (S4). Projected ENM of A. martianus added some areas of pine oak forest in northern Puebla and Veracruz and northern Guerrero as suitable areas for the species (Figure 3C). The variables with the highest contribution to the ENM of A. martianus were BIO2 (mean diurnal range) (43 %) and BIO18 (28.8 %). The variable with the highest permutation importance was BIO2 (62.5 %).

Niche analyses indicated that empirical niche overlap between A. flagelliformis and A. martianus was low for de D index (D = 0.261); and moderate for the I index (I = 0.654). The identity test indicated that the ENM between the two species was significantly different (D_H0 = 0.772 ± 0.038 vs. D_H1 = 0.261 and I_H0 = 0.947 ± 0.017 vs. I_H1 = 0.654) (Figure 4A). The background similarity test comparing A. flagelliformis ENM in the A. martianus background and vice versa showed that the observed values of empirical niche similarity (D = 0.261, I = 0.654) were lower than expected under the null distribution (Figure 4B, C), indicating that the niches of the two species were significantly different than expected by chance in the available background environments.

Figure 4 Niche conservatism inference in Aporocactus. A) Niche identity test, green bars: D index frequency from null distribution, pink bars: I index frequency from null distribution, green arrow: empirical niche overlap D index, pink arrow: empirical niche overlap I index. B) Niche similarity test of Aporocactus flagelliformis as focus species and A. martianus as background. C) Niche similarity test of Aporocactus martianus as focus species and A. flagelliformis as background. For B) and C), blue bars: D index frequency from null distribution, orange bars: I index frequency from null distribution, blue arrow: empirical niche overlap D index, orange arrow: empirical niche overlap I index. 

Discussion

Species delimitation in Aporocactus. Considering monophyly as a property to recognize species, as well as the geographic distribution and floral morphology of each clade, our results indicated that the two clades in Aporocactus represent two different species (separately evolving metapopulation lineages, ^{De Queiroz 2007}). The first clade is formed by the samples initially identified as A. flagelliformis, A. flagriformis, and A. leptophis, but no internal group is formed based on these putative names or by their geographic origin; therefore, in this study, we recognize that samples comprise one species. It is worth mentioning that samples corresponding to the name A. leptophis and A. flagriformis were collected in their respective type localities (Zimapán, probably los Mármoles, Hidalgo and San José del Oro, Hidalgo, respectively). However, those have morphological features corresponding to the variation reported for A. flagelliformis. All samples included in this clade from Querétaro, Hidalgo, and northern Veracruz present zygomorphic flowers and magenta tepals (Figure 5A, B, C, D). Aporocactus flagelliformis (L.) Lem. (≡Cactus flagelliformis L.) is the first published name of the three samples mentioned above, and according to the principle of priority (Art. 11, ^{Turland et al. 2018}), it is the correct name for this species. The second clade includes the samples previously identified as A. conzattii and A. martianus (Figure 1). All specimens were distributed from central Veracruz to Oaxaca and exhibited actinomorphic symmetry with red tepals (Figure 5E, F, G, H, I, J). In this case, the name Aporocactus martianus (Zucc.) Britton & Rose (≡Cereus martianus Zucc.) has priority. This result is in keeping with the proposal of ^{Hunt (1989)}, who discussed the recognition of a northern species with zygomorphic purplish pink flowers and a southern species with regular scarlet flowers, assigning names on base to the ancient name. ^{Hunt (1989)} considered A. flagriformis and A. leptophis as stem and flower variations of A. flagelliformis and considered that A. conzattii is a re-description of A. martianus. Notably, a subclade was recovered with the samples of A. conzattii (Figure 1) from the Sierra Madre de Oaxaca at the Sierra Madre del Sur province. However, no particular character was observed in those samples of A. conzattii (Figure 5G), and this group probably represents the population genetic structure of A. martianus. We did not observe any infraspecific entity in A. martianus. Our results agree with the current taxonomy of Aporocactus, which recognizes two species for the genus (see Taxonomic treatment section in ^{Korotkova et al. 2017}). Wide variation in flower colour and size was observed, ranging from pink to magenta and from 4 to 7 cm in A. flagelliformis and from light red to deep red and from 7 to 12 cm in A. martianus (Figure 5).

Figure 5 Aporocactus flowers and their variation in color and sizes. A-D) Aporocactus flagelliformis, pink to magenta flowers, all zygomorphic [A, S. Arias 1225, Hidalgo; B, I. Rosas 006, Querétaro; C, I. Rosas 022, Veracruz; D, I. Rosas 024, Hidalgo]. E-J) Aporocactus martianus, light red to deep red flowers, with short to long receptacular tube, actinomorphic [E, M. A. Cruz 09, Oaxaca; F, I. Rosas 17, Oaxaca; G, I. Rosas 14, Oaxaca; H, I. Rosas 15, Oaxaca; I, M. A. Cruz 02, Veracruz; J, I. Rosas 08, Oaxaca]. 

Phylogenetic relationships of Aporocactus. The results supported the monophyly of the genus Aporocactus (^{Cruz et al. 2016}, ^{Korotkova et al. 2017}) and rejected the hypothesis of some authors that Aporocactus is a member of Disocactus because of the similarity in the shape, colour, and diurnal anthesis of these plants, which are presumably pollinated by hummingbirds (Barthlott in ^{Taylor & Hunt 1991}, ^{Bauer 2003}, ^{Hunt et al. 2006}). These results indicated that Aporocactus and Disocactus are independent lineages in different clades and suggest that diurnal anthesis in bright-coloured flowers appeared independently at least two times in Hylocereeae. In the sister tribe Echinocereeae, hummingbird pollination syndrome independently evolved in Morangaya pensilis (K. Brandegee) G.D. Rowley, Echinocereus section Triglochidiati Bravo, Stenocereus alamosensis (J.M. Coult.) A.C. Gibson & K.E. Horak and S. kerberi (K. Schum.) A.C. Gibson & K.E. Horak (^{Sánchez et al. 2014}). ^{Martínez-Quezada et al. (2020)} postulated that Aporocactus has two anatomical synapomorphies in the stem: 1) a delay in fibre development in the wood and 2) cortical bundles with secondary growth. In field work, we observed that Aporocactus plants do not develop wood, as occurs in other genera, such as Disocactus or Selenicereus; instead, in the base of the oldest stem in Aporocactus, the roots release them and promote vegetative propagation.

Aporocactus was recovered as a sister to Selenicereus + Weberocereus in the hylocereoid clade. This result was significant, since ^{Korotkova et al. (2017)} did not recover these relationships by using cpDNA markers only. We noted that the addition of cpDNA markers in the present study results in a more resolved phylogeny. This sisterhood (Aporocactus (Selenicereus and Weberocereus)) was also achieved by ^{Martínez-Quezada et al. (2020)} by using the cpDNA markers from ^{Korotkova et al. (2017)} and a complement of morpho-anatomical characters. ^{Martínez-Quezada et al. (2020)} suggest that the hemiepiphytic condition and the presence of adventitious roots along the stem represent the synapomorphies of this clade. Nevertheless, other members of Hylocereeae, such as Disocactus and Epiphyllum (phyllocactoid clade), can develop this type of root frequently in different stages of growth (juvenile, adult); rather, this root represents a homoplasy, which in combination with other characters is useful to diagnose the hylocereoid clade. It is important to highlight that in the absence of more DNA sequences, the addition of morphological characters can be useful for obtaining a more resolved topology, as observed in other cacti (^{Sánchez et al. 2018}, ^{Vargas-Luna et al. 2018}, ^{Martínez-Quezada et al. 2020}).

Distribution of Aporocactus. The known distribution of Aporocactus (Figure 3A) was restricted to the old pine-oak and cloud forests. As suggested by ^{Hunt (1989)}, A. flagelliformis represents the northern species through the Sierra Madre Oriental and extends to central Veracruz in the Transmexican Volcanic Belt. Traditionally, the distribution of A. martianus was only reported in Oaxaca at the Sierra Madre del Sur; however, our results showed that this species is also distributed in central Veracruz, at the limit of the Transmexican Volcanic Belt. Although the distribution of both species converges in central Veracruz, a detailed analysis of this region indicated that A. flagelliformis and A. martianus present an allopatric distribution. Our results suggested that speciation of the ancestral Aporocactus lineage was influenced by the formation of the modern Transmexican Volcanic Belt in the eastern part during the late Pliocene-Quaternary (2.0-0.1 ma) (^{Rodríguez et al. 2010}). A similar biogeographic pattern is also observed in other epiphytic sister species, namely, Disocactus phyllanthoides and D. ackermannii (^{Cruz et al. 2016}). Even the vicariant consequence of the Transmexican Volcanic Belt can be observed in sister species, such as Cephalocereus senilis and C. columna-trajani (^{Tapia et al. 2017}), in the lower western parts of the Sierra Madre Oriental and Sierra Madre del Sur.

Niche modeling and niche conservatism. Temperature and precipitation are the main factors determining the altitudinal and longitudinal plant distribution (^{Archibold 1995}). It has been proposed that precipitation and humidity variations have a more prominent effect on epiphytic plants (^{Hernández-Ruíz et al. 2016}, ^{Zotz 2016}). Even in the globular cactus Thelocactus, precipitation (precipitation in the wettest quarter) constrains the ENM for most species (^{Mosco 2017}). This pattern coincides for the ENM of Aporocactus flagelliformis, in which precipitation of the warmest month (BIO18) and the driest month (BIO14), and the temperature seasonality defined the model. Also, for A. martianus, the mean diurnal range (BIO2) and precipitation of the warmest month (BIO18) defined the model. Temperature seasonality is considered important in growth and other phenological processes (^{Menzel & Sparks 2006}). The latter factor is critical for the conservation of Aporocactus and other epiphytic cacti in the context of climate warming. Although other regions with high suitability of distribution for A. flagelliformis were recovered, it is necessary to corroborate their presence in particular zones (e.g., Sierra Madre Oriental at San Luis Potosí) or to investigate biological factors limiting the actual distribution (e.g., pollinator availability).

Analyses suggested that the niche overlap is low and niches of the two species of Aporocactus are not identical and are significantly differentiated. Species of Aporocactus have specific environmental constraints and do not occupy niches that are similar as possible given what is available. Epiphytic plants in cloud forests are especially sensitive to climate changes (^{Foster 2001}), floristic and climatic differences have been documented for cloud forests in Hidalgo, Querétaro, and central Veracruz versus cloud forests in southern Veracruz and Oaxaca (^{Ruíz-Jiménez et al. 2012}). Comparative analysis of niche overlap and niche similarity has been addressed in other close related Mexican plants and cacti, and lead some authors to consider the existence of niche conservatism on those lineages (^{Suárez-Mota et al. 2015}, ^{Mosco 2017}, ^{Gutiérrez-Ortega et al. 2020}), however our results suggest niche divergence in these sister species. A critical review by ^{Münkemüller et al. (2015)} suggests that studies investigating niche conservatism should compare alternative evolutionary models, including multiple-optima OU models. A comparative niche evolution analysis, as previous authors recommend, including a wider sampling of the tribe Hylocereeae, will allow to corroborate phylogenetic niche conservatism and niche shift in the Mesoamerican epiphytic lineages of cacti. For now, base on the difference of the ecological niches, we suggest the possibility of niche divergence in Aporocactus, as is expected for allopatric species (^{Peterson et al. 1999}, ^{Warren et al. 2008}). Finally, the primary differences between both species of Aporocactus are established by the floral morphology; therefore, it is likely that the primary factor driving the evolution of these lineages is their association with pollinators. For many years, epiphytic cactus species have received scarce attention. Although Aporocactus is a small genus, it may represent an interesting model for research on such topics as the ecology of pollination, population genetics, and flower development to characterize the evolution of those specialized cacti.

Supplementary material

Supplemental data for this article can be accessed here: https://doi.org/10.17129/botsci.2893

Supplementary material Figure S1, Figure S2, Table S1 and Table S2.