Introduction

The genus Vanilla Plum. ex. Mill (Orchidaceae Juss.) consists of approximately 118 species worldwide, most remaining unexplored scientifically (Karremans et al. 2020). Perhaps because of this lack of information, only nine species have been included in the International Union of Conservation of Nature's List of Threatened Species (IUCN 2022). Vanilla planifolia Andrews, V. pompona Schiede, and the hybrid V. × tahitiensis J.W.Moore, are extensively studied due to their commercial value as sources of natural vanillin (Cameron 2011). However, wild and cultivated populations of V. planifolia are highly vulnerable to genetic erosion due to the over-collection, habitat destruction, and the clonal management of the crop (Bramel & Frey 2021, Soto Arenas 2006). Consequently, there is a need for conservation strategies and alternative vanillin sources (da Silva Oliveira et al. 2022, Flanagan et al. 2018). Vanilla chamissonis Klotzsch is a fragrant, hemiepiphytic orchid native to South America with terrestrial and climbing growth (Soto Arenas 2003, Soto Arenas & Dressler 2010). Taxonomically, it has been placed in the subgenus Xanata, section Xanata by Soto Arenas & Cribb (2010) and Besse et al. (2021), while others recognize only the section shared with 38 species with fragrant fruits (Karremans et al. 2020). Indeed, it has been recently demonstrated that V. chamissonis and V. bahiana Hoehne synthetize aromatic compounds with promising bioeconomic potential (da Silva Oliveira et al. 2022). Furthermore, a recent study unveiled the pollination strategy of Neotropical vanillas, including V. chamissonis. It was found that these species employ a pollination strategy that entices euglossine bees with nectar rewards (Pansarin 2022). On a physiological level, it has been shown that sunlight directly and proportionally affects the leaf area, succulence, and chlorophyll content (Lopes et al. 2019). The geographical distribution of V. chamissonis ranges from northern South America, including French Guiana, to southeastern Brazil, Bolivia, and northeastern Argentina (Biganzoli & Múlgura de Romero 2004, Jørgensen et al. 2013, Soto Arenas & Cribb 2010, Szlachetko 2016). In Brazil, it is among the most frequent Vanilla species found in the Atlantic Forest, extending from Espírito Santo to Rio Grande do Sul (Reis et al. 2011). In Argentina, V. chamissonis is restricted to specific areas of the province of Misiones, with recorded occurrences in Isla Caraguatay Provincial Park (located in the Paraná River), the Teyú Cuaré Provincial Park and the adjacent private Osununú Natural Reserve (Reserva Natural Osununú; RNO) (Biganzoli & Múlgura de Romero 2004, Munno et al. 2011, M. A. Munno et al. unpub. data). Due to the ongoing fragmentation of natural habitats and the restricted distribution of V. chamissonis, this species has been classified as threatened in Argentina, listed on the Appendix II of Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES 2022). Therefore, a significant project is underway at the RNO to conserve and sustainably manage V. chamissonis. The objective of the present work was to genetically characterize the wild unexplored germplasm of V. chamissonis from the RNO. This step represents the initial phase in identifying priority areas for preserving locally adapted variants. To date, such genetic characterization has not been carried out. Additionally, the genetic data obtained could be valuable for future breeding programs. Molecular genotyping techniques allow revealing the genetic background of organisms through the scrutiny of highly polymorphic loci. These techniques have been applied in studying several plant species, including many orchids (George et al. 2009, Hou et al. 2012, Ma & Yin 2009, Qian et al. 2013, 2014, Verma et al. 2009). Insights on population structure and species identification are basic knowledge in defining the candidate units that merit in situ conservation. Notably, Vanilla species from Argentina have been documented under various names (Biganzoli & Múlgura de Romero 2004, Johnson 2001, Keller et al. 2019, Munno et al. 2011, M. A. Munno et al. unpub. data), many of which have been synonymized (e.g., Soto Arenas & Cribb 2010). Therefore, as part of the characterization of this genetic resource, we undertook its molecular identification within a phylogenetic context.

Materials and methods

Plant material.-- The plant material was collected at the RNO in May 2015 and May 2019. Within the RNO, these climbing plants develop within sectors called nuclei. Four of the nuclei reached were outlined in the work of M. A. Munno et al. (unpub. data) (i.e., N1 to N4), and the additional nuclei were pinpointed by the RNO's personnel, who identified these sectors during their vigilant monitoring activities. Given that these plants may reproduce both clonally and sexually (Munno et al. 2011, M. A. Munno et al. unpub. data, Reis et al. 2011), each plant (stem) is considered -a priori- as an individual. The plants sampled from each nucleus are treated as a population (Appendix 1). Eight nuclei were reached, with a variable number of individuals collected from each nucleus. Leaf samples were collected from healthy and accessible plants at least 1.5-2.0 m apart, ensuring they were not interconnected, and were then preserved in silica gel. Out of the 72 individuals, those that dried without any apparent signs of oxidation were employed for total genomic DNA extraction (Appendix 1). This extraction process involved processing small leaf pieces with the DNeasy Plant Quick Extraction kit (QIAgen Inc., Duesseldorf, Germany) following the manufacturer’s instructions.

