Mimosa catharinensis is a plant species endemic to the Brazilian Atlantic Rainforest with an extremely restricted distribution, limited to one population located in the Rio Vermelho State Park (Parque Estadual do Rio Vermelho – PAERVE; Figure 1), Santa Catarina State, Southern Brazil (Burkart, 1979; Ferreira, 2010). The only documented population of M. catharinensis (~700m²) occurs in one of the most heavily impacted areas of the PAERVE (Figure 1) due to deforestation and environmental degradation for over 200years (Berenhauser, n.d.). This narrow endemic and critically endangered plant species (CONSEMA, 2014) is a lianescent shrub that produces masculine and hermaphrodite white, glabrous, and tetramerous flowers (Burkart, 1979). Pollen dispersal occurs via zoophily (e.g., Apis mellifera; Silva et al., 2005), although nectar-dependent visitors are infrequent due to a lack of floral nectar (Silva et al., 2005). The pods are of the indehiscent craspedium type and are linear-oblong with curved prickles on the margins (Medeiros and Stefani, 2018). Although the reproductive structures of M. catharinensis have been characterized morphologically, pods without seeds were observed during field collection in 2019. Based on herbarium records, the lack of seed production in M. catharinensis was first recorded in 1994 (Voucher 30,482; FLOR herbarium, UFSC). As no seedlings were found in its area of occurrence, M. catharinensis may be considered an ecologically sterile plant species (i.e., a plant species in which the recruitment rate over time is nil as a consequence of lack of sexual reproduction). As successful in vitro pollen germination had been reported for M. catharinensis (Silva et al., 2005), further pollination studies are needed to improve our understanding of its reproductive biology.

Although there are efforts by the Santa Catarina State government (Gasper et al., 2012) to generate genetic and ecological data to inform conservation efforts and to safeguard plant species threatened with extinction, no previous evolutionary studies – including phylogenetic analysis and genomic characterization – have been conducted on this critically endangered species.

Since 2002, the PAERVE has been recognized by UNESCO as one of the core areas of the Atlantic Forest Biosphere Reserve and consists of 1,530 hectares of Atlantic Rainforest biome (Ferreira, 2010; Bechara et al., 2013). The introduction of pine and eucalyptus into the park has brought about significant changes to its restinga (sandy plains) vegetation and is the main cause of degradation of the coastal ecosystem of Santa Catarina Island (Ferreira, 2010). Currently, it is estimated that 750ha of the park are covered by Pinus spp., of which about 33% are located in dunes and sandbanks that have been progressively invaded since pine was first introduced (Bechara et al., 2013). Nevertheless, approximately 400ha of dense rainforest and 250ha of restinga vegetation fragments remain protected in the park (Bechara et al., 2013), the latter located mainly at the southern and northern tips of the conservation area (Ferreira, 2010).

We collected leaf samples from all 33 identified adult individuals of unknown age that constitute the only remaining M. catharinensis population. Leaves were stored in silica gel for later extraction. One voucher (FLOR30482) is deposited at the FLOR herbarium at the Federal University of Santa Catarina – UFSC.

