This study focused on the Magnolia pacifica species complex, which comprises six species (M. granbarrancae, M. pacifica, M. pugana, M. talpana, M. tarahumara, and M. vallartensis), as well as M. iltisiana (Figs 1, 2). Magnolia iltisiana was chosen as an outgroup due to its morphological similarity and geographical proximity to the species within the M. pacifica complex (Vázquez-García 1994).

Morphologically, Magnolia pacifica s. str. and M. talpana are similar. They share ellipsoid glabrous fruits and elliptic to elliptic-lanceolate leaves. However, they differ in flower size, with M. talpana having smaller flowers and a compact pollination chamber involving the outer and inner whorls of petals, where the inner petals are nearly equal in size, unlike the varying sizes in M. pacifica (Vázquez-García et al. 2021). Additionally, M. talpana has pubescent peduncles, whereas M. pacifica has glabrous peduncles. Both species inhabit cloud forests, riparian forests, and ecotones with oak-pine forests, with M. pacifica found in north-western Jalisco and southern Nayarit, and M. talpana endemic to western Jalisco (Fig. 2). Magnolia vallartensis is characterized by smaller glabrous fruits with fewer carpels, broadly oblong to elliptic leaves with an obtuse apex, creamy-white flowers, and a smaller tree size (up to 15 m) (Vázquez-García et al. 2012). It is found in riparian forests, tropical sub-deciduous and subperennial forests, occurring in ecotones with tropical cloud forests at lower elevations in the coastal mountains of Puerto Vallarta and Cabo Corrientes (Fig. 2).

Magnolia pugana and M. granbarrancae are morphologically similar, with both species having glabrous fruits with glabrous peduncles and lanceolate to narrowly elliptic leaves. Unlike M. pugana, M. granbarrancae has smaller flowers and a tighter pollination chamber. Additionally, M. granbarrancae has smaller subglobose to broadly obovoid fruits with orange seeds, while M. pugana has oblongoid to ellipsoid fruits with scarlet red seeds. Both species thrive in seasonally dry climates, with M. granbarrancae strictly inhabiting gallery riparian forests surrounded by tropical seasonal dry forest endemic to the Rio Verde region, central Jalisco, and M. pugana inhabiting streams and springs surrounded by tropical seasonal dry forest and oak-pine forests in northern Jalisco and southern Zacatecas (Fig. 2).

Magnolia tarahumara exhibits a gradient of morphological characteristics across its broad latitudinal geographic distribution range. Its leaves vary from elliptic to lanceolate-elliptic, and its glabrous fruits range from small and rounded with few carpels to elliptical-oblongoid with more carpels, rarely showing pubescence. This species inhabits cloud forests from northern Nayarit to moist canyons with perennial streams, characterized by riparian vegetation and oak forests, in the southern part of Chihuahua and Sonora throughout the Sierra Madre Occidental (Fig. 2).

Magnolia iltisiana features elliptic to oblong-lanceolate leaves that are green and glabrous above, and pale green and pubescent beneath (Vázquez-García 1994). Its flowers are creamy-white to pale yellowish, and its fruits are sub-globose to oblong and densely yellowish pubescent. This species inhabits moist sub-deciduous cloud forests, primarily in ravines, ranging from south Jalisco to Michoacan (Fig. 2).

Plant material was collected in 25 localities throughout the Magnolia pacifica species complex range (Fig. 2; Table 1). Samples were stored in Silica Gel at -20°C until DNA extraction. Genomic DNA was extracted using a DNeasy® Plant Mini Kit (QIAGEN), and resultant concentrations were quantified using a NanoDrop™ One UV-Vis spectrophotometer.

We performed PCR amplification of three chloroplasts (cpDNA) intergenic spacers: trnT-trnL, rpl32-trnL (Azuma et al. 2011; Rico and Gutiérrez Becerril 2019), and trnH-psbA (Azuma et al. 1999; Kim and Suh 2013). Additionally, three low-copy nuclear Conserved Ortholog Set (COS) genes were amplified: LEAFY (LFY), Phytochrome A (PHYA), and Alanine-glyoxylate aminotransferase (AGT1). As described in previous studies (Mathews and Donoghue 1999; Nie et al. 2008; Kikuchi and Osone 2020), universal primers were used to amplify COS genes, which have proven useful in phylogenetic reconstructions of closely related species (Li et al. 2008; Zimmer and Wen 2013; Kikuchi and Osone 2020).