Sampling site. -- The private reserve RNO is located at the Municipio and Department of San Ignacio (Misiones, Argentina), and is adjacent to the Provincial Park Teyú Cuaré (Fig. 1); both the park and reserve are within the Paranaense Forest and show remnants of the Cerrado biome. The reserve involves 174 ha of undulating topography with heights up to 215 m; five distinct environments or landscape units were delineated within the RNO, encompassing wetlands, vegetation areas associated with exposed rocky outcrops or elevated sandstone formations, mixed forests, secondary successions, and forests situated in lowlands and slopes (Keller et al. 2019, Velazco et al. 2015). The mixed forest encompasses most of the area (79 ha). However, the Vanilla under study is associated with the arboreal vegetation of the Myrtaceae forests found in both lowlands and slopes. Additionally, it is associated with the Annonaceae trees located on the sandstone rocky outcrops slopes, where the soil is characterized by nutrient deficiencies (Keller et al. 2019, M. A. Munno et al. unpubl. data). The RNO exhibits a subtropical, humid temperate climate without a clearly defined dry season. The mean annual temperature is 21.86 ºC, while the mean annual rainfall is 1644 mm. The rainiest period in the region comprises October to February.

Molecular genotyping with co-dominant markers.-- Six heterologous simple sequence repeat (SSR) markers, originally designed for V. planifolia and reported showing amplification signals in V. chamissonis and V. bahiana by Bory et al. (2008a), were assayed on eight randomly chosen DNA samples. The SSR loci assayed were: 005, 015, 019, 025, 031 and 047. PCR amplifications were performed using 30-60 ng of genomic DNA, 60 ng of each primer, 1X PCR buffer (Invitrogen Thermo Fisher Scientific, California, USA), 2 mM MgCl₂, 0.2 mM of each dNTP and 0.5 unit of Taq DNA polymerase (Invitrogen Thermo Fisher Scientific), in a final volume of 25 µl. PCR cycles were as in Cascales et al. (2014), but the annealing steps were tested at 50°C, 55°C and 60°C. PCR products were checked by means of electrophoreses in 1.5% (w/v) agarose gels, stained with ethidium bromide, and photographed under UV light. Promising loci were amplified for all individuals and subjected to high-resolution electrophoreses. For this, SSR amplicons were mixed with an equal volume of 98% (v/v) formamide, heat-denatured, ice-cooled, and loaded into a 6% (w/v) polyacrylamide gel under denaturing conditions (8 M urea). Electrophoreses were carried out in 1X TBE buffer (0.09 M Tris, 0.09 M boric acid, 2 mM EDTA, pH 8.0) at 60 W for 4-5 hs. Bands were visualized by silver nitrate staining using the Silver Sequence Staining Reagents kit (Promega, USA) according to the manufacturer’s instructions. Air-dried gels were then digitalized, and band sizes were estimated from gel images. A molecular weight standard (30-330 AFLP DNA Ladder; Invitrogen Life Technologies) was included in every run as a size reference and staining control.

Molecular genotyping with dominant markers.-- A set of Inter-Simple Sequence Repeat (ISSR; Zietkiewicz et al. 1994) markers were screened on 13 randomly chosen genomic DNA samples; the primers assayed were: ISSR-P1: (AG)₈Y; ISSR-P3: CT(GA)₈; ISSR-P4: (GA)₈C; ISSR-P6: (GACA)₄; ISSR-P7: (GA)₈T; ISSR-P9: (GATA)₄; ISSR-P13: (AG)₈YT. PCR conditions were as in Paiva et al. (2020). Primers that produced more than five clear bands, between 300-1500 bp in size in agarose gels, were selected for further analyses. Aliquots of ISSR amplifications were prepared as above, separated through denaturing 5% (w/v) polyacrylamide gel electrophoreses, and visualized through silver staining procedures, as already stated. Fragment sizes were estimated from digitalized gel images by comparison against molecular weight standards (100 bp DNA Ladder -Productos Bio-lógicos, Argentina- and 30-330 AFLP DNA Ladder). For a given locus, each band was registered as present (1) or absent (0). These data were used to construct five binary matrices, one per ISSR primer plus a combined matrix comprising individuals with >50% of loci information. In addition, low-frequency bands (≤ 0.10) were removed to reduce errors generated during scoring. Genotyping error rates were determined for each ISSR marker by running duplicate, independent PCR amplifications for 4-5 individuals; errors were estimated as the ratio of the total number of loci with discordant scores (all individuals combined) to the product of the number of individuals by the total number of scored loci (Bonin et al. 2004).

