Young leaves from 107 individuals of A. colubrina (adults and seedlings) were collected from four different Argentinean ecoregions: Paranaense forest (n = 95), Yungas (n = 6), Humid Chaco (n = 2), and the Delta and Islands of the Paraná River (n = 4). In the southern region of the Paranaense forests (Santa Ana, Misiones; -27.43372198, -55.579419), where A. colubrina characterizes forest patches within grassland landscapes, two life stages were sampled: 31 reproductively mature trees and 64 seedlings resulting from the germination of seeds collected from the fruits of four mother trees. The seedlings thus represent different sample sizes of four half-sib families and may contain full sibs since a previous study of an A. colubrina population in the same region suggested a selfing rate of 51–56% estimated from the inbreeding coefficient s = 2F_IS/(1 + F_IS) (Hartl and Clark 2007; Goncalves et al. 2019). Trees were georeferenced by GPS (Geographic Position System) using a Garmin eTrex® 20× receiver (precision of ± 3 m) and identified by an individual code. Leaves were dried with silica gel in the field and stored at room temperature.

Total genomic DNA was extracted for each individual using the modified cetyl-trimethylammonium bromide (CTAB) method (Doyle 1991). We added 15 mg of dry leaves and two steel beads (4 mm in diameter) into a 2 mL Eppendorf tube. Liquid nitrogen was used to break the sample into a powder using a Geno Grinder 2010 m (Spex Sample Prep, New Jersey) tissue homogenizer. Then, we added 400 µL of CTAB buffer and put the tubes into a 65°C water bath for lysis. The subsequent steps followed Doyle (1991). The DNA quality was assessed using a Nanodrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, USA) showing that the genomic DNA concentration was at least 80 ng/µL.

SSRseq markers were developed from low-coverage shotgun sequencing of a single library, prepared using the Qiaseq FX DNA library kit (Qiagen, Hilden, Germany). This library was generated from four pooled samples, each representing a different ecoregion (YS113, F134, IT185, SB5). Sequencing was conducted using Illumina MiSeq v.3 (Illumina, San Diego, USA) 2 × 300 bp paired-end sequencing, generating 4,497,218 read pairs. Overlapping forward and reverse reads were merged using BBMerge v.38.87 (Bushnell et al. 2017) with a minimum quality of 25, an overlap greater than 100 bp without any authorized mismatch, resulting in 3,239,270 merged reads. The QDD pipeline v.3.1.2 (Meglécz et al. 2014) was used to detect sequences containing microsatellites (SSR mining), identify good quality sequences, and design candidate primer pairs flanking the identified SSR. From 991,696 sequences longer than 200 bp, a total of 54,149 sequences containing microsatellites were extracted by QDD pipe 1, from which 28,197 singleton sequences were retained for QDD pipe 2. Then, QDD pipe 3 designed primers in singleton sequences with the following filter parameters: amplicon size between 100 and 180 bp, primer size between 21 and 26 bp, melting temperature (Tm) of primers between 60 and 75°C, a maximum difference of Tm between primers of 10°C, and GC content between 40 and 60%. Several criteria were used to select 60 markers among the 10,372 candidates with designed primers: the selection of one primer set per sequence, the exclusion of motifs consisting exclusively of CG or AT, the retention of only perfect or compound microsatellites, ensuring minimal distance between primers and microsatellites, retaining loci with at least 10 repeats, and excluding loci whose primers contain multiple identical repeated bases. The 60 selected SSR primer pairs were tagged at the 5’-end with universal Illumina adapter overhang sequences: TCGTCGGCAGCGTCAGATGTGTATAAGAGACAG for forward primers, and GTCTCGTGGGCTCGGA GATGTGTATAAGAGACAG for reverse primers. These primers were ordered with standard desalted purification at Integrated DNA Technologies (Coralville, USA) and tested in simplex PCR. SSR markers that successfully amplified in simplex PCR were amplified in a multiplexed PCR following Lepais et al. (2020), but using a single-step multiplexed PCR before the indexing PCR, and sequenced on an Illumina iSeq 100 2 × 150 bp. The genotyping of 107 samples was blind-repeated to assess the genotyping error rate based on genotype reproducibility using a pipeline (Lepais et al. 2020) integrating FDSTools v.1.2.0 (Hoogenboom et al. 2016). We kept as validated loci those showing less than 20% of missing genotypes across samples and less than 3% of allelic genotyping error estimated based on repeated genotypes.