Amplifications were conducted in 18 µL reaction volumes, containing 150 ng of genomic DNA, 1X PCR buffer, 2–2.5 mM MgCl₂, 0.2 mM dNTPs mix, 0.4 μM of each primer, and 0.8 U Taq polymerase (Invitrogen), using a PTC-100 thermal cycler (MJ Research, Inc.). PCR programs were tailored for each primer (Suppl. material 1). Amplification products were stained with GelRed (Biotium, Fremont, California) and separated by electrophoresis on 1.2% agarose gels in 1X TBE buffer, run at 80 V for approximately 1.5 h. The gels were then visualized and photographed under UV light using a D-Digit™ Gel Scanner (LI-COR, Bioscience).

Following visualization, PCR products were sequenced by the University of Arizona Genetics Core (Arizona, USA). Sequence quality was subsequently reviewed and edited using Sequencher v.4.1.4 (Gene Codes Corporation, Ann Arbor, MI, USA), while sequence alignments were manually constructed using PhyDE-1 (Müller et al. 2010). The resulting sequences were deposited in GenBank (Suppl. material 2).

For the nuclear COS genes, the genotypes were phased into haplotypes using DnaSP v.6.12.03 (Rozas et al. 2017) with the PHASE algorithm (Stephens et al. 2001; Stephens and Donnelly 2003). This process was applied to each gene using default MCMC settings, assuming no recombination within a gene. Additionally, we included extra sequences from GenBank (Suppl. material 3).

To determine the model of molecular evolution that best fits each marker, we applied corrected Akaike Information Criterion (AICc) in jModelTest v.2.1.6 (Darriba et al. 2012). Subsequently, a concatenated dataset comprising the cpDNA intergenic spacers and nDNACOS genes was conducted using Mesquite v.2.75 (Maddison and Maddison 2011).

To determine whether the Magnolia pacifica complex forms a monophyletic group within M. sect. Magnolia, a phylogenetic analysis was conducted including 13 representative taxa from the section and eight individuals from the M. pacifica complex. Two species from M. sect. Macrophylla were included as outgroup to root the tree (Suppl. materials 2, 3). To further analyse the relationship and classification of the M. pacifica complex species, a focused reconstruction was conducted, including 67 individuals representing all sampled populations (Suppl. material 2), with Magnolia macrophylla Michx. retained as the outgroup. The outgroup selection was based on studies by Nie et al. (2008) and Dong et al. (2022). We employed two analytical methods: Bayesian inference (BI) and maximum likelihood (ML).

The BI analyses were run on the CIPRES Science Gateway (Miller et al. 2010) using MrBayes v.3.2.7 (Ronquist et al. 2012). Each run consisted of four Markov Chain Monte Carlo (MCMC) of 20,000,000 generations, with tree sampled every 1,000 generations. The first 25% of the sampled trees were discarded as burn-in. Posterior probabilities (PP) were calculated using a 50% majority-rule consensus of the retained trees after burn-in. The harmonic means of likelihood were estimated using the ‘sump’ command (Ronquist and Huelsenbeck 2003). MCMC chains were evaluated for each independent run in Tracer v.1.7.1 (Rambaut et al. 2018) to verify that the effective sample sizes (ESS) were > 200 (Drummond and Bouckaert 2015). The 50% majority-rule consensus tree was visualized in FigTree v.1.4.3 (Rambaut 2016). Concurrently, ML analyses were performed using RAxML (Randomized Accelerated Maximum Likelihood) via raxmlGUI 2.0 software (Edler et al. 2020). The parameter function was set to “ML + rapid bootstrapping” with 1000 replicates.

To estimate divergence times, we constructed a ML tree based on 94 individuals (Suppl. materials 2, 3), including 78 ingroup and 16 outgroup samples. The ingroup comprised representative individuals of the Magnolia pacifica complex, including 23 individuals of M. tarahumara, 21 of M. pacifica s. str., four of M. talpana, 12 of M. pugana, eight of M. vallartensis, five of M. granbarrancae, and five of M. iltisiana. Additionally, 13 individuals from other taxa within M. sect. Magnolia were included as outgroups. Two taxa from M. sect. Macrophylla were also incorporated as outgroups, and M. mexicana DC. (M. sect. Talauma) was designated as the most external taxon to root the tree.