Analyses of genotyping data.-- The performance of the markers was evaluated through the polymorphism information content (PIC) calculated as in Roldán-Ruiz et al. (2000). Also, the power of the set of loci to discriminate between unique individuals (referred to as multilocus genotypes) was assayed through the genotype accumulation curve, with the function genotype_curve of the poppr package v. 2.9.4 (Kamvar et al. 2014, 2015) implemented in R4.3.1 (R Core Team 2023), setting the resampling to 1000. Partitioning of genetic variation within and among pre-defined populations was assessed by the analysis of the molecular variance (AMOVA; Excoffier et al. 1992) using GenAlEx 6.5 (Peakall & Smouse 2012); statistical significance of each variance component was assessed based upon 999 permutations in GenAlEx. To identify genetic clusters without any a priori group delimitation (for instance, populations), we conducted Discriminant Analysis of Principal Components (DAPC; Jombart et al. 2010) on the five matrices, using the adegenet 2.1.0 package (Jombart & Collins 2015) implemented in R 4.2.1 (R Core Team 2022). Genetic clusters were assessed using the find.clusters function with n.iter=100,000 and n.start=10, evaluating a range from 1 to 6, by the K-means algorithm (Legendre & Legendre 1998). The optimum K value and number of discriminant functions were determined for each matrix.

Relationships among individual plants were investigated, as well, through a clustering distancebased analysis for the combined matrix. Distances were obtained applying the complementary of the Dice similarity index, in FAMD 1.31 (Schlüter & Harris 2006); a Neighbor-Joining unrooted phylogram was obtained and the bootstrapping was performed with 1000 pseudoreplicates with the same program. Resulting topology was visualized with FigTree v1.4.3 (Rambaut 2016).

A correlation analysis between pairwise genetic and geographical distance matrices was performed through a Mantel test in the Alleles In Space program (AIS; Miller 2005). The spatial genetic structure was approximated through the Genetic Landscape Shape interpolation procedure, as implemented in AIS, employing the geographical distances among the nuclei and the residual genetic distances among individuals to construct a connectivity network among the sampled points based on Delaunay triangulation. The X and Y coordinates are the midpoints of each edge in the triangulation, and the ‘heights’ in the landscape are reflections of the genetic distances (i.e., residual genetic distances) between the observations found at the vertices of the triangles. An inverse distance weighted interpolation (weighting parameter= 1.5) across a uniform grid size of 100 × 100, was set.

Phylogenetic analyses.-- Nucleotide sequences of the entire nuclear Internal Transcribed Spacer region of the ribosomal DNA (namely, ITS1-5.8S-ITS2) and partial sequences of plastid genes coding for the Rubisco large subunit (rcbL) and the maturase K (matK) were obtained from three individuals representing different plant sectors or nuclei (namely, N1P374, N3P397, and N7P376). The ITSs were generated using primers ITS-4 and ITS-5 (White et al. 1990); for rbcL the Z1 and Z1351R primers (Soltis et al. 1992) were employed, and for matK the forward 5'-CGATCAACATCTTCTGGAGTGT-3'and reverse 5'-TGAACAGTTATGAATGAGCCACT-3' primers, which were designed in this study with Primer3web (Untergasser et al. 2012). The sequencing service was provided by Macrogen Inc. (South Korea). Chromatograms were inspected in BioEdit (Hall 1999), and the consensuses between forward and reverse sequences were obtained with the same program.

The ITS sequences obtained here (OQ080023- OQ080025) were analyzed together with 77 ITS sequences retrieved from the Genbank belonging to 55 species of Vanilla, and 11 sequences from three genera that were used as outgroups (Appendix 2). Plastid genes were analyzed separately from the nuclear partition; the sequences generated here for rcbL and for matK (OR472978-OR472983) were analyzed with 38 Vanilla accessions derived from 28 species and two Epidendrum species used as outgroups (Appendix 2).