DNA was extracted from leaf samples of all collected individuals employing the NucleoSpin® kit (Machereney-Nagel GmbH & Co. KG) following manufacture’s guidelines. After extraction, the quality of each sample was verified using Thermo Scientific NanoDrop 2000 Spectrophotometer (Thermo Fisher Scientific Inc.) and the concentration of double-strand DNA (dsDNA) was obtained by using Qubit dsDNA Assay Kit (Invitrogen). The genomic library was prepared using a double-digest restriction site-associated DNA sequencing (i.e., ddRADseq) protocol (Peterson et al., 2012) with modifications proposed by Nazareno et al. (2017). Digestions reactions were performed using 0.5μg of genomic DNA, 5,000 units of MseI, 5,000 units of HF-EcoRI, and 1x CutSmart buffer (New England Biolabs) in a total 20μl reaction volume for 3h at 37°C, ending with a 20min deactivation step at 65°C (Nazareno et al., 2017). Reactions products were cleaned using Agencourt AMPure XP system (Beckman Coulter) according to manufacturer’s recommendations. The dsDNA was quantified using Qubit and the amount of DNA was standardized before to procedure with the ligation. We used 80ng DNA, 1M of MseI adaptor, 0.33M of EcoRI double-strand adaptor unique for each sample, 1U of T4 DNA ligase (New England Biolabs) and 1.40μl of T4 ligase buffer in 30μl ligation reaction that was stored at a room temperature (23°C) for 30min After, heat-killed the reaction at 65°C for 10min and then, slowly cooled until reach room temperature. Details of barcode sequences can be found in Nazareno et al. (2017). Products of the ligation reaction were purified using the Agencourt AMPure XP system and used in PCR reactions. PCR reactions were carried out in a total volume of 20μl, containing 13.5μl of ligation product, 0.2M dNTPs, 1.0μM MgCl₂ 0.5U iProof™ High-Fidelity DNA Polymerase (Bio-Rad) and 2X of iProof buffer. The PCR protocol (98°C for 30s, 20cycles of 98°C for 20s, 60°C for 30s and 72°C for 40s, followed by a final extension at 72°C for 10min) was performed in an Eppendorf System. The amplicons were cleaned with the Agencourt AMPure XP system and quantified using Qubit dsDNA assay Kit. We used an automated size-selection technology (i.e., Pippin Prep; Sage Science, Beverly, MA, United States) at 2% agarose cartridge to select DNA fragments at a target range size of 375–475bp. The library was sequenced (100-bp single-end) on a lane of Illumina HiSeq 2,500 flow cell (Illumina Inc., San Diego, CA, United States) at The Centre for Applied Genomics in Toronto, Canada.

Data quality was checked using the program FastQC version: 0.11.8.¹ The file containing raw sequence reads was analyzed in Stacks 2.41 (Catchen et al., 2011, 2013; Rochette et al., 2019) using de novo assembly. Initially, we used the process_radtags program in Stacks to examine individual barcodes, enzyme cutsite integrity and to demultiplex the data. For barcodes rescue we admitted at most two mismatches (−barcode_dist 1 2). We filtered raw reads using the phred score 33 and used the option −t to trim all of it at 85 base pairs (−phred 33, −t85). With the process_radtags output, we run the Ustacks program, that uses a maximum likelihood framework to aligns short-read sequences and create putative alleles (stacks). The maximum distance allowed between “stacks” was two nucleotides. The minimum depth of coverage required to create a “stack” and the maximum number of stacks at a single de novo locus was set as three. We enabled the deleveraging algorithm (−d) to solve merged tags, and a bounded-error model to identify SNPs. An alpha value of 0.05 and an upper bound of 0.1 were used. In the next step, we performed Cstacks to merge alleles of all samples and create a catalog of consensus loci. Three mismatches were allowed between loci in catalog build. Each individual sample loci were then compared against the catalog through the Sstacks program. We ran tsv2bam program to organize the single-end reads by locus, instead by sample, creating a BAM file that was used as input to the Gstacks program. Gstacks uses the single-reads to build contigs and merges them into loci. It also aligns reads to the loci using Ukkonen’s algorithm, identifies SNPs for each locus and each individual genotype and converts SNPs into phased haplotypes. Finally, we used POPULATIONS in Stacks (Catchen et al., 2011, 2013; Rochette et al., 2019) to generate and export the SNP data set, as well as FASTA files containing the per-locus consensus sequences, and individual loci sequences, that were applied in the further analyses. To assess the effects of missing data on genetic estimates, we run POPULATIONS several times with a Minor Allele Frequency (MAF) of 5% (−min_maf 0.05), a maximum heterozygosity (−max_obs_het) of 0.65, and the percentage of missing data varying from 0 to 30%. All data sets include one random SNP per locus. To avoid potential bias sources in the forward genetic analyses, Hardy–Weinberg (H–W) equilibrium tests was done using the adegenet and pegas packages² (Jombart, 2008; Paradis, 2010; Jombart and Ahmed, 2011) implemented in R. In addition, linkage disequilibrium (LD) between loci was tested using Arlequin 3.5.2 (Excoffier and Lischer, 2010). Type I error rates for these tests were corrected for multiple comparisons using the sequential Bonferroni procedure (Rice, 1989), and SNPs that failed the H–W equilibrium test and/or SNP pairs in LD were excluded from further analyses. To exclude non-nuclear loci, the per-locus consensus sequences were aligned against reference chloroplast and mitochondrial genomes using Acacia dealbata (NCBI accession number KX852435) and Acacia ligulata (NCBI accession number MH933866), respectively. We used the BLASTn program³ to identify loci that presented identity greater than or equal to 80%.