Divergence times were estimated along branches on the concatenated ML tree (comprising cpDNA intergenic spacers + nDNACOS genes) using the RelTime-ML algorithm in MEGA v11.0. This algorithm employs ML-optimized branch-length estimation and represents a robust statistical approach for divergence time inference (Mello 2018). To calibrate our findings, we utilized a normal distribution with a mean age of 32 Ma and a specified standard deviation for the root node of M. sect. Macrophylla. This calibration was based on the calibrated multi-species coalescent summary tree for Magnoliaceae by Veltjen et al. (2022). The summary tree was visualized using FigTree v.1.4.3.

Species delimitation within the Magnolia pacifica complex was assessed under the multispecies coalescent (MSC) framework using Bayesian Phylogenetics and Phylogeography (BPP) v.4.2 (Yang 2015). BPP uses the MSC model to compare alternative species delimitation methods within a Bayesian framework using molecular data, without relying on reciprocal monophyly as a criterion for species recognition (Yang and Rannala 2010; Rannala and Yang 2013). This approach is particularly suitable for recently diverged lineages and gene tree discordance.

Species delimitation and species tree estimation were first conducted using algorithm A11, which jointly estimates the number of species and their phylogenetic relationships under a reversible-jump Markov Chain Monte Carlo (rjMCMC) framework. Analyses were based on 67 individuals of the M. pacifica complex (Suppl. material 2) and six loci, including three plastid intergenic spacers and three low-copy nuclear COS genes. Nuclear COS genes were treated as diploid, and phased haplotypes were used as input, while cpDNA loci were modelled as haploid. Locus-specific inheritance scalars were set to 1.0 for nuclear genes and 0.25 for plastid markers (Hey and Nielsen 2004). The species tree topology was allowed to vary during MCMC using nearest-neighbour interchange (NNI) and subtree-pruning-regrafting (SPR) moves, thereby enabling exploration of alternative species delimitation models and topologies without assuming species monophyly a priori.

To further evaluate support for the species tree inferred under A11, a second analysis was conducted using algorithm A10, in which the species tree topology with the highest posterior probability obtained under the A11 analysis was fixed. Under this framework, rjMCMC moves were restricted to species delimitation models compatible with the fixed topology, while maintaining the same prior distributions for population size and divergence time parameters. This approach allowed for the estimation of posterior probabilities for internal nodes and lineage relationships conditional on the inferred species tree, providing an explicit assessment of support for relationships among delimited lineages.

For both phylogenetic reconstructions, the molecular evolution models produced by jModelTest for cpDNA intergenic regions were: trnH-psbA = GTR + G, trnT-trnL = TPM1uf + G, rpl32-trnL = TVM + I; meanwhile, for the COS genes, models were: LFY = TPM1 + G, PHYA = TIM1ef + G, and AGT1 = HKY + I.

The final dataset used to test the monophyly of the Magnolia pacifica complex within M. sect. Magnolia included 23 samples; the ingroup comprised 13 representative taxa of M. sect. Magnolia and eight individuals of the M. pacifica complex, while the outgroup was composed of M. macrophylla and M. dealbata Zucc. from M. sect. Macrophylla. The number of aligned base pairs (bp) was 433 bp for trnH-psbA, 809 bp for trnT-trnL, 700 bp for rpl32-trnL, 797 bp for LFY, 965 bp for PHYA, and 1,109 bp for AGT1. The concatenated cpDNA + nDNA matrix was 4,818 bp in length and contained 413 polymorphic sites. Partitioned molecular reconstruction of the concatenated cpDNA + nDNA matrix using Bayesian inference (BI) and maximum likelihood (ML) retrieved similar topologies. The ML tree is depicted in Fig. 3. The phylogenetic tree showed M. sect. Magnolia is a monophyletic group rooted by Magnolia dealbata and M. macrophylla (both belonging to M. sect. Macrophylla). Within M. sect. Magnolia, the M. pacifica complex forms a monophyletic group with high posterior probabilities (PP = 1) and parametric bootstrap support (BS = 99) values.