Sequences were aligned in MEGA 7.0.25 (Kumar et al. 2016) using the MUSCLE algorithm (Edgar 2004) with default parameters. A pairwise p-distance matrix was generated in MEGA to survey for duplicated sequences. Evolutionary model selection and maximum likelihood (ML) phylogenetic analyses were conducted on IQ-TREE v. 2.0.6 (Minh et al. 2020, Nguyen et al. 2015). Each segment of the ITS, namely ITS1, 5.8S, and ITS2, and each plastid gene had a distinct substitution model, selected through the ModelFinder algorithm (Kalyaanamoorthy et al. 2017) in IQ-TREE; when needed, the best partition scheme was selected with the same program. The command path line for the nuclear matrix was: Iqtree2 -s inputfile.fas -keep-ident -p partition_file.txt -m TESTMERGE -pre outputfile-bnni -nt AUTO. The command line was as above for the plastid genes, but the command -s was excluded and -p was used to call both matrices. Thus, the best model determined for ITS1 and ITS2, was TamuraNei with unequal base frequencies (F) and site rate heterogeneity accounted through a Gamma distribution with four categories (G4) that is, TN+F+G4 (alpha parameter= 1.91); the 5.8S partition received the model K2P+G4 (Kimura two parameters, alpha parameter= 0.26). The fittest model selected for rbcL was K3Pu+F+I (Kimura three parameters with unequal base frequencies and a proportion of invariant sites, I=0.8) and for matK was TVM+F+G4 (transversion model, alpha parameter= 0.56). Ten independent runs were carried out using the same multiple sequence alignment obtained for each marker, but with different orders of the taxa, to escape from local optima. In each run, support for internal nodes was estimated by ultrafast bootstrapping (Hoang et al. 2018) and by the Shimodaira-Hasegawa approximate likelihood ratio test (SH-aLRT; Guindon et al. 2010) with 10,000 pseudoreplicates. The command path line used for each ITS's run was: Iqtree2 -s inputfile_runX.fas -keep-ident -p partition.best_scheme.nex -pre outputfile_runX -bnni-bb 10000 -wbt -alrt 10000 -nt AUTO. The command -p inputfile_runX.fas was used for running the plastid markers, and the command -p partition.best_scheme. nex was excluded. Trees were visually inspected with Figtree. Afterwards, to detect differences among runs, the topology test of Shimodaira and Hasegawa (1999) and the Approximately Unbiased test (AU; Shimodaira 2002) were performed for ITS, with the command line: Iqtree2 -s inputfile_runX.fas -p partition.best_scheme. nex -z 10_trees4test.trees -zb 10000 -zw -au -pre output_Tree_test -nt AUTO. For the plastid genes, the command line was adjusted as outlined above. A full bootstrapping was performed after the best runs were selected (with command: -b 10000). These took 96 hr for ITS and 16 hr for the plastid genes, using 20 CPU cores. The final phylograms were rooted with Epidendrum species, following Pérez-Escobar et al. (2021), and visualized in Figtree.

Results

Ninety-three percent of the plants gathered yielded genomic DNA of good quality (N=67). Of the six heterologous SSR markers assayed, solely markers 019 and 031 yielded banding patterns compatible with co-dominant loci. The other four markers produced several bands (between 100-1000 bp in size) even under stringent amplification conditions (i.e., 60°C), and thus were excluded. When markers 019 and 031 were assayed in high-resolution gels, no variation was evidenced among >60 individuals tested. Then, the plants were characterized through ISSR genotyping. Four of the seven examined primers exhibited reliable banding patterns (Appendix 3) and were successfully amplified across 58-65 individuals (Table 1). Out of 296 bands registered, 54% were excluded due to their appearance at a frequency below the preestablished threshold (i.e., <0.10). Each primer produced high PIC values (>0.43; Table 1), and solely for primers ISSR-P3 and ISSRP7 it was detected a genotyping error (1.0-1.2%). The genotype accumulation curve displays that within the conditions of this study, complete discrimination (100%) of multilocus genotypes can be achieved by utilizing a minimum of 27 ISSR loci (see Appendix 4). The AMOVA results showed that most of the genetic variation appeared within populations ranging from 77% to 84% for single ISSR primer datasets and reaching 81.3% for the combined matrix (Table 1). The DAPC discriminated between four to six genetic clusters, according to the matrix considered, without any prior group indication (Table 1). The combined matrix showed that the individuals were assigned to four nonoverlapping genetic clusters or gene pools (numbered 1-4; Fig. 2A), consisting of 11, 12, 23, and 13 plants, respectively. DAPC results also evidenced that 81% of the individuals showed >0.90 probability of membership to the cluster assigned, that is, most individuals have in their genetic makeup predominance from a single gene pool (Fig. 2B). Solely 11 individuals were ambiguously assigned, exhibiting membership probabilities ranging from 0.483 to 0.815, showing signals of genetic admixture (Fig 2B; Fig. 3, pie charts). Then, for the plants from nuclei N1, N3, and N4, at least three distinct gene pools have been estimated to contribute with high probability. In contrast, for the plants from N5, N7, and N8, the estimation indicates a genetic input from two pools with high probability. For N2, there is a predominant assignment to a single gene pool (Fig. 1, pie charts; Fig. 2B).