To assess the influence of missing data on genetic diversity estimates, we used the “BasicStats” function in the DiveRsity package in the R software environment⁴ (Keenan et al., 2013). We estimated unbiased expected genetic diversity (uH _E; Nei and Roychoudhury, 1974), observed heterozygosity (H _O), and the inbreeding coefficient (Wright’s Fixation Index F _IS; Wright, 1943). Population genetics statistics were averaged across loci using the DiveRsity package in R (Keenan et al., 2013). Confidence intervals at 95% were obtained to evaluate differences among means of genetic parameters estimated for all data sets. We used the final data set to calculate minor allele frequencies using the adegenet package in R (Jombart, 2008; Jombart and Ahmed, 2011).

We used the kinship coefficient (Loiselle et al., 1995) implemented in the SPAGeDi program (version 1.5; Hardy and Vekemans, 2002) to estimate random outcrossing rates (1-s, where s indicates the estimated selfing rate) based on standardized identity disequilibrium for all data sets. Significance for the identity disequilibrium coefficient was obtained with 1,000 permutations, and a jackknife over loci approach was applied to calculate the standard error of outcrossing estimates.

We evaluated the power of discrimination (PD) for the complete single nucleotide polymorphism (SNP) data set (i.e., the data set without missing data; MD) using the equation 1-PI, where PI represents the probability that two individuals drawn at random from a population will have the same genotype at multiple loci (Waits et al., 2001). This analysis was performed in GeneAlEx 6.5 program (Peakall and Smouse, 2012). Taking into account that clonal reproduction has been reported for plant species of the Fabaceae family (e.g., Dev et al., 2010; Gui et al., 2013; Gaddis et al., 2014; Mori et al., 2018; Amor et al., 2020), including M. catharinensis (Rogalski et al., 2005), a clonality test was applied for M. catharinensis using the “poppr” package (Kamvar et al., 2014, 2015) in R version 3.5.3 (R Core Team, 2019). The “bitwise.dist” function was used to calculate a pairwise genetic distance matrix between individuals. Then, we ran the “mlg.filter” function to apply a threshold that defines the minimum distance to determine distinct multi-locus genotypes (MLGs; see Kamvar et al., 2015). Applying a threshold when using SNP markers is important in order to deal with the limited observable genetic differences between genets due to sporadic somatic mutations or issues associated with high-throughput sequencing technologies (e.g., genotyping errors and missing data; Kamvar et al., 2015; Mastretta-Yanes et al., 2015). We chose the average neighbor as the clustering algorithm with a threshold of 0.04, which has been shown to be sufficient to explain intra-clonal variation when SNPs derived from ddRADseq are used (Amor et al., 2020). In addition, a less conservative threshold of 0.06 was used to infer the numbers of MLGs in M. catharinensis. The threshold of 0.06 was chosen in order to minimize the limitation of our experimental design, which did not include replicates to estimate genotyping errors. To visualize the number of putative lineages, we ran the “upgma” function (Phangorn package; Schliep, 2011) to construct ultrametric trees.