Size homoplasy was calculated for 25 high-quality validated loci as the number of alleles differing in sequence minus the number of alleles differing by length and divided by the number of alleles differing by their length. Data analyses were conducted per locus on the entire set of individuals (n = 107) and based on two different datasets, for which alleles were coded according to the amplicon length, and the sequence identity. The genetic variability of each dataset was characterized by locus by the number of alleles (N_A), the effective number of alleles (N_E), the allelic richness (R), the observed heterozygosity (H_O), and the expected heterozygosity (H_E). Rarefied allelic richness for a random subsample of gene copies (k = 166) was calculated based on the minimum sample size per locus (n = 83). The genetic diversity estimates were computed in SPAGeDi v.1.5a (Hardy and Vekemans 2002). Significant differences for N_A, N_E, R, H_O, and H_E between the amplicon length and the sequence identity datasets were tested using the paired Wilcoxon signed-rank test in R v.4.4.1 (R Core Team 2024).

Biologically informed statistical analyses were performed on population samples to determine the reliability of the SSRseq data for subsequent population genetic analyses. Genotyping errors and null alleles were assessed using Micro-Checker v.2.2.3 (Van Oosterhout et al. 2004). Hardy-Weinberg equilibrium and linkage disequilibrium were tested using the MCMC algorithm with the following parameters: 1000 iterations per batch, 100 batches, and 1000 dememorization steps, using Genepop v.4.7.5 (Raymond and Rousset 1995). The genetic diversity parameters (N_A, N_E, R, H_O, H_E) and the inbreeding coefficient F_IS were estimated and compared between the two life stages of the forest from the Paranaense region. To standardize comparisons, the rarefied allelic richness (R) was calculated using the minimum sample size (n = 26). Significant differences between the life stages were also tested as described above.

Forty-six developed loci out of 60 were successfully amplified in the simplex PCR test and 25 were validated as good-quality loci after amplification of the 46 loci in a single multiplex and sequencing (Table 1; Supplementary material 1). These 25 SSR loci produced high-quality genotypes for 107 A. colubrina samples (4.45% of missing genotypes and 0.8% of genotyping error).

The microsatellite loci were highly polymorphic in both sequences and length for most of them (Table 1). A high number of alleles per locus was detected in both datasets. Significantly more alleles were observed per locus based on sequence identity (16) than based on allele length (9) and expected heterozygosity was also higher in the sequence-based dataset (H_E = 0.744 vs H_E = 0.655) (Table 1). Size homoplasy reached a mean of 97.85% across loci, or, in other words, for each allele observed based on length, next to two (1.98) alleles were observed based on sequence. The number of alleles, allelic richness, and both observed and expected heterozygosity showed significant differences between allele calling methods (p = 4.645 × 10^-7 for N_A, p = 2.751 × 10^-5 for N_E, p = 2.019 × 10^-7 for R, p = 4.768 × 10^-5 for H_O, and p = 8.108 × 10^-6 for H_E) (Table 1).

The SSRseq adult trees dataset from the Paranaense forest showed no evidence of scoring errors due to stuttering or large allele dropout across all 25 loci in the genotyping error analysis. However, null alleles were detected at three loci (SSRseq.A06, SSRseq.A20, SSRseq.A43), potentially causing an excess of homozygosity. Consequently, these loci were excluded from subsequent population analyses. No statistically significant deviations from Hardy-Weinberg equilibrium were detected in none of the two life stages, despite expectations of higher relatedness due to family structure in the seedlings. No evidence of linkage disequilibrium was observed between pairs of loci.

High population genetic diversity was detected in adults and seedlings from the Paranaense forest (based on 22 loci). The mean and effective numbers of alleles per locus (N_A, N_E), observed and expected heterozygosity (H_O, H_E), and the inbreeding coefficient (F_IS) did not show significant differences between life stages. However, rarefied allelic richness (R) values were significantly higher in adults than in seedlings (p = 0.023) (Table 2).

High-throughput sequencing allowed the de novo development of an NGS-based multiplex marker panel for Anadenanthera colubrina. The effectiveness of SSRseq strengthens the advantage of using modern NGS platforms for qualitatively and quantitatively increasing data in large-scale population genetic studies.

The quality of DNA extraction is a crucial starting point for genetic studies. For Leguminosae forest tree species, the CTAB method is highly recommended because it is an inexpensive and rapid protocol that continues to be widely used in several species of the family (Riahi et al. 2010; Silva et al. 2023; Caycho et al. 2023; da Rocha et al. 2024).

Microsatellites have been used as highly polymorphic markers in previous population genetic studies of Anadenanthera colubrina, although the number of nuclear and plastid SSRs was limited (Barrandeguy et al. 2014; Goncalves et al. 2019; Feres et al. 2021). The de novo development of SSRs using traditional methods constitutes a time-consuming and costly endeavour (Edwards et al. 1996), particularly for non-model species such as A. colubrina, of which genomic resources are scarce. Conversely, the newer SSRseq technology is an efficient method for obtaining high-quality and informative species-specific nuclear microsatellite markers in non-model species. Additionally, multiplexing samples in a single sequencing run using individual-specific barcoding reduces the cost and time required for marker development and large-scale genotyping (Guo et al. 2020).