The average genetic distance estimated among all individuals equals 0.3648. Plants from nucleus N1 yielded the highest average distance (d=0.4197), and those from N3 showed the smallest distance (d=0.2159). The population comparisons indicated that N1-N8 are the most dissimilar pair (d=0.4404), being linearly among the closest nuclei (154 m; Appendix 5). In contrast, N2-N3 appear as the most similar nuclei (d=0.2725), being 388 m apart. The unrooted neighbor-joining dendrogram showed groupings concordant with DAPC clusters 1, 2, and 3; while individuals assigned to cluster 4 appeared dispersed (Fig. 3). Still, very few nodes received bootstrap support values worth considering.

A slight but significant positive correlation between genetic and geographic distances was detected (r= 0.1724; P< 0.001). The genetic landscape shape interpolation analysis showed that the highest positive peaks, i.e., the greatest genetic distances, are found among individuals from the northeastern region of the RNO (Appendix 6). Instead, valleys and negative peaks are apparent among the nuclei located towards the southwest. This landscape reveals a genetic discontinuity between both regions.

The alignment of 88 ITS sequences spanned 755 bp, of which 31.5% were constant sites. ITS sequences obtained from plants N1P374, N3P397, and N7P376 resulted identical. Likewise, the Genbank ITS sequences of V. odorata (accession numbers MN902028 and MN902032) on the one hand, and of V. mexicana Mill. (accession numbers MW829691 and MW829692) on the other, also resulted in being identical. As a result, only a single representative was subjected to further analysis for each of these species. Thirteen accessions belonging to different nominal taxa showed identical ITS sequences. Specifically, the pair V. chamissonis (MW829674) and V. bahiana (EU498163); the quartet V. pompona (GQ867234), V. pompona subsp. pittieri (Schltr.) Dressler (MN902062), V. pompona subsp. grandiflora (Lindl.) Soto Arenas (MN902058), and the hybrid V. × tahitensis (MN902063); and the triplet V. planifolia (U66819), V. pompona subsp. pompona (MN902053) and V. cf. planifolia (FJ425832). The pairs V. norashikiniana Go & Raffi (MH777755) and V. griffithii Rchb.f. (MH777726), as well as V. shenzhenica Z.J.Liu & S.C.Chen (JF796930) and V. somae Hayata (KY966687), also had identical ITS sequences but were kept for analyses. The ten independent runs yielded phylogenetic trees showing slight topological differences, although none in the clade containing V. chamissonis; statistical tests failed to discriminate among the topologies (SH p-value ≥ 0.381; AU pvalue >> 0.05). The phylogram from run_7 was selected based on the highest logL, and lowest BIC and AIC values (Appendix 7); it shows that all Vanilla sequences formed three highly supported clades (Fig. 4), with their relationships poorly supported (30-64% bootstrap support value, BSV). One highly supported clade (94% BSV) comprises 44 sequences representing 18 taxa from North, Central, and South America, including V. × tahitensis and V. hirsuta M.A.Clem. & D.L.Jones from Oceania. The ITS sequences obtained in this study formed a monophyletic group (100% BSV) with V. chamissonis MW829674, V. bahiana EU498163, and the three sequences of V. calyculata Schltr. A second highly supported clade (94% BSV) involved 26 sequences from 20 Asian or African species plus three species from Central America (V. barbellata Rchb.f., V. claviculata Sw., and V. dilloniana Correll) whose position within this clade is uncertain. The third clade involves three American taxa (98% BSV). Noteworthy is the placement of V. mexicana (MN902027) and V. planifolia (MN902045); these accessions appeared nested among Epidendrum species, with full support (100% BSV).