The demographic history of M. catharinensis was inferred using a composite likelihood-based approach implemented in FastSimCoal v. 2.6 (Excoffier et al., 2013). VCF files were converted into the *.arp format using the program PGDSpider 2.1.1.5 (Lischer and Excoffier, 2012). Assuming genetic distances of 0.04 and 0.06, we obtained two observed folded Site Frequency Spectrums (SFS based on minor allele frequency) for our filtered SNP data set using the Arlequin 3.5 program (Excoffier and Lischer, 2010). We tested four evolutionary scenarios: constant population size, bottleneck, population decline, and population expansion (Figure 2). We performed 100 replicates for each tested model and folded SFS. The parameters for each run were estimated based on 100,000 simulations and 40 ECM optimization cycles. We used an overall substitution rate of 7×10⁻⁹ per site/generation as reported for Arabdopsis thaliana (Ossowski et al., 2010). The maximum estimated likelihood for each demographic scenario and SFS were used to identify the model with the best fit, which was chosen based on Akaike’s Information Criterion (AIC; Akaike, 1973, 1974) estimated as ΔAIC. All models were ranked, with the model with the lowest ΔAIC being considered the most plausible. It is important to note that models with ΔAIC ≤2 have substantial support (Burnham and Anderson, 2001). We also assessed the probability that each model is the best fit by estimating Akaike weights (w _i). Once the best model was selected, its estimated parameters were used to simulate 100 SFS data sets to built confidence intervals based on bootstrap distribution.

The Illumina sequencing generated a total of 79,759,227 raw reads, of which less than 1% was discarded due to low quality. The average (±SE), minimum, and maximum retained reads were 2,384,161 (±115,118 SE), 1,470,849, and 4,011,918. The total number of genotyped RADtag loci was 274,113 with a mean coverage depth of 12.7 (±2.7 SD). After Bonferroni adjustment, significant deviations from Hardy–Weinberg equilibrium were observed for all data sets (Table 1). No linkage disequilibrium was observed between loci after Bonferroni correction for k tests (Table 1). After alignment, one locus matched the plastome of A. dealbata (96.47% identity) and it was excluded from the data set. The final numbers of unlinked SNPs used in each data set are shown in Table 1.

Genetic diversity parameters (uH_E and H_O ) and Wright’s Fixation Index (F _IS) for data sets with different percentages of missing data are shown in Table 1. For all genetic indices (uH_E , H_O , and F _IS) obtained for the M. catharinensis population, no statistical differences were observed among the data sets (Table 1). Further analyses were performed using the data set without missing data, which presented a minor allele frequency (MAF) averaged at 0.125 (±0.093 SD).

Based on mating system analysis for the data sets with different percentages of missing data (Table 1), the outcrossing rates varied from 0.883 (SE±0.0483) to 0.909 (SE±0.011). Although the outcrossing rate was inflated with an increase in MD, no statistical differences were observed among the data sets (Supplementary Table S1).

The results of the PD analysis showed that at least 13 loci are enough to discriminate closely related individuals with high accuracy, indicating that the data set with 0% missing data (i.e., 128 SNPs) has high discriminatory power to identify M. catharinensis individuals. The defined threshold of 0.04 did not cluster MLGs into multilocus lineages (MLLs; Figure 3A). While the genetic distance of 0.04 did not indicate the presence of clonal genotypes for M. catharinensis, the threshold of 0.06 indicates 23 putative MLLs (Figure 3B).

The highest likelihoods are summarized in Table 2, as well as the results for the Akaike estimates (ΔAIC and w _i). The model with the best fit for the complete (N=33) and the reduced (n=23) data sets was the bottleneck (Table 2). The inferred parameters for the bottleneck model and the 95% CIs obtained for the data sets with different numbers of individuals are shown in Supplementary Table S2. The parameter estimates under the empirical SFS suggests that population contraction occurred only a few generations ago (TBOT_compl=82; TBOT₃₃=10), with a continuous decline in population size after the bottleneck. However, the 95% CIs were wide and overlapping.