For the focused phylogenetic reconstruction of the M. pacifica complex, the final dataset consisted of 68 samples, including 67 ingroup and M. macrophylla as outgroup. The number of aligned base pairs (bp) was 430 bp for trnH-psbA, 788 bp for trnT-trnL, 702 bp for rpl32-trnL, 797 bp for LFY, 905 bp for PHYA, and 906 bp for AGT1. The concatenated cpDNA + nDNA matrix was 4,533 bp in length and contained 185 polymorphic sites. The partitioned molecular reconstruction of the concatenated cpDNA + nDNA matrix using ML and BI retrieved similar topologies. The ML tree is depicted in Fig. 4. The phylogenetic analysis recovered a moderately supported tree with four clades. Clade 1 (shown in yellow in Fig. 4): composed of seven out of eight M. vallartensis individuals, one individual of M. pacifica, and one individual of M. tarahumara. Clade 2 (shown in red in Fig. 4): composed of all individuals of M. pugana and M. granbarrancae. Clade 3 (shown in blue in Fig. 4): composed of all individuals of M. pacifica, M. iltisiana, and all individuals of M. talpana. Clade 4 (shown in green in Fig. 4): composed of 15 out of 23 individuals of M. tarahumara. The remaining individuals are outside the four clades.

Based on the concatenated matrix (cpDNA intergenic spacers + nDNACOS genes) and the calibration applied in this study, the split of M. sect. Magnolia was estimated to have occurred around 25 Mya during the late Oligocene. The divergence of the species within the section was estimated at ~11 Mya (95% confidence interval: 9.75–11.31 Mya) during the late Miocene. Furthermore, the M. pacifica complex is grouped within M. sect. Magnolia with a split at ~7.5 Mya and a very recent divergence date of 1.27 Mya (95% confidence interval: 0.64–2.53 Mya) during the early Pleistocene (Fig. 5).

Species delimitation analysis using the A11 algorithm in BPP (under the multispecies coalescent), strongly supports a model recognizing four independently evolving lineages: corresponding to Magnolia tarahumara, M. iltisiana – M. pacifica (i.e. M. iltisiana + M. pacifica + M. talpana), M. pugana (i.e. M. pugana + M. granbarrancae), and M. vallartensis. The most frequently sampled model (PP = 0.738) recovered a species tree topology where M. tarahumara is sister to a clade comprising M. iltisiana – M. pacifica, as well as M. pugana and M. vallartensis. Models recognizing fewer than four species received negligible support (PP ≤ 0.004), providing strong evidence against lineage collapse within the complex.

Complementarily, the A10 analysis, which evaluates species delimitation on a fixed guide tree, provided unequivocal support for the same four-species hypothesis (Fig. 6). The model in which all three internal nodes represent speciation events received a posterior probability of 1.0, while all alternative models were rejected.

The guide tree recovered M. tarahumara as the most divergent lineage, with M. pugana and M. vallartensis forming a strongly supported clade sister to M. iltisiana – M. pacifica, and all nodes showing posterior probabilities of 1.0. Together, these results indicate that, despite limited resolution in concatenated phylogenetic reconstructions, multilocus coalescent-based analyses consistently support the recognition of four distinct evolutionary lineages within the Magnolia pacifica complex.

Here, we present the first phylogenetic analysis of all described species in the M. pacifica complex (M. granbarrancae, M. pacifica, M. pugana, M. talpana, M. tarahumara, and M. vallartensis). These species, along with M. iltisiana, were recovered as a monophyletic group within M. sect. Magnolia. This result is in line with previous phylogenetic and phylogenomic analyses involving certain species of the M. pacifica complex (Azuma et al. 2001; Kim et al. 2001; Nie et al. 2008; Wang et al. 2020; Rico et al. 2021; Veltjen et al. 2022; Aldaba Núñez et al. 2024). The placement of M. iltisiana within the M. pacifica complex is supported by its close relationship, as evidenced in previous phylogenetic reconstructions (Azuma et al. 2001; Kim et al. 2001; Nie et al. 2008; Veltjen et al. 2022; Aldaba Núñez et al. 2024, 2026). Moreover, these species exhibit morphological similarities (Fig. 1) and geographical proximity (Fig. 2).