The alignments of the plastid genes for 28 Vanilla species and two Epidendrum L. species (43 sequences, in total; Appendix 2) reached 1264 bp for rbcL and 840 bp for matK, and exhibited 93% and 79.3% of constant sites, respectively. The sequences obtained here for plants N1P374, N3P397, and N7P376 resulted identical for both genes. Among the other nucleotide sequences included, those attributed to V. bahiana from voucher specimens CR0099 and CR0668, exhibit identical sequences for both rbcL and matK; as well, the rbcL of V. bahiana is indistinguishable from those of V. insignis and V. phaeantha. Similarly, V. madagascariensis Rolfe and V. humblotii Rchb.f., on the one hand, and the five accessions of V. planifolia and V. sotoarenasii M.Pignal, Azof.-Bolaños & Grisoni (CR0068), on the other, show identical plastid sequences. No discernable topological variations regarding sister-group relationships among the ten phylogenetic trees were identified. Furthermore, the statistical tests did not discriminate among the topologies (SH p-value ≥ 0.224; AU pvalue > 0.05). The phylograms from run_0 and run_6 emerge equally likely based on their highest logL and lowest BIC and AIC values (Appendix 7). The plastid phylogram (Fig. 5) depicts a highly supported clade (98% BSV) comprising 12 taxa (25 sequences) from North, Central, and South America, and V. × tahitensis ; within this clade, the sequences obtained in this study relate to the reference voucher V. chamissonis CR0666 (81% BSV). Another strongly supported clade (91% BSV) involves 15 Asian and African species that conform the sister group of the latter clade (92% BSV). Then, with full support, V. mexicana CR2144 appears as the sister taxon to all other vanillas. As expected, specimens with identical sequences formed clades with zero-length branches.

Discussion

The present work contributes basic and novel genetic knowledge for the understudied yet valuable germplasm of the fragrant orchid Vanilla chamissonis. The findings offer constructive insights for conserving and managing this species within the Reserva Natural Osununú (RNO), given that the natural habitat of V. chamissonis has been under anthropic reduction, mainly because of progressive deforestation. The utilitarian expression “nucleus”, used to denote distinct groups of V. chamissonis plants, emerges as a conceptually relevant unit in planning management strategies. The RNO comprises a remnant of the Cerrado biome, unique within Argentina and intermingles with rupiculous vegetation; its biological diversity harbors a set of distinct local adaptations that serve as valuable genetic reservoirs. For instance, 23 to 38 plant species have been reported with a distribution constrained to the geographical boundaries of the Teyú Cuaré Provincial Park and the Reserva Natural Osununú (Keller et al. 2019, Velazco et al. 2015). To genetically characterize the V. chamissonis germplasm, firstly, we attempted the transferability of SSR markers that were documented as successful in this species (Bory et al. 2008a); however, our efforts failed to reproduce their results on V. chamissonis. Four of those SSR markers yielded banding patterns incompatible with microsatellite loci, even under highly stringent amplification conditions, and the other two markers produced monomorphic patterns in high-resolution polyacrylamide gels. Gigant et al. (2012) developed 19 microsatellite markers for V. humblotii and V. roscheri Rchb.f., successfully tested on 14 individuals of four American species. However, it is noteworthy that none was V. chamissonis. Then, Gigant et al. (2016), in trying to transfer SSR markers from the African vanillas to V. mexicana, had a mild success, ending with four reliable markers. Additionally, their attempt to transfer 14 SSRs of Bory et al. (2008a) from V. planifolia to V. mexicana proved unsuccessful. Since species-specific co-dominant markers are an ideal tool for conducting populationalbased studies, developing novel homologous markers could significantly enhance such studies across the distribution range of V. chamissonis.

To gather genetic data that could unveil the genetic diversity, we genotyped the plants by applying anonymous dominant markers at hand, employing highly resolving and sensitive procedures, and following stringent inclusion criteria. The unsupervised DAPC genetic analysis allowed the detection of at least four gene pools that contribute to the genetic diversity of these plants. Also, each nucleus of vanilla plants varies in its genetic constitution, so that most of the total genetic variation (81%) was due to differences among individuals within each nucleus. Even though small populations of V. chamissonis have been reached in this study, present results suggest that future efforts should focus on preserving as many individuals as possible from each nucleus to maintain the existent diversity. Vanilla species form small populations in Colombia due to infrequent flowerings and rare pollination events; there is also a dependency on the forest dynamics' (Flanagan et al. 2018). Studies on the floral biology of V. chamissonis from Brazil revealed a spontaneous fruit production rate of 21.2%; this rate was deemed high, considering the inability to ascertain the pollinating agents (Reis et al. 2011). At the RNO, the modestly sized populations of V. chamissonis exhibited reproductive strategies involving regeneration through seed germination and agamic reproduction (M. A. Munno et al. unpubl. data). Although that study needs formalization, at the time, those authors examined 30 plants from four nuclei, of which only nine showed ripened pods. This suggests a higher fruit production rate (namely, 30%) than the rate documented by Reis et al. (2011). Thus, other factors may well be restricting the establishment of larger populations at RNO. For instance, the absence of suitable soil microbial communities wherever a seed is deposited (Mújica et al. 2009). The present study exposed the heterogeneity of the genotypes surveyed (yielding an average intraspecific genetic distance= 0.3648). This distance value is considerably higher than that reported by Verma et al. (2009), which estimated an overall interspecific mean distance of 0.45 among eight Vanilla species using ISSR markers. In addition, our investigation did not find evidence of vegetatively reproduced individuals (clones), as identical genotypes were not detected. Therefore, present results would suggest a predominance of sexually reproduced individuals, at least among the ones sampled so far.