Despite the general expectation of low intrapopulation genetic variability and high susceptibility to inbreeding depression associated with pronounced effects of genetic drift in small populations (e.g., Kimura and Crow, 1964; Hamrick and Godt, 1989; Barret and Kohn, 1991; Ellstrand and Elam, 1993; Honnay and Jacquemyn, 2007), the results for M. catharinensis do not support the hypothesized loss of genetic diversity, which should reflect high levels of inbreeding, as a consequence of reduced population size and narrow endemism. Although contradictory to theoretical predictions, unexpectedly high levels of genetic diversity are not unusual for narrow endemics (e.g., Fernández-Mazuecos et al., 2014; Jiménez-Mejías et al., 2015; Turchetto et al., 2016; Forrest et al., 2017; Goetze et al., 2018; Sękiewicz et al., 2020; Bard et al., 2021; Garcia-Jacas et al., 2021), creating an indistinct pattern for this kind of plant species. For instance, in the studies by Fernández-Mazuecos et al. (2014) and Jiménez-Mejías et al. (2015), the authors highlight the paradox of genetic diversity levels in narrow and extremely narrow endemic plant species from the Mediterranean. While Fernández-Mazuecos et al. (2014) found moderate levels of genetic diversity for Naufraga balearica, Jiménez-Mejías et al. (2015) reported high levels of genetic diversity for Pseudomisopates rivas-martinezii. In comparison to previous genetic assessments for narrow endemics in the same regions, the authors (Fernández-Mazuecos et al., 2014; Jiménez-Mejías et al., 2015) found that, although some species have low genetic variability, most species showed moderate to high levels of genetic diversity (e.g., Sales et al., 2001; Coppi et al., 2008; Mameli et al., 2008; Mayol et al., 2012; De Castro et al., 2013; Forrest et al., 2017; Sękiewicz et al., 2020; Garcia-Jacas et al., 2021). Thus narrow endemism does not always imply limited genetic diversity. However, such contrasting results may be due to differences in population size, distribution range, ecological traits, and evolutionary history (Fernández-Mazuecos et al., 2014).

In addition, several studies have indicated that life-history traits (e.g., lifecycle, growth form, mating and breeding systems) strongly influence the amount and distribution of genetic variation in natural populations (e.g., Hamrick and Godt, 1996; Ellis et al., 2006; Elegren and Galtier, 2016; Goetze et al., 2018; Bard et al., 2021; De Kort et al., 2021). For Aechmea kertesziae, a narrow endemic species of Brazilian Coastal restingas (sandy vegetation), traits such as self-incompatibility, long-term persistence, clonal reproduction, and consistent population size may explain the high levels of genetic diversity observed for the species (Goetze et al., 2018). Specifically for M. catharinensis, the predominance of historical outcrossing, as shown in our analysis, seems to be an important factor that contributes to the moderate levels of intrapopulation genetic diversity. It is important to note that outcrossing plants are likely to present higher levels of intrapopulation genetic diversity than self-compatible species (Hamrick and Godt, 1996; Glémin et al., 2006; Ball et al., 2020). Moreover, a lack of successful sexual reproduction in M. catharinensis (i.e., production of pods without seeds) suggests that the level of intrapopulation genetic diversity has remained static for at least two decades, even though our demographic inferences indicated a recent decline in population size.

As life-history traits play a significant role in shaping the patterns of genetic diversity, it is important to compare the levels of genetic diversity among congeneric species in order to mitigate phylogenetic effects (Bevill and Louda, 1999; Simon and Hay, 2003). This is essential to better understand how features such as population sizes and geographic distribution ranges affect reproductive biology and, as a consequence, the genetic variability of closely related plant species. For example, population genetics studies indicate that Petunia secreta and Petunia exserta, both of which are narrow endemic species, have high genetic diversity indices; however, these indices are lower than those reported for its congener Petunia axillaris, which has a wider geographic distribution range (Turchetto et al., 2016). Beyond dissimilarities in the area of occurrence of these Petunia species, different floral syndromes and reproductive systems have been described as the probable causes of the observed genetic variation (Kokubun et al., 2006; Turchetto et al., 2016).

Although Mimosa is one of the largest genera from Mimosoideae (Barneby, 1991), there are few population genetics studies for species of this genus (Moreira et al., 2011; Pramual et al., 2011; Arruda et al., 2019; Araújo et al., 2020; Morales et al., 2020), and none of these previous studies used SNP markers. Nevertheless, genetic diversity has been estimated using markers such as RAPD (Mimosa pigra; Pramual et al., 2011), AFLP (Mimosa subser. Dolentes–Brevipedes; Morales et al., 2020), ISSR (M. caesalpiniifolia Benth.; Araújo et al., 2020), and allozyme (Mimosa scabrella; Moreira et al., 2011; Arruda et al., 2019). In this context, M. catharinensis presented lower levels of intrapopulation genetic diversity than that reported for its widespread congener M. scabrella (H_E =0.362 to 0.469; Moreira et al., 2011; Arruda et al., 2019).

