The leaf samples for the morphological and genetic analyses were collected between 2015 and 2020 from individuals representing of Magnolia sect. Talauma subsect. Talauma throughout their entire distribution range in the mountains of Nipe-Sagua-Baracoa and Sierra Maestra in eastern Cuba (Fig. 2). In the field, individuals were identified based on tree and leaf shape according to the morphological criteria outlined by Bisse (1988), because this author defined the highest number of species units (Fig. 1). León and Alain (1951) considered T. truncata an independent species. However, in the present work, individuals that could have been considered as T. truncata, were considered as part of the variability of M. orbiculata, as has been recognized by Bisse (1988), Imkhanitzkaja (1993), and Palmarola et al. (2016). To confirm species identity, 43 herbarium vouchers were collected or reviewed (Table 1). All herbarium vouchers were deposited in the Herbarium Johannes Bisse (HAJB, herbarium acronyms follow Thiers 2022) at the National Botanic Garden (University of Havana). The number of samples per species and localities is shown in Table 1.

Figure 2. 

Geographic distribution of sampling locations of Magnolia subsect. Talauma in Cuba.

For the morphological analyses, 4–8 healthy leaves from 200 individuals were randomly collected, across the entire range of taxa within each locality. A leaf was considered healthy if the full outline of the leaf was undamaged. Leaves were photographed with a Nikon camera on a white background with a fixed ruler. The petiole of the leaf was removed before taking pictures, and the camera was mounted on a tripod to standardize the angle and distance of the photographs. To expand the geographic scope of our study, we also included leaf samples from 43 herbarium specimens (deposited in HAC, HAJB, and B). Hence, in total 243 individuals of Magnolia sect. Talauma subsect. Talauma in Cuba were morphologically analyzed.

For the genetic analyses, young leaf samples of a total of 461 individuals, belonging to 26 of 30 known localities, were stored in self-sealed bags with silica gel for DNA extraction. The resulting number of DNA samples represented 52% of the known individuals of Magnolia subsect. Talauma in Cuba (close to 900 individuals).

Analyses based on morphological variables were aimed at comparing the relevance of each of the three CS previously proposed: the two taxa CS, Magnolia minor and M. orbiculata, of León and Alain (1951); the four taxa CS, M. minor, M. orbiculata, M. oblongifolia, and Talauma ophiticola, of Bisse (1988); and the most recent, the three taxa CS, M. minor, M. orbiculata, and M. oblongifolia, of Palmarola et al. (2016) (Fig. 1). Two types of morphological analyses were carried out on three independent datasets: 1) multivariate morphometry analysis: a linear and angular measures dataset, 2) geometric morphometry analysis: an outline dataset and a landmarks coordinates dataset.

In the multivariate morphometry analysis, linear and angular measures of leaf characters were automatically taken from the digital photographs using the R v. 3.4.1 (R Core Team 2017) package FOLIOMETRIK v. 0.2.2 (Ramírez-Arrieta and Denis 2020). Eleven leaf variables were measured: central axis length (Length), maximum width, width at the three main quartiles (25, 50, and 75 quartiles), the perimeter of the contour (Perimeter), surface area (Area), and internal angles (v1 = angle of the base, v2 = angle of the apex; m1 and m2 = lateral angles at the maximum width) (Fig. 3A). Additional to the eleven measured variables, we calculated the maximum width/length ratio, named Calculated Index of Bisse (B_ci) for each leaf. The eleven variables were recorded for each leaf. Subsequently, the twelve variables were averaged per individual for the 4–8 leaves available per individual. These averages of the twelve variables were used for all the subsequent statistical analyses.

Figure 3. 

A the 11 morphological variables measured on leaves of Magnolia subsect. Talauma in Cuba; v1 = angle of the base, v2 = angle of the apex; m1 and m2 = lateral angles in the maximum wide B quadratic grid with six lines and the position of the 14 landmarks (type 1: points 1 and 8; type 2: the other 12 points), placed on the leaves of Magnolia subsect. Talauma in Cuba.

