1) Genetic diversity
The higher genetic diversity of Magnolia cubensis subsp. cubensis compared to M. cubensis subsp. acunae (Table 2) could be explained by the larger geographic extent of the population (Fig. 1b), the presence of more habitat heterogeneity (Capote and Berazaín 1984; Simón 2020), the habitat being conserved better (Palmarola et al. 2012; González-Torres et al. 2016), or a combination of the three. Gitzendanner and Soltis (2000) state that many fragmented habitats have lost the ability to support plant populations large enough to maintain a mutation-drift equilibrium. Therefore, the fragmentation of the habitat of M. cubensis subsp. acunae increases the risk that the population will reach a threshold below which it is not viable. Fortunately, its longevity guarantees that the time between generations is extended and therefore, the loss of alleles due to genetic drift is moderated (Honnay and Bossuyt 2005).
The mean values of allelic richness, observed and expected heterozygosity (Table 2) found for M. cubensis subsp. acunae exceeded that reported by Hernández et al. (2020a) (Topes de Collantes: AR = 4.10 Ho = 0.50 He = 0.54; Lomas de Banao: AR = 2.93 Ho = 0.36 He = 0.41), most likely because of the higher number of microsatellite markers used in the present study. The diversity measures were lower than the diversity values reported by Veltjen et al. (2019) in M. cubensis subsp. cubensis in Pico Turquino (Ho = 0.60 He = 0.61), but similar to the values of M. cubensis subsp. acunae in Topes de Collantes (Ho = 0.60 He = 0.59). These authors used a higher number of markers (30 loci), even when the number of samples used per subspecies was less (20 individuals). A similar behavior was reported by Hall and Beissinger (2014), who compared diversity values in several studies based on sample size and the number of markers used. Both results show that in population genetic studies, the cost-benefit relationship should guarantee an equilibrium between the number of markers and the number of samples.
When compared with other similar studies in non-Caribbean Magnolia, in general, the average values of the number of alleles per locus, allelic richness, observed and expected heterozygosity are lower to those reported for Magnolia tomentosa Thunb. (Setsuko et al. 2004), M. stellata (Siebold & Zucc.) Maxim. (Setsuko et al. 2007; Tamaki et al. 2008), M. shiluensis (Chun & Y.F.Wu) Figlar (as Michelia shiluensis Chun & Y.F.Wu; Deng et al. 2020) and five Magnolia species of Veracruz, Mexico (Aldaba-Núñez et al. 2021); but higher than in M. sieboldii K. Koch (Kikuchi and Isagi 2002), M. coriacea (Hung T.Chang & B.L.Chen) Figlar (as Michelia coriacea Hung T.Chang & B.L.Chen; Zhao et al. 2012) and M. wufengensis L.Y.Ma & L.R.Wang (Wang et al. 2017).
When we compare the genetic diversity of the localities and their SNAP management category – which serves as a proxy for the degree of habitat disturbance, we do not retrieve that the habitats categorized to be in a higher degree of disturbance harbor less genetic diversity. In Magnolia cubensis subsp. acunae, the Topes de Collantes locality was the most diverse, while Lomas de Banao showed lower values for all the measures. However, Hernández et al. (2020a) found that, the Ecological Reserve Lomas de Banao has more conserved forests and less fragmentation than the Protected Natural Landscape Topes de Collantes, in correspondence with its management category (SNAP 2009). In M. cubensis subsp. cubensis, the most diverse locality was the National Park La Bayamesa and the least diverse were the Ecological Reserves Pico Caracas and El Gigante; yet all three have a similar degree of management category (SNAP 2009).
Our genetic diversity results may be determined by the number of individuals and the sample size of each locality. Poudel et al. (2014) demonstrated that the sample size can have a significant effect on the estimate of expected heterozygosity in microsatellite markers. Furthermore, Vergeer et al. (2003) state that the population size is positively correlated with the heterozygosity observed. However, most likely, the genetic diversity results are explained by demographic causes. On the one hand, Cuban magnolias were highly exploited for timber for many years (Carreras and Dechamps 1995; Álvarez 2003). On the other hand, low natural regeneration has been reported by Palmarola et al. (2018) in M. cubensis subsp. acunae, and by Molina-Pelegrín et al. (2014) and Testé et al. (2019) for M. cubensis subsp. cubensis in El Gigante and Gran Piedra, respectively. The four less genetically diverse localities (i.e. Lomas de Banao, El Gigante, Pico Caracas and Gran Piedra) hence may partly be explained by their demographic history.
2) Genetic structure
The PCoA analysis (Fig. 2a) and the structure analysis (Fig. 3a) showed two well-defined groups, corresponding to each subspecies. Currently, these taxa are considered subspecies (Imchanitzkaja 1991, 1993; Figlar and Nooteboom 2012; Greuter and Rankin 2017). However, the large distance between their populations (ca. 400 km), morphological differentiation (Hernández 2014), ecological features (Simón 2020), phylogenetic evidence (Veltjen et al. in review), as well as our results, support the decision to consider them as different species.
In the analyses per subspecies, all spatial genetic patterns (genetic structure: Fig. 2b, c and Fig. 3b, c, isolation by distance and spatial autocorrelation: Fig. 4) show similar results. In Magnolia cubensis subsp. acunae there is no substructure, no isolation by distance and no spatial autocorrelation. In M. cubensis subsp. cubensis, there is substructure, with significant isolation by distance and significant values in the first five distance classes of the spatial autocorrelation analysis. The structure (Fig. 3c) and PCoA (Fig. 2c) analyses clearly separate the GRP (Gran Piedra) locality from the other three localities. Next, the REC (Pico Caracas) locality is genetically detectable as a randomly mating population; however, there are migrants of this genetic cluster detectable in the TUR (Pico Turquino) populations. In the PCoA, REC is differentiating a little bit from the rest of the localities, but still overlaps.
3) Functional connectivity
The fixation index (FST) is a measure of the genetic divergence between subpopulations. Values close to zero, such as those found in Magnolia cubensis subsp. acunae (Table 3), indicate that there is similar heterozygosis in each locality compared to the value of the two localities together, so genetic differentiation is low (Hartl and Clark 1997). High values of this coefficient (FST > 0.15), such as those detected in M. cubensis subsp. cubensis (Table 3), show there has been less gene flow. However, between pairs of geographically contiguous localities such as Pico Caracas – Pico Turquino, Pico Turquino – La Bayamesa and La Bayamesa – El Gigante (Table 3), the genetic differentiation was moderate. If we compare our results with the FST values for Caribbean Magnolia (Veltjen et al. 2019), the FST values between some localities of M. cubensis subsp. cubensis (Pico Caracas – El Gigante: FST = 0.23, Gran Piedra and the rest: 0.20 < FST < 0.30), are comparable to those found for sister species pairs more than for population pairs. In general, the genetic differentiation between localities supports the results of the genetic patterns. The locality of Pico Caracas presents values of genetic differentiation with Pico Turquino, similar to (or lower than) other pairs of localities that belong to the same genetic group (e.g., La Bayamesa – El Gigante). The results of genetic differentiation, together with those obtained in the Structure analysis and the PCoA, suggest to consider Gran Piedra as a separate genetic entity and not to appraise Pico Caracas as an independent genetic unit.
The clear genetic cluster in the Structure barplots for the Gran Piedra individuals (Fig. 3c), as well as their differentiation with the rest of the localities (Table 3), can be the cause of the large number of private allelic variants at this locality (Table 2), resulting from the geographical distance and isolation (Fig. 1b) or the diversification history of the Cuban magnolias. The geographical distance between localities can regulate the gene flow between them and has a direct effect on the existing genetic differentiation. The characteristic red-orange arilloid of the mature Magnoliaceae seeds (Judd et al. 2016) is easily visible and ingested by birds. As we do not find an indication of (recent) gene flow between Gran Piedra and the other localities, it could be that the distance of 88.49 km is either too far to overcome by the seed dispersers, which causes the barrier in gene flow; or it may be because of the landscape limiting the birds, even though they could overcome this distance.
In the study by Veltjen et al. (in review), the precise relationship among the two Magnolia cubensis subspecies and a third Cuban species from that subsection (M. cristalensis Bisse) was unclear, and the estimated divergence time of that clade was 5–2 million years ago. Hence, considering the simulations of Iturralde-Vinent (2006) and the work of Borhidi (1996), it is a possibility that Gran Piedra was the first Cuban locality to which Magnolia members of the subsection Cubenses migrated after which the rest of the Cuban localities were reached and diversification occurred. Our study shows clear genetic differences between the individuals of M. cubensis subsp. cubensis of Gran Piedra and those of the rest of the Sierra Maestra mountain massif. This was not unexpected based on its complex taxonomic history. Magnoliacubensis subsp. cubensis is one of the Cuban magnolias that underwent numerous taxonomic changes, mainly at the infraspecific level. The current subspecies was treated by Imchanitzkaja (1991, 1993) as two different taxa: Magnoliacubensis subsp. cubensis with distribution in the Gran Piedra mountain range and Magnolia cubensis subsp. turquinensis Imkhan. from the Sierra Maestra mountain range. Recently, Greuter and Rankin (2017) synonymized M. cubensis subsp. turquinensis within M. cubensis subsp. cubensis based on the similarity between taxa in the leaf shape and leaf areole size.
Considering the different data (Fig. 2, 3, Table 2, 3), and according to the criteria of Ryder (1986) and Moritz (2002), we can define one evolutionary significant unit in Magnolia cubensis subsp. acunae, and two significant evolutionary units in M. cubensis subsp. cubensis: Pico Caracas – Turquino – La Bayamesa – El Gigante and Gran Piedra. Defining the conservation units within each subspecies makes it possible to preserve their evolutionary processes, previously not considered in conservation (Moritz 2002), and it is essential for the long-term survival of these threatened magnolias.
4) Structural connectivity
The coincidence between genetic and geographical barriers shown by Monmonier's algorithm, should not be taken as evidence of a relationship between gene flow patterns and the studied landscape features. This can be explained in a relatively simple way in Magnolia cubensis subsp. cubensis: if the landscape characteristics reflected the gene flow patterns, only large genetic differences would be found between Gran Piedra and the rest of the localities, due to the existence of Santiago de Cuba basin (CNNG 2000), which represents an area of differences in altitude, topography, vegetation and land use, between both mountain ranges (Sierra Maestra and Gran Piedra) (Fernández and Pérez 2008). However, the genetic fixation between Pico Caracas – El Gigante (both in the Sierra Maestra mountain range) was greater than between El Gigante – Gran Piedra (localized in different mountain ranges: the firth in Sierra Maestra and the second one in Gran Piedra).
From the results of the partial Mantel tests, it can be inferred that the geographical, rather than the environmental, distances explain the genetic differences found between localities. However, statistically significant correlations were detected between geographic and environmental distances for the two taxa, by controlling for genetic effect. Environmental differences between the distribution areas of Cuban magnolias are reported by Simón (2020). The latter study found that the greatest niche similarity in the family did not correspond to the taxa that were phylogenetically closest, but instead to those that live in the same mountain range. However, this analysis was based on the comparison of climatic niche between species and did not include the possible microclimatic differences in the whole distribution range for each taxon.
The correlation between geographical and environmental distances can be explained by microclimatic differences. In general, the mountain systems where these magnolias live have higher precipitations and lower temperatures than those of the rest of the country (Durán 2012). In turn, a single mountain at a higher elevation corresponds to an increase in the average precipitation and a decrease in the average temperature; therefore, individuals found at different elevations will also develop under different microclimatic conditions (Díaz 1988).
The differences in the relation of the geographical distances, the cost distances and the fixation index (FST) between localities for each subspecies, not necessarily indicate that they have different dispersal ability or seed dispersers. Their seeds are very similar (Judd et al. 2016), the possible dispersers in Cuba (Hernández et al. 2020a) are distributed over the whole island and are permanent residents (Garrido and Kirkconnell 2011). Hence, it is most likely that these differences are because the Least-cost-path models reflect the present environmental data, and not the environmental data from when the dispersal happened.
Although the Least-cost-path models may be related to a resistance to dispersal in each mountain massifs, it is necessary to clarify that the models were based on the habitat requirements of Magnoliacubensis. It is possible that the seed dispersers have similar requirements, in which case the model would represent the resistance to the movement of the disperser, and therefore, the gene flow of these species. Although the method has some limitations, the Least-cost-path models can be used for the establishment of biological corridors (Li et al. 2010; Alexander et al. 2016).
It should be noted that the analyses presented above (partial Mantel tests, Least-cost-path), attempted to evaluate the effect of landscape characteristics and some environmental variables on migration, dispersal and gene flow of Magnoliacubensis, not the possible interaction between the environment and the genetic variation of these taxa. Although some studies of landscape genetics in plants (Vandepitte et al. 2007; Schmidt et al. 2008) have focused on evaluating whether genetic diversity is related to soil type, elevation, humidity, temperature or vegetation, the type of marker used in the present study excludes this type of analysis. Neutral genetic diversity (such as that quantified by microsatellite markers) is not directly related to adaptive genetic variation and is only indirectly influenced by local environmental factors if the latter affect processes such as gene flow (Ouburg et al. 2010). To evaluate the direct effect of environmental factors on the genetic diversity of these magnolias, adaptive genetic diversity would have to be quantified through another type of molecular marker such as Single Nucleotide Polymorphism (SNPs). They make it possible to test which SNPs are in genes or close to genes under selection and then to correlate this to environmental variables.