When comparing studies on narrow endemics assessed using SNPs, and considering the theoretical maximum heterozygosity for biallelic markers (H_E =0.5), we are unable to make generalizations about the genetic diversity of narrow endemics. Indeed, M. catharinensis displays higher levels of genetic diversity than self-compatible (e.g., Ball et al., 2020) and clonal and functional sterile species (e.g., Amor et al., 2020). However, some tropical plant species (e.g., Lanes et al., 2018; Ball et al., 2020) showed the most similar genetic diversity to M. catharinensis.

Our findings indicate an excess of heterozygosity for this plant species, counteracting the theoretical expectation that reductions in population size will result in high levels of inbreeding. In terms of genetic variability, the levels of inbreeding in narrow endemic plant species with small population sizes vary according to individual features of each species/population (Angeloni et al., 2011). Reports for plant species with very small populations vary from high levels of inbreeding (Jiménez et al., 2014; Finlay et al., 2017; Rodrigues et al., 2019; Wang, 2020), to no evidence of inbreeding (Edwards et al., 2014, 2021; Spoladore et al., 2017), to an excess of heterozygotes (Cabrera-Toledo et al., 2008; Radosavljević et al., 2015; Amor et al., 2020) as found herein. Barret and Kohn (1991) emphasize that the mating pattern is a prime determinant of inbreeding levels in natural populations regardless of their size. In fact, self-compatible species can be more susceptible to inbreeding considering that self-incompatibility is likely to have evolved to prevent inbreeding depression (Charlesworth and Charlesworth, 1987). On the other hand, predominantly outcrossing species can suffer from inbreeding depression due to mating between relatives (Husband and Schemske, 1996).

The (historical) outcrossing rate observed in M. catharinensis suggests that it has a mixed mating system, which is characterized as ranging from 5 to 95% depending on environmental conditions and pollination frequency (Karasawa, 2015). Such a wide range in outcrossing rates for mixed-mating species are common and can change over time and space due to environmental conditions and intrinsic population features, such density and demographic structure (e.g., Wright et al., 2013). As with M. catharinensis, high outcrossing rates have been reported for many Neotropical plant species (e.g., Ward et al., 2005; Feres et al., 2012, 2021; Nazareno and Reis, 2012; Spoladore et al., 2017; Godoy et al., 2018; Montagna et al., 2018; Sujii et al., 2021), including two congener Mimosa species (Moreira et al., 2011; Arruda et al., 2020). Studies on populations of M. scabrella (Arruda et al., 2020), for instance, showed a similar reproductive pattern with high outcrossing rates (t=0.925/0.845) and low estimates of selfing (s=0.075/0.155). For M. catharinensis, the interpretation of negative F _IS values coupled with high t values should be viewed with caution due to the distinct biology of the species (i.e., production of seedless pods and absence of seedlings in the study area). Such biological characteristics observed for M. catharinensis suggest that a significant selective pressure against homozygotes must have occurred in previous generations, especially when taking into account the excesses of heterozygosity observed in its small population. As M. catharinensis population is constituted only by adult trees, a further study investigating inbreeding depression seems not to be a feasible task as a consequence of the lack of different ontogenic stages in its population.

Although pollen grains have been reported as viable for M. catharinensis (Silva et al., 2005), the apparent reproductive failure seems to be due to an extremely low frequency of floral visitors (Silva et al., 2005), implying low pollination efficiency. This phenomenon has been also described for other Mimosa species, such as Mimosa bimucronata (Seijo and Neffa, 2004). As with M. catharinensis, M. bimucronata does not produce nectar, however, its newly opened flowers exude a slight fruit odor to attract pollinators (Silva et al., 2011). Despite poor fruit set due to inefficient pollination, M. bimucronata produces viable fruits and seeds (Seijo and Neffa, 2004). In contrast, Mimosa pudica showed high pollination efficiency, mediated by A. mellifera, resulting in high fruit and seed yields and quality (Taimanga and Fohouo, 2018). Further pollination biology studies are needed to clarify the causes of reproductive failure in M. catharinensis.