Also, our survey allowed ascribing most individuals to one of the four distinct genetic clusters, or gene pools, with a high probability (>0.90). From a conservation standpoint, the nucleus N1 is an interesting unit because it receives input from three gene pools (DAPC clusters 1, 3, and 4); a fourth genetic input is also present but in a much lesser proportion (DAPC cluster 2; Fig. 1, 3). Furthermore, our research has revealed that a genetic discontinuity could exist between the nuclei in the northeastern sector of the RNO, where N1 is located, and the populations placed towards the southwest of the reserve (Appendix 6). Among the latter, the germplasm from N3 and N4 appear constituted by a different combination of three gene pools, and, therefore, its preservation might also be deemed valuable. As per the floristic survey of Keller et al. (2019), the northeastern RNO is dominated by a Myrtaceae forest, in which Eugenia psidiiflora O.Berg abounds. As well, the landscape features riparian “ararí” forests (Calophyllum brasiliense Cambess.), “pacurizales” (Garcinia brasiliensis Mart.), and rupicolous shrubby vegetation, although these are present in a much lesser extent. Moving southwards of RNO, the prominent plant community shifts to a xylopial forest (Xylopia brasiliensis Spreng.), as documented by Keller et al. (2019). Our findings suggest that sexual reproduction has shaped the observed genotypes among the examined plants. The distribution of this genetic variation may be attributed to pollen movement and ectoand endozoochorous seed dispersal (Karremans et al. 2023, Pansarin 2021, 2022, Pansarin & Suetsugu 2022, Reis et al. 2011), which could be compatible with the tenuous isolation by distance and the low genetic differentiation (PhiPT =0.187) evidenced among the populations. Thus, the present results highlight that the geographical proximity of V. chamissonis nuclei may not accurately reflect the plants’ genetic closeness. Similar results were obtained for outcrossing Zingiber species (Huang et al. 2019), and for the primarily allogamous populations of seven wild Vanilla species from Madagascar (Andriamihaja et al. 2021), among other cases. According to Bory et al. (2008a), even rare sexual events have the potential to contribute to the genetic diversity within Vanilla, and diversity may play a role in succeeding by settling vegetatively. Further investigation is required to elucidate how the unique organismal communities (for instance, bats) and the distinctive topography jointly found at the RNO have shaped the historical and current distribution of V. chamissonis.

To attain the molecular identification of the plants under study, we employed the biparentally inherited nuclear ITS region and the maternally inherited plastid genes rbcL and matK. Pérez-Escobar et al. (2021) found high congruence between nuclear and organellar phylogenetic inferences in orchids. Nevertheless, due to interspecific hybridization and polyploidy in Vanilla (Hu et al. 2019), we carried out separate phylogenetic analyses to keep track of both evolutionary histories. Flanagan et al. (2018) highlighted the importance of the ITS region and the matK as barcoding loci, particularly for this genus. More recently, Besse et al. (2021) demonstrated that the ITS region resulted far more diverse than the plastid markers tested (namely, matK, psaB, rbcL, or psbC), and the most appropriate phylogenetic marker for the genus Vanilla. The ITS phylogenetic analysis presented here intertwined sequences derived from the seminal works of Soto Arenas and Dressler (2010) and Besse et al. (2021) (32 and 18 accessions, respectively), which, to the best of our knowledge, were not analyzed jointly before. The ML phylogram shows that the orchids from the RNO share a most recent common ancestor with V. chamissonis reference voucher CR0666 (Besse et al. 2021), V. calyculata (=V. columbiana Rolfe; Karremans et al. 2020) and a V. bahiana representative (Pansarin et al. 2008), with full support. Still, a 100% identity was evidenced between the Genbank sequences used here for V. chamissonis and V. bahiana. Likewise, other ITS accessions showed 100% sequence identity and formed well-supported clades. Another clade was formed by seven accessions of V. pompona, V. × tahitensis, and three V. planifolia. As to the latter three accessions, when a Blast search was performed, it showed that, in each case, the best 30 hits matched with V. pompona (96.53-99.41% identity; e-value=0.0). Similar results were obtained for the accession of V. × tahitensis. Noteworthy is the placement of V. mexicana MN902027 and V. planifolia MN902045 within the Epidendrum clade, as evidenced in our phylogenetic analyses. These “vanilla” ITS sequences were derived from Soto Arenas and Dressler (2010), where these suspicious sequences appeared as the sister group of the remaining Vanilla sequences since solely one Epistephium Kunth and one Lecanorchis Blume sequences were used as outgroups. The calling for paralogs as an explanation for that location (Soto Arenas & Dressler 2010) seems unlikely given that Blast searches with MN902027 and MN902045 as queries retrieved the first 27 best hits with Epidendrum accessions (99% identity; e-value range: 2e^-170-1e^-152). An unnoticed mislabeling or even contamination with Epidendrum DNA could explain both results. The vouchers of these two “vanilla” accessions should be revisited to clarify their identification. Hence, it appears that several vanilla ITS accessions deposited in Genbank require careful curation.