The phylogenetic relationships within the complex remain somewhat inconclusive in this analysis. Magnolia pacifica, M. talpana, M. tarahumara, and M. vallartensis are polyphyletic, potentially due to introgression or incomplete lineage sorting. The inclusion of additional nDNA and cpDNA regions in future analyses could help further refine the phylogeny.

In our study, we estimated the split of M. sect. Magnolia occurred around 25 Mya, with the divergence within the section happening approximately 11 Mya. These estimated times are consistent with previous research, which suggested an earlier split of M. sect. Magnolia around 30 Mya in the late Oligocene, and a later divergence of M. sect. Magnolia at approximately 10 Mya during the late Miocene (Dong et al. 2022; Veltjen et al. 2022; Guzmán-Díaz et al. 2025). Furthermore, the divergence of the M. pacifica complex from the rest of M. sect. Magnolia is estimated to have occurred around 7.5 Mya, followed by rapid radiation around 1.27 Mya during the Calabrian age of the Pleistocene, the most recent period of repeated glaciations (Hewitt 2004). The Pleistocene climate fluctuations, combined with the highly heterogeneous mountainous environments in western and northwestern Mexico, likely created climatic and orographic barriers, promoting isolation and hindering gene flow among Magnolia populations (Gugger et al. 2011; Ornelas et al. 2013). This scenario could have facilitated the allopatric speciation of the M. pacifica complex. For instance, previous studies have suggested that climatic events have influenced the geographical distribution of species like M. pedrazae A.Vázquez and M. schiedeana Schltdl., leading to range shifts, expansions, and contractions (Rico et al. 2021). This hypothesis is supported by various phylogeographic studies in different plant species (Ruiz-Sánchez and Ornelas 2014; Ornelas et al. 2019; Ortiz-Rodriguez et al. 2020). Further research should consider incorporating population genetics processes (Maddison and Knowles 2006), such as those assessed by phylogeography, as well as incorporating additional nDNA and cpDNA regions or even phylogenomics. These approaches could provide further insights into the biogeographic and climatic factors influencing the speciation processes within the complex.

The Magnolia pacifica complex has undergone a complex and dynamic taxonomic history with repeated changes in species circumscription based on morphological, geographic, and ecological characteristics. Originally, M. pacifica was treated as a single, morphologically variable species comprising multiple subspecies with broad, partially overlapping geographic distributions, and closely associated with M. iltisiana due to their morphological similarity and geographic proximity (Vázquez-García 1994). Subsequent taxonomic revisions progressively elevated several of these lineages to species rank: M. pugana, M. tarahumara, M. vallartensis, M. talpana, and M. granbarrancae, based primarily on combinations of vegetative, floral, and fruit characters, as well as ecological and geographic differentiation without a formal statistical framework (Vázquez-García et al. 2002, 2012, 2013, 2021). However, the recognition of these taxa has remained challenging due to overlapping morphological variation and the recent divergence of lineages, raising long-standing questions regarding species boundaries and evolutionary independence within the complex. This taxonomic uncertainty provided the central motivation for applying a multispecies coalescent framework to explicitly test species limits and evolutionary relationships within the M. pacifica species complex using multilocus molecular data.

Under the multispecies coalescent model in BPP, we recovered four independently evolving lineages: Magnolia iltisiana – M. pacifica (i.e. M. iltisiana + M. pacifica + M. talpana), M. pugana (i.e. M. pugana + M. granbarrancae), M. tarahumara, and M. vallartensis, reflecting shared ancestry and distinct coalescent histories among populations.