In the geometric morphometry analysis, the outline dataset was obtained through a semi-automated shape analysis performed in FOLIOMETRIK v. 0.2.2 (Ramírez-Arrieta and Denis 2020). We set the program outputs to the Elliptic Fourier Descriptors (EFDs) (Jensen 2003), to obtain the first 25 harmonics (Chuanromanee et al. 2019). The harmonics were normalized to eliminate the differences in size, position, rotation, and starting point. This allowed removing the undesired experimental source of random variation and analyzing true differences of leaf shape between individual measurements (Jensen 2003). The landmarks coordinates dataset was obtained as follows. The positions of landmarks were determined by placing a quadratic grid with six lines on each leaf. In between the intersections of the grid and the border of the leaf we set 14 landmarks, two of them anatomical (type 1, i.e. apex and base) and the other 12 mathematically defined (type 2) (Fig. 3B). All analyses were carried out in FOLIOMETRIK v. 0.2.2. The landmarks X and Y coordinates were standardized using a Generalized Procrustes Analysis in PAST v. 2.14 (Hammer et al. 2001). Next, two variables: the Sum EDMA (Euclidean Distance Matrix Analysis) and centroid size were calculated.

The statistical significance of the differences among taxa for each measured variable (linear and angular variables, Sum EDMA, and centroid size) was assessed by a MonteCarlo analysis in PopTools v. 3.23 (Hood 2010) with 10 000 random permutations. The variability of the whole sample was described by using a normalized Principal Component Analysis (PCA). Differences among taxa were tested according to a one-way nonparametric MANOVA, using Euclidian distance, with 10 000 randomizations. Correction of p-values for multiple testing was done using the Bonferroni method. The multivariate comparisons were done independently for each dataset (linear and angular measures dataset, outline datasets, and landmarks coordinate dataset). All statistical analyses were conducted in R v. 3.4.1 (R Core Team 2017) and PAST v. 2.14 (Hammer et al. 2001), and the threshold used to decide for statistical significance was a p-value of 0.001.

A Bayesian clustering approach based on Gaussian finite mixture models was carried out using each of the three datasets of morphological variables using the “mclust” R package (Scrucca et al. 2016). The method tests the number of clusters and different mixture models that best fit the data according to the number of clusters (G) chosen a priori. The method allows comparing the quality of the discrimination among clusters based on the Bayesian Information Criteria (BIC) allowing to choose the best value(s) of G, without any information about individual assignation to the different clusters. The default “mclust” setting was used to assess the 14 types of models which all differ in the covariance matrix landscape (see Scrucca et al. 2016; Zhang and Di 2020, for further details about the models). We varied G values between 1 and 9 (default option). Three independent analyses were performed using the three datasets; to compare the power of those three groups of variables to discriminate among taxa. Because the analyses sometimes provided clusters with only one individual, those clusters were considered “ghost” clusters and not considered as true clusters.

DNA was extracted from dried leaf tissue using a modified cetyltrimethylammonium bromide (CTAB) extraction protocol (Doyle and Doyle 1990) with MagAttract Suspension G solution mediated cleaning (Xin and Chen 2012). DNA quality was assessed using a spectrophotometer NanoDrop 1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Individuals were genotyped using 21 microsatellite markers (simple sequence repeats, SSR) (Suppl. material 8) developed on four Neotropical Magnolia species: M. lacandonica A. Vázquez, Pérez-Farr. and Mart.-Camilo (MA39), M. mayae Vázquez and Pérez-Farrera (MA40), M. dealbata Zucc. (MA41) and M. cubensis subsp. acunae Imkhan. (MA42) (Veltjen et al. 2019), using four-primer PCR multiplex method (Vartia et al. 2014). PCR conditions and primer labeling followed Veltjen et al. (2019). The combination and parameters of the four multiplex reactions are given in Suppl. material 8. The lengths of the DNA fragments were detected using an ABI 3130XL fragment analyzer, quantified with a GeneScanTM 500 LIZ size standard (Thermo Fisher Scientific), and analyzed in Geneious v. 8.0.5 (Kearse et al. 2012) with the microsatellite plugin.