Our findings regarding plastid markers are consistent with Besse et al. (2021) research, demonstrating that the variation exhibited by matK surpasses that of rbcL. In line with that work, the ITS sequences displayed the highest nucleotide variation among the three markers in the present study. As to the plants studied here, the plastid markers indicate that the organellar genome of V. chamissonis from Argentina may share a most recent common ancestor with the reference voucher's sequences of V. chamissonis. The species compared comprise several lineages of maternally inherited plastid genes. Among these, the North, Central, and South American lineages appeared monophyletic, as do African and Asian species' plastid genes; this clade conforms the sister group of the latter, while V. mexicana appears as the sister taxon to all other vanillas. This general pattern of relationships resembles that of Bouetard et al. (2010) based on four plastid genes and that of the ITS phylogram presented here, supporting the findings of Pérez-Escobar et al. (2021).

It is known that intraspecific phenotypic diversity, interspecific hybridization, and polyploidy prevail among species of Vanilla, possibly causing difficulties in species recognition (Bory et al. 2008b) and the instability of the taxonomic names associated with sequence data. This, of course, applies to V. chamissonis and allied species. Karremans et al. (2020) emphasized the necessity for in-depth investigations into living materials of V. argentina Hicken, V. carinata Rolfe, V. chamissonis, and V. vellozii Rolfe. These studies are crucial for clarifying their identity, relationships, and distribution. Similarly, further studies are needed for V. edwallii Hoehne and V. rojasiana Hoehne (=V. angustipetala Schltr.), previously reported in Northeastern Argentina (Johnson 2001, Soto Arenas & Cribb 2010). Unfortunately, sequence data is presently unavailable for the referred species, except for V. chamissonis (Genbank accessed August 2023).

Following a conservative approach, we deem it appropriate to continue considering the RNO-studied germplasm as V. chamissonis until more molecular data can be gathered from related species. This proposal is consistent with Biganzoli & Múlgura de Romero (2004), Munno et al. (2011), and M.A. Munno et al. (unpub. data), yet it diverges from Keller et al. (2019), who mentioned this orchid as V. vellozi. We strongly support Karremans et al. (2020), who highlighted the negative effects of inaccurate species identification on various biological aspects. This includes phylogenetic studies, species distribution and abundance estimation, and conservation and biodiversity assessment, all of which are affected.

Acknowledgements

The authors are grateful to the Ministerio de Ecología y Recursos Naturales Renovables from Misiones province for providing collection permits number 099/14 and 01/19, and Fundación Temaikèn for providing logistic and financial support during field trips. Special thanks to A. Ávalos, J. Guastavino, M. P. Bertolini, A. de la Colina, R. Flores, A. Rodríguez, R. Ramírez, and J. Torrezin, from the Fundación Temaikèn, for their kind assistance at RNO during field work and transportation; E. Stolle and Y. Wurm, from the Museum Koenig, Leibniz Institute for the Analysis of Biodiversity Change, Center for Molecular Biodiversity Research (Bonn, Germany) for their generous donation of supplies. Financial support from UBACyT 20020170100614BA is gratefully acknowledged. We thank the developers of Sci-hub for allowing unrestricted access to all scientific knowledge. The authors have no conflicts of interest to declare that are relevant to the content of this article.

Author Contributions: Carolina I. Paris designed and prepared the sample collection and map; Jimena Cascales and Alexandra M. Gottlieb conceived and performed the experiments, analyzed the data, and prepared all figures; Alexandra M. Gottlieb wrote the original draft. All authors read, edited, and approved the manuscript.

Funding: Partial financial support was received from Fundación Temaikèn (a non-profit organization). The research leading to these results received funding from Secretaría de Ciencia y Técnica, Universidad de Buenos Aires (UBACyT) under Grant No 20020170100614BA. The Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) and the Universidad de Buenos Aires (UBA) paid the salaries.

Literature cited

Publication Dates

History