The Magnolia iltisiana – M. pacifica lineage, includes individuals currently assigned to M. iltisiana, M. pacifica, and M. talpana. These populations are distributed along a latitudinal continuum across the mountainous regions of the Western Sierra Madre del Sur subprovince (Morrone 2017) and partially share ecological conditions, suggesting historical or ongoing gene flow among populations. Notably, M. talpana (TAL locality; Table 1) was originally included within M. pacifica (Vázquez-García 1994), and its subsequent recognition as a distinct species was based on consistent floral differences, including reduced flower size, a more compact inner petal whorl, and a differentiated pollination chamber (Vázquez-García et al. 2021). However, previous molecular evidence from six nuclear ISSR markers revealed only weak genetic differentiation between M. talpana and M. pacifica s. str. (Muñiz-Castro et al. 2020). Consistent with previous molecular evidence, our multispecies coalescent analyses support the interpretation that these taxa form a single, independently evolving lineage. The strong support for a single coalescent lineage (PP = 1.0 in the A10 analysis) indicates that, despite their phenotypic divergence, these populations have either diverged too recently for ancestral polymorphism to have sorted completely across the six loci assayed or that low levels of gene flow have been sufficient to maintain genetic cohesion. However, it is possible that we are underrepresenting M. iltisiana, as we included individuals from only one population, located close to M. pacifica. Magnolia iltisiana differs morphologically from M. pacifica in having densely pubescent bracts, peduncular internodes, petioles, and flower buds, and yellowish-greenish velutinous to felted fruits, in contrast to the glabrous to rarely pubescent and usually completely glabrous fruits of M. pacifica (Vázquez-García 1994). Following Morrone (2017), these species inhabit different biogeographic districts within the Western Sierra Madre del Sur subprovince. While M. pacifica occurs in the Jaliscan-Tuito district, M. iltisiana inhabits the Jaliscan-Manantlan district. This biogeographic pattern has been observed in other plant species: Abies jaliscana (Martínez) Mantilla, Shalisko & A.Vázquez, Daphnopsis mexiae Nevling, Juniperus jaliscana Martínez, Pinus jaliscana Pérez de la Rosa, Quercus cualensis L.M.González, Q. tuitensis L.M.González, occurring in the Jaliscan-Tuito district; and Abies flinckii Rushforth., A. religiosa (Kunth) Schltdl. & Cham., Beilschmiedia manantlanensis Cuevas & Cochrane, and Arachnothryx manantlanensis (Lorence) Borhidi in the Jaliscan-Manantlan district. This transition zone is considered a potential ecological barrier between adjacent areas due to differences in climatic conditions such as humidity and seasonal precipitation (Vázquez-García et al. 2010, 2014; Santiago-Alvarado et al. 2016). Their biological and taxonomic significance will require further evaluation using additional molecular data, as well as detailed morphological, geographic, and ecological evidence analysed within a statistical framework. Therefore, our results suggest the synonymization of M. talpana under M. pacifica, while M. iltisiana is herein maintained as a distinct species.

The M. pugana lineage includes individuals currently assigned to M. pugana and M. granbarrancae. The population of M. granbarrancae was previously considered part of M. pugana. Evidence from six nuclear ISSRs showed genetic differentiation among allopatric populations of M. pugana and revealed two subgroups (Muñiz-Castro et al. 2020). A portion of one subgroup (locality RV) was later described as M. granbarrancae based on morphological differences such as smaller flowers and fruits, fewer carpels, and orange seeds vs scarlet red. These distinctions led to the description of M. granbarrancae as a distinct species from M. pugana (Vázquez-García et al. 2021). However, our multispecies coalescent analysis reveals that this morphological and moderate genetic differentiation (as detected by dominant ISSR markers) does not correspond to a break in coalescent history. The model that recognizes M. pugana and M. granbarrancae as a single independently evolving lineage received unequivocal support (PP = 1.0, A10 analysis). This suggests that the observed differentiation likely represents either: 1) incipient divergence with insufficient time for complete lineage sorting across the loci sampled, or 2) phenotypic divergence driven by local adaptation. Consequently, our genetic data do not support the recognition of M. granbarrancae as a distinct lineage from M. pugana. Based on these results, we propose the synonymization of M. granbarrancae under M. pugana.

We confirmed that M. vallartensis represents a separately evolving lineage, consistent with the morphological distinctiveness of Magnolia populations inhabiting the coastal mountains of the municipality of Puerto Vallarta (Vázquez-García et al. 2012). It is distinguished from M. pacifica by a combination of morphological characters: outer petals with a colourful abaxial base; leaves broadly oblong to elliptic (vs elliptic-lanceolate) and larger in size; fewer carpels; and more numerous petals. Furthermore, ecological evidence suggests that the tropical river canyons separating M. vallartensis and M. pacifica populations act as barriers to effective gene flow. Although birds and beetles can potentially disperse pollen and seeds across these canyons, the microclimatic conditions at lower elevations – characterized by tropical deciduous forest rather than cloud or evergreen forests – prevent seedling establishment. This creates an ecological filter that, despite long-distance dispersal events, results in genetic isolation and differentiation between the two species (Muñiz-Castro et al. 2020). These genetic differences may also be attributed to parapatric ecological differentiation (Nosil 2012), given that M. vallartensis and M. pacifica occur in distinct habitats. Magnolia vallartensis thrives in the ecotone of tropical montane cloud forest and tropical sub-evergreen forest at elevations ranging from 100 to 1,120 m a.s.l. experiencing a tropical maritime climate with higher relative humidity, mean annual temperature, and rainfall. In contrast, M. pacifica grows in montane cloud forests at elevations of 750 to 2,250 m a.s.l., in a temperate montane climate characterized by lower relative humidity, mean annual temperature, and precipitation. Additionally, these species differ in their flowering phenology: M. vallartensis exhibits no evident flowering seasonality, whereas M. pacifica shows a highly seasonal pattern (Dahua-Machoa 2018; Muñiz-Castro et al. 2020).