Genetic diversity values were calculated for each taxa using GeneAlex v. 6.5 (Peakall and Smouse 2012) and Genepop v. 4.7.5 (Rousset 2008). Genetic differentiation between taxa was estimated through pairwise comparisons of F_ST (Weir and Cockerham 1984) and D_JOST values (Jost 2008) using the fastDivPart function of the R package diveRsity (Keenan et al. 2013). The identification of genetic clusters and the assignment of individuals was performed using STRUCTURE v. 2.3.4 (Pritchard et al. 2000), which uses a Bayesian clustering approach using MCMC for posterior distribution sampling. STRUCTURE analyses were conducted using a model that assumes admixture, correlated allele frequencies, and without prior population information. First, 10 replicates were run for each genetic clusters (K), with K varying between 1 to 20 and a burn-in period of 50 000 iterations followed by a run-length of 150 000 iterations of the Markov Chain. The most probable number of groups was determined according to the method of Evanno et al. (2005) as implemented in STRUCTURE HARVESTER v. 0.6.94 (http://taylor0.biology.ucla.edu/structureHarvester) (Earl and von Holdt 2012). Then, 100 new repetitions of the MCMC method were run for the best K value. CLUMPP v. 1.1.2 (Jakobsson and Rosenberg 2007) was used to estimate similarities between runs and to average the membership probabilities. Final bar plots displaying individual admixture coefficients were obtained thanks to Structure Plot v. 2.0 (Ramasamy et al. 2014). An individual was considered a member of a genetic group when its probability of belonging to that group was higher than or equal to 0.9. A second STRUCTURE analysis was executed (using the same configuration) without considering the individuals of M. oblongifolia (sensu Bisse 1988).

Because the MCMC method implemented in STRUCTURE is based on a population genetic model, the results of genetic clusters and assignment of individuals, may be affected by the potential low model fit to data. Thus, a non-model-based multivariate clustering analysis was also performed. A DAPC analysis (Discriminant Analysis of Principal Components) was executed in R v. 3.6.1 (R Core Team 2017) using the adegenet R package (Jombart et al. 2010). Firstly, a PCA was run on the whole dataset for which the first 200 Principal Components (PCs) were retained. Secondly, a discriminant analysis was executed using the number of genetic clusters defined in the previous step. Parallel to the STRUCTURE analysis, a second DAPC analysis was done without M. oblongifolia (sensu Bisse 1988). An individual was considered a member of a genetic group when its probability of belonging to that group was higher than or equal to 0.9.

For all analyses, graphical representations of outputs were built using the four taxa CS to have a representative overview of the correspondence between genetic clusters and each already defined taxon.

Because 138 individuals were analyzed both at the morphological and the genetic level, the correspondence between the groups inferred from both type of characters was assessed. The distributions of individual assignment to each morphological (mclust) and genetic (STRUCTURE) clusters were compared with a Chi² test carried out on the Contingency Assignment Table using PAST v. 2.14 (Hammer et al. 2001). A heatmap was made, in R v. 3.4.1, to analyze the variation in the cluster assignation inside each taxon.

The results of the multivariate morphometry analysis are summarized in Figs 4, 5, Suppl. materials 1–3, 9. The Calculated Index of Bisse (B_ci) showed significant differences between the defined species whatever the CS tested (p < 0.0001). Magnolia orbiculata and Talauma ophiticola displayed the highest and the lowest mean B_ci, respectively. Most of the eleven other variables showed significant differences between taxa, whatever the CS (Suppl. materials 1–3). There were three exceptions. The leaf perimeter did not show significant differences between Magnolia minor and Talauma ophiticola (p = 0.211), and between M. oblongifolia and T. ophiticola (p = 0.132) when the four taxa CS was considered. Likewise, the leaf area between M. minor and M. oblongifolia (p = 0.115) did not show differences when the three taxa CS was used. When following the two taxa CS (Fig. 4A, Suppl. material 1), eight variables (Maximum width, B_ci, Width-quartiles 50, Width-quartiles 75, Internal angles-v1, Internal angles-m1, Internal angles-v2, Internal angles-m1) showed an intra-taxon bimodal pattern within M. minor.

Figure 4. 

Graphic representation of Calculated Index of Bisse (B_CI) measured in the individuals of Magnolia subsect. Talauma in Cuba. Two taxa CS (A), four taxa CS (B), three taxa CS (C).

Despite the clear morphological differentiation between taxa, overlap in the multivariate distributions of leaf morphology variables was observed (Fig. 5). The internal angle of the base (-0.324) and the leaf perimeter (0.5627) displayed the highest weight in the first two principal components, respectively (Suppl. material 9). The NPMANOVA showed significant statistical differences (p < 0.0001) between taxa for each of the CS (Suppl. material 10). The comparison between groups, based on Sum EDMA and centroid size, showed significant differences for most comparisons (Suppl. material 4). The exceptions were: the Sum EDMA between M. minor and M. oblongifolia (p = 0.316) and between M. orbiculata and T. ophiticola (p = 0.406), when referring to the four taxa CS (Suppl. material 4). Fig. 6 illustrates PCAs on the outline dataset (Fig. 6A, C, E), and the landmark dataset (Fig. 6B, D, F) for the two (Fig. 6A, B), three (Fig. 6C, D) and four (Fig. 6E, F) taxa CS. Based on PCA for elliptic Fourier descriptors and Landmark, the different taxa had little overlap in the ordination space (Fig. 6). However, a clearer distinction among taxa was obtained with landmark positions than with other quantitative variables. This was especially obvious with M. orbiculata, which was strongly differentiated from other taxa when using landmark positions, no matter the CS considered.

The NPMANOVA showed significant statistical differences (p < 0.001) between taxa for each of the CS in the linear and angular measures dataset (Suppl. material 10). Similarly, the NPMANOVA showed significant statistical differences (p < 0.001) between the groups in the outline and landmark datasets (Suppl. material 10).

Figure 5. 

Principal Component Analysis for the multivariate morphometric variables measured in the individuals of Magnolia subsect. Talauma in Cuba. Two taxa CS (A), three taxa CS (B), four taxa CS (C).

Figure 6. 

Principal Component Analysis for the Elliptic Fourier Descriptors (A, C, E) and Coordinates of the landmark (B, D, F) which characterized the leaves of Magnolia subsect. Talauma in Cuba. Two taxa CS (A, B), three taxa CS (C, B), four taxa CS (E, F).

The clustering analysis based on morphological variability showed differences in the number of groups inferred by the best models, according to the different datasets (Fig. 7; Suppl. material 11). The highest BIC scores were retrieved for G = 4 for linear and angular dataset, G = 2 for the Elliptic Fourier Descriptors dataset (with other 3 ghost clusters), and G = 6 for the Landmark dataset (with other ghost clusters) (Fig. 7; Suppl. material 11). It was noticeable that for each data set, the probabilities of assignment of each individual were higher than 0.9 in all cases based on the Elliptic Fourier Descriptors. In the case of the linear and angular variables and matrix of landmarks, only 5 and 22 individuals showed probabilities of an assignment less than 0.9, respectively (data not shown). The linear and angular variables allowed a clear discrimination between M. orbiculata, T. ophiticola and M. minor, the latter taxa being split into two clusters. One of these two clusters was shared only with the majority of M. oblongifolia individuals.

Figure 7. 

Graphic representation and classification matrix obtained after the cluster analysis using the morphological data of Magnolia subsect. Talauma in Cuba. A Linear and angular variables B elliptic Fourier Descriptors C matrix of landmarks. * Ghost Cluster.

The clustering analysis based on Elliptic Fourier Descriptors provided only two clusters (Fig. 7). The assignment of individuals was therefore different from that obtained with linear and angular variables. Indeed, all individuals of M. oblongifolia, and most individuals of T. ophiticola and M. minor, were assigned to the same cluster (cluster 1), while most individuals of M. orbiculata were assigned to a different cluster (cluster 4). Therefore, Elliptic Fourier descriptors were efficient to discriminate between M. orbiculata on the one hand and the 3 other taxa on the other hand. Finally, the analysis carried out on the matrix of landmarks showed a similar pattern to that obtained with the linear and angular variables for M. minor, T. ophiticola and M. orbiculata. The main difference between these two analyses (matrix of landmarks and linear and angular dataset) was that in the first one, M. oblongifolia was split into two clusters, one of which was shared with M. minor and the other one with T. ophiticola (Fig. 7).

Thus, despite a continuous variation of leaf morphology across taxa, a clear delimitation of M. orbiculata is shown by our analyses whichever data set was used. In cases where individuals of the same taxon were assigned to different clusters, no obvious correspondence between the assigned clusters and the geographic origin of those individuals was found. Indeed, many individuals of the same taxon/locality were assigned to different clusters (data not shown).

The species with the greatest genetic diversity were Magnolia minor and Talauma ophiticola, while the lowest diversity was found in Magnolia orbiculata. The expected heterozygosity was similar in the four taxa (Table 2). The genetic differentiation among taxa was relatively high (global F_ST = 0.10, D_JOST = 0.23). Magnolia orbiculata contributed mainly to this result since it was highly differentiated from the three other taxa, while M. minor and M. oblongifolia were the less differentiated taxa (Table 3). The Bayesian clustering analysis clearly provided three genetic clusters as the unambiguously best solution in the two analyses (with and without M. oblongifolia) (Fig. 8A, B, Suppl. material 5: fig. S5A, B). In the following, an individual was considered to be correctly assigned to a unique genetic cluster if the ancestry coefficient of this individual to this cluster was higher than or equal to 0.9. One of those clusters corresponded obviously to M. orbiculata (red cluster in Fig. 8C). The 88.8% (32/36) of individuals from M. orbiculata were assigned to this cluster, while the 4/36 M. orbiculata individuals were considered unclear. The second cluster (green cluster in Fig. 8C) consisted mainly of the majority of M. minor (171) and M. oblongifolia (16) individuals, but also included some individuals (14) of T. ophiticola. (Fig. 8C). The third cluster (blue cluster on Fig. 8) was predominantly composed of T. ophiticola with only one individual of M. minor. We will therefore refer hereafter to the “orbiculata”, “minor-oblongifolia” and “ophiticola” genetic clusters, keeping in mind that ancestry coefficients within each taxon of these genetic clusters still varied. Indeed, despite a clear delimitation between three genetic clusters, a significant proportion of individuals (130/461) displayed genetic admixture (on the basis of a 0.9 admixture coefficient value as a threshold). Based on these “admixed” individuals, the level of genetic admixture varied according to taxa. The mean value of probability to belong to their a priori cluster (defined from their taxonomic status) was 0.624 (±0.104) for M. orbiculata, 0.678 (±0.172) for M. minor (including M. oblongifolia), and only 0.477 (±0.262) for T. ophiticola.

Figure 8. 

Structure results of Magnolia subsect. Talauma in Cuba for the complete dataset A Delta K plot B the mean Ln(K) plot C representative bar plot (out of 100 en replicates) for K = 3.

Magnolia orbiculata was strongly homogeneous pertaining to ancestry coefficient values with only four individuals displaying genome admixture with the “minor-oblongifolia” cluster (Fig. 8C). This is strongly in agreement with what was observed for leaf characteristics. Magnolia minor and M. oblongifolia displayed a high level of genome admixture with the “ophiticola” cluster. 35 individuals (16.1%) of the individuals of M. minor showed genome admixture with the “ophiticola” cluster. In M. oblongifolia, 51% of the individuals exhibited an ancestry coefficient over 0.9 to the “minor-oblongifolia” cluster, the rest showed high admixture levels. Moreover, it is noticeable that one individual of M. oblongifolia displayed a very high ancestry to M. orbiculata. The localities of Cupeyal del Norte (CN), Monte Fresco (MF), Piedra la Vela (PV), and Cayo Guam (CG) show the highest levels of misclassification of M. minor and M. oblongifolia into the “ophiticola” cluster.

For T. ophiticola, 56.4% (97/172) of individuals could be assigned to the “ophiticola” genetic cluster while 8.14% (14/172) could be assigned to the “minor-oblongifolia” genetic cluster (referred to as “misclassified” individuals hereafter). Similar to M. minor and M. oblongifolia, many individuals of T. ophiticola (61/172) also displayed signals of genetic admixture, mainly with the “minor-oblongifolia” cluster, but also, for a few of them, with the “orbiculata” cluster. The localities of Subida a la Melba (SM), Mina Iberia (MI), and Sur de las Delicias del Duaba (SDD) showed the lowest levels of misclassification. Four individuals from Cupeyal del Norte (CN) were clustered with the group of M. orbiculata. Most individuals from La Melba (MEL), Pico Cristal (PC), and Monte Fresco (MF) showed an ancestry coefficient similar to the “minor-oblongifolia cluster”. The clustering analysis without individuals of M. oblongifolia also provided K = 3 as the best solution (Suppl. material 5: fig. S5A, B). Moreover, it was striking that this analysis provided an ancestry pattern very similar (Suppl. material 5: fig. S5C) to the analysis including this taxon (Fig. 8C). This demonstrated the very good stability of inferences on individuals’ ancestry coefficients which could be explained by the strong genetic delimitation between the three identified genetic clusters.

The PCA analysis on the whole SSR data set showed that the 200 first principal components explained 99.3% of the variation, which were therefore kept for the discriminant analyses. Based on the number of taxa that have been defined across the history of Cuban Talauma taxonomy, but also on the STRUCTURE results, two solutions for the number of genetic clusters were considered in the following discriminant analysis (DAPC) K = 3 and K = 4. When K = 3, individual assignment displayed a pattern very similar to that found with the Bayesian clustering approach; with one cluster predominantly composed by M. minor and M. oblongifolia, the other cluster with T. ophiticola, and the third one with the individuals of M. orbiculata. In the three clusters, some level of misclassification was found. Many individuals “misclassified” in the DAPC analysis were the same that were “misclassified” based on the STRUCTURE analysis. The DAPC analysis confirmed the correspondence of M. orbiculata to a unique genetic cluster as expected because of its high genetic differentiation from the three other taxa (Suppl. material 6: fig. S6A). For K = 3 only one individual of T. ophiticola showed an assignment probability value less than 0.9.

K = 4 (Suppl. material 6: fig. S6B), seems to be a less meaningful solution. In this case, three clusters were predominantly composed of M. minor, T. ophiticola and M. orbiculata respectively, confirming the main pattern found with K = 3, with the difference that a higher proportion of M. minor and T. ophiticola, but also a majority of M. oblongifolia were not assigned to unique clusters. When K = 4 seven and three individuals of M. minor and M. oblongifolia, respectively, showed probabilities values under 0.9. As for structure, the analysis without considering M. oblongifolia with K = 3 displayed very similar results to the analysis including this taxon (Suppl. material 6: fig. S6C); in this case, only one individual of T. ophiticola showed probabilities values under 0.9.

Overall, the morphological and genetic classifications were highly congruent (χ² = 173.69, p < 0.0001). The concordance between the two classifications (genetic and morphology) was especially high for Magnolia orbiculata and M. minor, and to a lesser extent for M. oblongifolia and Talauma ophiticola (Fig. 9). In this last taxon, the classification of several individuals based on genetic markers on one side and leaf traits on the other side were not congruent. Only a few genetic and morphological inconsistencies were also observed in M. minor and M. oblongifolia.

Figure 9. 

Heatmap with the congruence between morphological (MT: Multivariate, OUT: Elliptic Fourier Descriptors, LM: Matrix of Landmark) and genetic (Structure) cluster probabilities, inside each taxon of Magnolia subsect. Talauma in Cuba (A–C) Magnolia minor (D–F) M. oblongifolia (G–I) Talauma ophiticola (J–L) M. orbiculata. The blue color represents the number of individuals (less individuals: light blue; more individuals: dark blue).

The observed leaf morphological variability for Cuban magnolias was higher than that described by previous studies. According to the two taxa CS, the values of leaf length and width were higher than those reported by León and Alain (1950, 1951) for Magnolia minor and M. orbiculata. Likewise, in the four taxa CS, these values were higher than what was previously reported by Bisse (1974, 1988), except for M. oblongifolia. This difference with previous studies is due to the larger sample size used in the present work and its wider geographic representativeness. In the three taxa CS, Palmarola et al. (2016) reported similar values of length and width for M. minor and lower values for M. oblongifolia and M. orbiculata. The average values of Bisse Index (B_CI) were similar to those reported by Bisse (1974, 1988) for M. orbiculata and M. oblongifolia. For M. minor and T. ophiticola, the average values of B_CI are slightly lower and slightly higher, respectively, than those reported by Bisse (1974, 1988). The high level of morphological differentiation between taxa observed in this study reinforces the value of leaf characteristics in taxonomic studies of Cuban magnolias (León and Alain 1950, 1951; Alain 1969; Bisse 1974, 1988; Imkhanitzkaja 1991, 1993; Hernández-Rodríguez 2014; Palmarola et al. 2016). Leaf morphological data are key traits for species delimitation (Jensen et al. 2002; Jensen 2003). This study confirmed they are highly relevant in groups like Magnolia, where very little variation is observed in flower and fruit characters (Treseder 1978).

In our study, Magnolia orbiculata was clearly distinguished from the other taxa of Magnolia subsect. Talauma in Cuba based both on morphology and genetic markers. The previously observed large variation of leaf morphology across subsection Talauma in Cuba, although based on the observation of only a few specimens, has been the basis for several authors to consider a unique species in this subsection, therefore including M. orbiculata within M. minor (Howard 1948; Alain 1969; Borhidi and Muñiz 1971; Lozano-Contreras 1994). In contrast, the present study, as well as lines of evidence already brought by molecular phylogeny of the subsection Talauma (Veltjen et al. 2022) and by studies on the ecological niches of Cuban Talauma (Testé et al. in press), strongly supported that Magnolia orbiculata should be considered as a well-delineated species.

However, in our study, a few cases of confusion with M. minor (sensu Bisse 1988 and Palmarola et al. 2016) on the basis of leaf morphology traits were observed. This confusion may be explained by the similar rounded shape and relation width-length present in both taxa. Different specialists have erroneously identified some herbarium specimens of M. orbiculata as M. minor in the past (personal observation in herbarium records). Moreover, our data showed that very few M. orbiculata individuals displayed genetic admixture with M. minor. Similarly, a few M. oblongifolia and T. ophiticola individuals displayed genome admixture with M. orbiculata. The levels of genetic differentiation among species are influenced by the time of separation and the amount of gene exchange (Hey and Pinho 2012). Genetic variation shared between closely related species may be due to the retention of ancestral polymorphisms because of incomplete lineage sorting (ILS) and/or introgression following secondary contact (Zhou et al. 2017).

Distinguishing between those two causes from observed patterns is challenging, although coalescence modeling can help (e.g. Zhou et al. 2017; Meleshko et al. 2021). However, in the case of M. orbiculata relative to other taxa, regular gene flow seems to be unlikely. The very clear morphological and genetic differentiation of M. orbiculata with other taxa in Cuba strongly suggested that the lowland between the Sierra Maestra (habitat of Magnolia orbiculata) and Nipe-Sagua-Baracoa (habitat of the other species) may have acted and still acts as a barrier to gene flow by strongly limiting pollination and seed dispersal. Hernández-Rodríguez (2022) reported high levels of genetic differentiation between Magnolia cubensis Urb. subsp. cubensis (from the Sierra Maestra) and Magnolia cristalensis Bisse (from Nipe-Sagua-Baracoa), both from subsection Cubenses. Vázquez-García et al. (2016) stated that allopatric speciation seems to be a major driver of Magnolia diversification in the Neotropics. Therefore, it seems more likely that the admixture signal between M. orbiculata and the other taxa could rather be explained by shared ancestral polymorphism with other Cuban talaumas due to the likely recent diversification of the subsection in Cuba, that is less than 5 mya according to Veltjen et al. (2022), and the recent separation of M. orbiculata from the other taxa. However, the possibility of rare events of inter-taxa hybridization involving M. orbiculata as one parent cannot be totally ruled out, especially because individuals displaying admixed genome involving M. orbiculata have intermediate ancestry coefficients, which is compatible with a hypothetical first- or early-generation hybrid status. Testé et al. (in press) have also shown that the ecological niche of M. orbiculata is differentiated from that of the other taxa considered in this study. This may suggest that selection against first- or early-generation hybrids due to local adaptation could also contribute to preventing genetic exchanges between that taxon and the other taxa of Magnolia subsect. Talauma in Cuba.

Undoubtedly, our data confirmed that the main taxonomic issues concern the north-eastern Cuban populations of Magnolia subsect. Talauma. León and Alain (1950) have stated that individuals of M. minor with more oblong leaves, considered by them as Talauma minor var. oblongifolia, may belong to a different species. However, the authors did not assign the species rank to this group because of the absence of reproductive structures in the available specimens. On the other hand, Bisse (1974, 1988) proposed to divide Magnolia minor (sensu León and Alain 1950, 1951) into three separate species (M. minor, M. oblongifolia and T. ophiticola). Our morphological and genetic data did not support those two proposals. Indeed, concerning M. oblongifolia (sensu Bisse 1988), the foliar phenotype observed in this taxon appears to be intermediate between M. minor and T. ophiticola. A recent diversification process or natural hybridization might explain the intermediate characteristics of M. oblongifolia, as has been observed for Quercus species (Burgarella et al. 2009; An et al. 2017) and the genus Rhizophora (Francisco et al. 2018). Rather, considering M. minor and M. oblongifolia as separate taxa is supported neither by morphological (see Figs 4–7) data nor by genetic data (Fig. 8) of the present study. On the other hand, the existence of a single species, including those three taxa, (Magnolia minor sensu León and Alain 1951) was supported neither by our morphological results, nor by genetic markers, which both showed a clear differentiation between M. minor and T. ophiticola. Yet, our results did not support either the combination of Talauma ophiticola and Magnolia oblongifolia (sensu Bisse 1974, 1988) in a unique taxon, as recently proposed by Palmarola et al. (2016) on the basis of the specimen HFC 5358 from the coast of Moa, which shows both oblong and elliptical leaves. Nevertheless, the delimitation of T. ophiticola is still challenging. In the present study, a significant proportion of individuals that were assigned to this taxon based on leaf morphology was unambiguously assigned to the “minor-oblongifolia” genetic cluster, while only one individual of M. minor was assigned to the “ophiticola” genetic cluster. Also, for each taxon, a significant proportion displayed high genetic admixture between the two genetic clusters identified (and as discussed above rare cases of admixture with M. orbiculata). This could be explained by a recent diversification of the three taxa that led to numerous genetic loci with incomplete lineage sorting and to overlaps in the distribution of morphological traits. In trees, factors such as long generation time, and large effective population sizes, increase the opportunity of sharing ancestral polymorphisms through incomplete lineage sorting which makes species identification based on neutral markers even more problematic (Zhou et al. 2017).

The taxa from the north-eastern part of Cuba live in the same habitats and in similar ecological conditions (Testé et al. in press), a situation that is not favorable for the emergence of reproductive barriers. Moreover, those three taxa are also found in sympatry in several locations. The phylogenetic closeness between those three taxa has recently been reported by Veltjen et al. (2022). Therefore, the high admixture level observed in these taxa with SSR markers, as well as the few cases of reciprocal “miss-assignment“, suggest gene flow between the taxa of northeastern Cuba has occurred recently and may still be occurring, producing recombinant and therefore intermediate genotypes and phenotypes. This hypothesis is reinforced by the observation that reciprocal genetic admixture between the two genetic clusters, corresponding mainly to M. minor and T. ophiticola, is more frequent in the localities where both taxa occur. According to Callaway (1994), hybridization is a common process in magnolias and is more common when the distribution ranges of two or more highly related taxa overlap (Soltis and Soltis 2009).