The selection of the best demographic model was the main focus of the coalescent-based inference for M. catharinensis, considering a comparison of two data sets with different numbers of individuals and excluding putative clonal individuals based on their genetic distance. The fact that confidence intervals were wide and overlapping did not affect our results and can be related to the uncertainty of the estimated parameters and/or the small SFS (Rödin-Mörch et al., 2019). Nonetheless, the values presented in Supplementary Table S2 should be interpreted with caution. The generated SFS for both data sets suggest that a recent bottleneck is the best explanation for the current levels of genetic diversity in M. catharinensis. This corroborates our findings of excess heterozygosity for the species, which is a characteristic of populations that have experienced a recent bottleneck (Luikart et al., 1998). Of the factors that may have resulted in a drastic decrease in species population sizes in the restinga, the fragmentation of Santa Catarina Island’s native vegetation is particularly relevant, especially the degradation that occurred within the PAERVE before it was deemed a protected area. The implementation of reforestation projects in the park using exotic, invasive species also resulted in fires, changes to the soil, and deforestation of native vegetation, all of which accentuated the damage occurring to such a heterogeneous vegetal formation (Heberle, 2011).

In contrast to the evidence of population contraction for M. catharinensis, the majority of studies on past demographic changes in tropical plant species from South America have reported indications of expansion during the Quaternary period (e.g., Turchetto-Zolet et al., 2013). On the other hand, demographic stability has been suggested for the period for species from the Southern Brazil coastal restinga (Goetze et al., 2018). Therefore, further studies are needed to infer historical changes in population sizes of plant species exclusive to Santa Catarina Island (Hassemer et al., 2015), whose recent demographic changes may be similar to that reported herein for M. catharinensis.

The extremely small number of individuals (N=33) in the only known M. catharinensis population, coupled with strong evidence for a lack of sexual reproduction, suggests that the conservation perspectives for this critically endangered species are concerning. Although no inbreeding was detected, the species is facing imminent risk of extinction because of its reduced population size, making it more susceptible to stochastic events. As a matter of fact, ex situ conservation activities must be implemented in order to safeguard the remaining genetic diversity of M. catharinensis. One potential strategy suggested in studies on biodiversity conservation of threatened species consists of the preservation of tissue culture, an in vitro cultivation of isolated live tissue, enabling the propagation of species facing difficulties with natural reproduction (Paiva and Paiva, 2001; Santos et al., 2019). Protocols to apply tissue culture techniques for conservation of endangered species have been successfully developed (e.g., Sherif et al., 2018; Choudary et al., 2020; Lerin et al., 2021; Mishra et al., 2020), and such an approach has been used to conserve the narrow endemic species Styphelia longissima (Thomas et al., 2021).

In addition, as in vitro germination tests show pollen grain viability (Silva et al., 2005), additional pollination biology studies (e.g., pollination treatments such as manual cross-pollination and self-pollination) should be performed to test for viable seed production. The information available through pollination biology studies can help us to solve the puzzle of M. catharinensis reproduction, allowing us to better plan effective conservation measures for this rare plant species. In terms of in situ conservation, the fact that the species is found in a protected area is a significant starting point, although to date this fact has been insufficient for its conservation. Today, the management of invasive, exotic species (e.g., Pinus spp. and Eucalyptus spp.), one of the main issues threatening the flora of PAERVE, should be prioritized. Furthermore, it is widely reported that adequate restoration of the contaminated areas of the park are urgently needed (Ferreira, 2010; Bechara et al., 2013), and efforts to achieve this goal have been implemented according to the park’s Management Plan (IMA, 2020). Moreover, although a robust floristic survey of Santa Catarina State has been conducted (Vibrans et al., 2012), not all the biodiversity present in the state was sampled. As such, biodiversity inventories within Brazil’s protected areas are crucial (Oliveira et al., 2017), and efforts on field surveys must be undertaken to verify the narrow occurrence of this critically endangered plant species. Likewise, we are committed to monitoring the M. catharinensis population size over time to avoid further losses in genetic diversity.