Our coalescent delimitation strongly supports M. tarahumara as a distinct evolutionary lineage (PP = 1.0, A10 analysis), despite its paraphyletic presentation in the concatenated phylogenetic tree. This apparent contradiction is a known phenomenon in species delimitation and can be explained by the effectiveness of multispecies coalescent models to distinguish between incomplete lineage sorting (ILS) and true evolutionary independence (Jiao et al. 2021). Magnolia tarahumara exhibits a broad latitudinal distribution spanning approximately 700 km along the Sierra Madre Occidental. It inhabits a variety of environments and elevations, ranging from cool, moist canyons with perennial streams and nearly evergreen riparian vegetation at lower elevations in the Pacific Low Lands (850 to 1200 m a.s.l.) in northern localities, to cloud forests at higher elevations on the Sierra Madre Occidental (1400 to 2100 m a.s.l.) in the southern localities (Vázquez-García 1994). This lineage also displays a gradient of morphological traits, with small, rounded fruits and few carpels in the northern localities, transitioning to elliptical-oblongoid fruits with more numerous carpels in the southern localities. The southernmost part of its range, in northern Nayarit, is morphologically recognized as the “Huajicori group” (Vázquez-García et al. 2021), which appears genetically divergent in the concatenated phylogenetic analysis. The Huajicori group could represent an incipient divergence or a locally adapted population, but it has not yet crossed the threshold of evolutionary independence recognized by the coalescent model. Therefore, our data validate M. tarahumara as a separately evolving lineage.

Our species delimitation hypothesis results should be interpreted with caution considering both intrinsic characteristics regarding speciation (1–3), and methodological concerns (4–6): 1) speciation is not instantaneous, but rather a continuum process (De Queiroz 2007), 2) morphological characters may evolve faster than the formation of species, or may not accompany genetic divergence (Chan and Grismer 2019), 3) a divergence time as recent as 1.27 Mya may not be sufficient for the formation of separately evolving lineages (Hudson and Coyne 2002; Knowles and Carstens 2007), 4) small sample size can impact the accuracy of the BPP results, warranting cautious interpretations (Carstens et al. 2013; Chan and Grismer 2019), 5) it is argued that the MSC model reflects genetic structure, which includes both population structure within species as well as between species (Sukumaran et al. 2021), 6) BPP has been shown to over-delimit species (Leaché and Fujita 2010; Sukumaran et al. 2021). De Queiroz (2007) introduced the concept of a “grey zone” in which lineage separation and divergence are not yet fully complete, leading to conflicts in species concepts and making species boundaries difficult to define. Despite these important considerations, our strong statistical support from the MSC analyses provides a robust and testable hypothesis for species boundaries within the Magnolia pacifica complex. We therefore propose the recognition of five independently evolving lineages as the most coherent taxonomic framework, given the current morphological and molecular data: Magnolia pacifica (encompassing M. pacifica and M. talpana), M. iltisiana, M. pugana (encompassing M. pugana and M. granbarrancae), M. tarahumara, and M. vallartensis. This species delimitation has direct implications for conservation, as it redefines the evolutionary units to prioritize, ensuring that management efforts align with the genetic diversity of the group. Future studies employing independent lines of evidence, such as genomic-scale datasets, morphometric, ecological, expanded sampling, and phylogeographic approaches, will be crucial to further test this hypothesis, and to elucidate the evolutionary dynamics within this fascinating group.