The 25 new SSRseq markers provide an additional advantage over traditional fragment-length genotyping. Here, the NGS-based method allowed us to detect a high percentage of size homoplasy and to resolve almost double the number of alleles than expected based on fragment length accessible by capillary electrophoresis-based SSR. The generated datasets increase the resolution for more accurate eco-evolutionary inference in A. colubrina populations. The percentage of increase in the detected number of alleles due to sequence analysis resulted as high as 97.85% in A. colubrina, while in other tree species, it was 36% in chestnut (Castanea sativa; C. crenata; Laurent et al. 2020), and 55% on oaks (Quercus faginea and Q. canariensis; Lepais et al. 2022).

Beyond detecting size homoplasy, the sequencing approach expands the scope of SSRs due to the use of compound marker systems integrating linked polymorphisms with different mutational dynamics, such as a microsatellite and its flanking sequences, allowing improvements in the estimation of population structure and inferences of demographic history (Payseur and Cutter 2006). This is especially relevant in the case of complex population dynamics (Lepais et al. 2022). Hence, the development and analysis of SSRseq provide an effective means to enhance our understanding of the processes affecting genetic diversity patterns, offering a valuable tool for informing conservation and management decisions.

High levels of population genetic diversity were detected in adults and seedlings from a Paranaense forest. The mean number of alleles per locus (N_A = 12.50) and the allelic richness after rarefaction (R = 10.49 in adults) were particularly high. These estimates exceed those reported in previous studies of A. colubrina at various spatial scales, which employed traditional nuclear microsatellites. For example, a survey of two populations located in northeastern São Paulo state, Brazil, detected a mean N_A = 7 (Feres 2013), while an analysis of four populations from the Yungas and Paranaense regions of Argentina found N_A = 9.72 and R = 6.75 (Barrandeguy et al. 2014). More recently, the genetic diversity of 79 populations representing four A. colubrina varieties from Brazil and northeastern Bolivia was characterized using 12 SSR loci, revealing a range of N_A = 3.36–5.59 (Mangaravite et al. 2023). Additionally, in other tree species from Seasonally Dry Tropical Forests, such as Tibouchina papyrus (R = 1.50–2.10, n = 66; Collevatti et al. 2012), and Myroxylon peruiferum (R = 2.73–3.01, n = 37; Silvestre et al. 2018), the estimated allelic richness was lower than that reported here.

Forest tree species generally exhibit a high genetic diversity, primarily distributed within populations, a pattern often attributed to life history traits such as longevity and outcrossing mating systems (Hamrick and Godt 1996; Petit and Hampe 2006). In our sample, adults had a higher allelic richness than seedlings (standardized for the same sampling size), which was expected considering that seedlings represented four half-sib families, i.e. 50% of their alleles corresponding to the maternal copies came from four mothers only, whereas adults represented a sample of naturally regenerated trees without a priori family structure. The overall high genetic diversity detected in this study may result from synergistic factors, including elevated multilocus outcrossing rates (Goncalves et al., manuscript in preparation), which suggest substantial pollen-mediated gene flow, and high mutation rates at SSR loci, further amplified by the polymorphisms identified within the SSRs and their flanking sequences.

The low inbreeding coefficients detected contrast with those reported in a previous study on fine-scale population genetic structure in A. colubrina, where both adults and saplings exhibited high F_IS coefficients (Goncalves et al. 2019). These discrepancies may arise from biological differences in the mating system, from differences in sampling design and/or differences in the nature of molecular marker variation. In the present study, distances between mother trees were kept above 50 m to minimize the likelihood of sampling closely related individuals, consistent with previous findings on a mixed mating system and relatively short-range seed dispersal in A. colubrina populations (Goncalves et al. 2019; Feres et al. 2021). Non-significant F_IS values in combination with high polymorphism in our study suggest high outcrossing rates and gene flow in Paranaense populations. The approximate selfing rate, s, estimated from the inbreeding coefficient using s = 2F_IS/(1 + F_IS) (Hartl and Clark 2007), suggests 4.5% selfing in this forest site. Therefore, the influence of the mating system and pollen dispersal on the genetic diversity of natural populations of A. colubrina in remnant forests can now be studied using SSRseq. This type of molecular marker may help foster research into non-model species, such as a recent study in other Neotropical forest tree species (Corvalán et al. 2023). This approach emphasizes the importance of understanding evolutionary and ecological processes and their impacts on ecosystem responses, thereby promoting the sustainability of South American fragmented native forests.

Recent advances in molecular and computational techniques, combined with transdisciplinary research that integrates ecology, evolution, and genetics, are crucial for understanding and conserving processes that support plant genetic diversity in a changing world. Overall, the development of SSRseq is a powerful tool for genetic analysis and can be used to identify genetic variation and diversity within and among populations, which can be useful for sustainable management and conservation policy.