Development of novel SSR molecular markers using a Next-Generation Sequencing approach (ddRADseq) in Stetsonia coryne (Cactaceae)

Abstract

The Cactaceae family is native to the American continent with several centers of diversity. In South America, one of these centers is the Central Andes and many species are considered to be threatened or vulnerable according to the International Union for Conservation of Nature (IUCN). Stetsonia coryne is an emblematic giant columnar cacti of the Chaco phytogeographic province. It has an extensive geographical distribution in many countries of the continent. However, to date there are no specific molecular markers for this species, neither reports of population genetic variability studies, such as for many cactus species. The lack of information is fundamentally due to the lack of molecular markers that allow these studies. In this work, by applying a Genotyping by Sequencing (GBS) technique, we developed polymorphic SSR markers for the Stetsonia coryne and evaluated their transferability to phylogenetically close species, in order to account for a robust panel of molecular markers for multispecies-studies within Cactaceae.

Key words
Cactaceae; Stetsonia coryne; GBS; ddRADSeq; Microsatellite

INTRODUCTION

The Cactaceae family, native to the American continent shows a high diversity of species in México and Southwestern United States, while in South America three main diversity centers have been recognized: 1). The Central Andes region, 2). East Brazil, 3). West and South Brazil, Paraguay, Uruguay (Ortega-Baes et al. 2010). Cacti are a group with special conservation interest as ornamental plants, in the food and medical industry, and as a source of wood (Casas & Barbera 2002, Inglese et al. 2002, Nerd et al. 2002, Nefzaoui & Salem 2002, Saenz-Hernandez et al. 2002, Wright et al. 2007, Ortega-Baes et al. 2010). Consequently, they have strong pressure for illegal collection, which along with the change in land use and climate are the most important threats for these species (Ortega-Baes et al. 2010). A recent global assessment of the family has indicated that approximately 31 % of its species present some degree of extinction risk (Goettsch et al. 2015) and many are considered to be threatened or vulnerable according to the International Union for Conservation of Nature (IUCN).

In the Central Andes, Northwest Argentina is one of the main centers of cactus diversity in the world (Ortega-Baes & Godínez-Alvarez 2006, Ortega-Baes et al. 2010). Studies of genetic structure in South America, however, are mainly concentrated in areas of Venezuela and Brazil (Moraes et al. 2005, Figueredo et al. 2010, Fava et al. 2016, Fernandes et al. 2016, Khan et al. 2018), some studies in central zone of Chile (Ossa et al. 2016, Larridon et al. 2018) and a single report assessed the genetic diversity of Echinopsis terscheckii, a species with Argentinean distribution (Quipildor et al. 2017, 2018). In other words, the study of the genetic variability of this family in Argentina is just beginning. This is mainly due to the absence of molecular markers that allow its study (Guerrero et al. 2019). For this reason, an exhaustive study of the genetic structure of different species of this family is necessary to evaluate the possible effect of these factors on the genetic variability of these species. This would allow the proposal of conservation actions and priority areas with the greater genetic variability of this family. Of the more than 2,000 cacti described, only around 28 wild species have been genetically analyzed to date (Cornejo-Romero et al. 2013).

Stetsonia coryne (Salm Dick Britton & Rose) is an emblematic giant columnar cactus (2n= 22, Sosa-Pivatto et al. 2014) appreciated for its use in ornamental, live fences, for its edible fruits and as a source of wood for handicrafts (Arenas & Scarpa 1998, Anderson 2001, Arenas & Scarpa 2007). This species its endemic of the Chaco phytogeographic province, which has an extensive geographical distribution that includes northern Argentina, southeastern Bolivia, and Paraguay. The Gran Chaco region is one of the areas with the highest deforestation rates globally (Gasparri & Grau 2009). Particularly, in the northwest of Argentina, this phytogeographic province has experienced the greatest expansion of the agricultural frontier over the last 25 years (Volante et al. 2016). Despite the wide distribution of the species, and its importance, to date there are no specific molecular markers for it, so there are no reports of studies of population genetic variability.

Microsatellites, or Single Sequence Repeats (SSRs), are repeated motifs of 1–6 nucleotides, highly distributed throughout the genomes of eukaryotes (Li et al. 2002). Even though SSR markers were firstly introduced in the late 1990s, they are still highly used for applications involving genome mapping, forensics, parentage identification, population studies, conservation genetics and phylogeography (Hodel et al. 2016). Indeed, SSR markers are the most commonly used molecular markers to study the genetic variability of natural plant populations, where sampling a high number of individuals, even with a small number of markers, is preferred over sampling a small number of individuals with a high number of markers (Bonatelli et al. 2015). The attractiveness of SSR markers relies on their abundance in the genome, high levels of polymorphism, high reproducibility, codominance and cost effectiveness (Li et al. 2002). Moreover, because SSR are PCR-based markers, they can be amplified using low concentration of DNA, even of low quality, and accordingly they are the markers of choice in studies involving ancient DNA and forensics (Zalapa et al. 2012, Hodel et al. 2016). Another important characteristic of these markers is the transferability between plant species of the same genus or family (Zalapa et al. 2012, Bombonato et al. 2019). Even though the success of amplification is inverse to the evolutionary distance between species, SSR allow multispecies studies (Barbará et al. 2007, Bombonato et al. 2019).

Until recently, the main drawback of these technique was the high cost associated with initial marker development, which required a large process of cloning and sequencing of the SSR loci (Guo et al. 2007, Davey et al. 2011). In this regard, SSR markers have been mostly developed for model organisms, or species of agronomic interest (Bonatelli et al. 2015), and then transferred to related species. Nowadays, this disadvantage has been overcome by the development of Next Generation Sequencing technologies (NGS), which facilitates the time and costs associated with the generation of genomic libraries (Davey et al. 2011). NGS derived genotyping strategies, as genotyping-by-sequencing (GBS, Elshire et al. 2011) and double digest Restriction Associated DNA sequencing (ddRADseq, Peterson et al. 2014) among others, have been widely used in recent years. In particular, ddRADseq is a strategy characterized by the use of two digestion enzymes, in which there is a reduction in the complexity of the genome, and it is possible to discover markers and genotype at the same time. Even though they were firstly developed for SNP marker identification, their use for SSR identification has been also reported (Qin et al. 2017, Aguirre et al. 2019) and allowed the study of the biodiversity of numerous plant species for which no specific markers were available (Deschamps et al. 2012, Egan et al. 2012, Zalapa et al. 2012, Bonatelli et al. 2015, Qin et al. 2017). Here it is important to note that, although both SSR and SNP markers share some favorable features, such as their codominant nature, high abundance and reproducibility, SSRs tend to be more variable and, consequently, more informative than SNPs. Given their multiallelic status, SSRs are still widely used in genetic diversity studies, especially those involving large sample sizes, with their associated high costs. Therefore, microsatellites continue to play an important role in population genetic studies, even in the current genomic era (Mason 2015, Hodel et al. 2016).

This work presents the first application of an NGS-derived strategy in the emblematic columnar cacti S. coryne. Using the ddRADseq approach, in this work we developed polymorphic SSR markers for the species and evaluated their transferability to phylogenetically close species, in order to account for a robust panel of molecular markers for multispecies-studies within Cactaceae.

MATERIALS AND METHODS

Plant material and DNA extraction

Ten individuals of Stetsonia coryne were collected in the northwestern of Argentina, from two populations, La Unión (Salta Province) and Dean Funes (Córdoba Province). The individuals were sampled from geographically distant locations to maximize the variability between them. Genomic DNA (gDNA) extraction was carried out from the parenchyma of the stem using a CTAB protocol (Bornet & Branchard 2001). The DNA quality was verified by Nanodrop (Thermo Fisher Scientific, Waltham MA, USA) and 1 % agarose gel electrophoresis analysis. The DNA samples were quantified in Qubit 2.0 fluorometer (Thermo Fisher Scientific) and the two samples with the highest DNA concentrations and best quality were selected for the ddRADseq protocol.

Additionally, the transferability of the generated SSRs was evaluated in three individuals of each one of four selected phylogenetically close species from the Browningeae-Cereae-Trichocereae clade (BCT: Cereus forbesii, C. aethiops, C. stenogonus, Echinpsis terscheckii), one from Notocacteae clade (Parodia microsperma) (Nyffeler 2002, Hernandez-Hernandez et al. 2011) and from two species of the Opuntiodea subfamily (Tunilla coruguata and T. soherensis). Isolation of gDNA was carried out as detailed before.

Evaluation and selection of restriction enzymes for ddRADseq

The first step for ddRADseq consisted of selecting digestion enzymes by performing both in silico and in vitro strategies. To carry out this selection process, the procedure proposed by Aguirre et al. (2019) was used, however, the enzymes combinations evaluated were different. The in silico strategy involved simulation of the digestion profiles expected for the enzyme combinations SphI-MboI and PstI-MboI with the SimRAD package (Lepais & Weir 2014). This package allows in silico predictions not only with genomes whose sequence are available, but with genomes with partial information (size in base pairs and percentage of GC content). That information is not available for S. coryne and thus partial genome information from the Carnegiea gigantea (2n=22, Copetti et al. 2017) available in https://www.ncbi.nlm.nih.gov/genome/39095, was used for the analysis. These simulations were carried out ten times to verify their validity. Additionality, regarding the in vitro strategy, DNA digestions were performed following the proposals of Aguirre et al. (2019). The fragmentation profile was visualized in agarose gel and capillary electrophoresis in a 5200 Fragment Analyzer System (Advanced Analytical Technologies, Inc., Santa Clara, CA, USA) using the DNA high-sensitivity kit (Agilent Technologies, Santa Clara, CA, USA).

Once the best enzyme pair was chosen, two selected samples were genotyped by ddRADseq, according to the protocol for small number of samples (Protocol 1) proposed by Aguirre et al. (2019). Briefly, this protocol has six steps. Digestion: two DNA samples were digested with the selected enzyme pair. Ligation: common adapters (double-stranded oligonucleotides), designed by Peterson et al. (2014), modified by adding sticky ends complementary for the recognition sites for the restriction enzyme pair were used. PCR: the dual-indexed primers designed by Lange et al. (2014) were used for PCR reactions and sample identification. These oligonucleotides have an index (8 bp), which allows sample identification. Pooling: Both samples were pooled in equimolar concentration after performing PCR. Size selection: the DNA fragments of interest (450-550 pb) were manually selected form the low-melting 1.5 % agarose gel. These fragments were purified from the gel using QIAquick Gel Extraction kit. Sequencing: the generated ddRADseq libraries were sequenced using MiSeq Paired-End (PE 2x250bp) reads (Illumina Inc., SanDiego, CA, USA) in the Genomic Unit of the Biotechnology Institute, IABIMO, INTA Castelar.

NGS data analysis, SSR identification and primer design

The processing of the NGS data was performed using Stacks v 1.42 (Catchen et al. 2013). Because there is no reference genome available, the “denovo_map.pl” function was implemented for the analysis. Raw Illumina reads were demultiplexed and quality checked using the process_radtags routine implemented in Stacks (v1.42, Catchen et al. 2013). All reads were trimmed to 200 bp after barcode removal. Sequence with an average Phred quality score lower than 30, or ambiguous restriction sites were eliminated from the analysis. Even though PE sequencing was performed, R1 and R2 reads were analyzed separately.

SSR markers were identified using the software MIcroSAtellite (MISA, Institute of plants Genetics and Crop Plant Research, Gatersleben, Germany, Beier et al. 2017). The parameters set for SSR calling were a minimum of five repeats for dinucleotide motifs, four repeats for trinucleotide motifs and three repeats for tetra, penta and hexanucleotide motifs. Even though all the repetitive sequences are reported here, only those loci that resulted polymorphic between or within samples were considered for primer development and posterior amplification.

Given that the NGS chemistry used herein was PE 2x250bp, when the repeated SSR motif was located approximately in the middle of the sequencing read, the same read was used for both primers design (i.e. forward and reverse), leaving an expected PCR fragment of 100-250 bp length. However, in such cases where the entire repeated SSR motif was located close to the 3’ end of the sequencing read, thus precluding the design of the pair of primers within the same read, the PE partner was included for primer design. We were not able to design the pair of primers only in those situations in which the repeated SSR motif was located 5’ of the sequenced read. PRIMER3 v4.0.0 (Untergasser et al. 2012) was used for PCR primer design. The repetitive region was set as the target region and the expected PCR product sizes were set from 100 to 250 bp. Minimum, optimum and maximum primer sizes were set as 18, 20, and 25, respectively; optimum and maximum melting temperature (Tm) were set at 58 °C, 63 °C, whereas minimum, optimum and maximum GC content were kept as default settings. The self-complementarity and hairpins for all primers were checked with Oligocalc (Kibbe 2007, http://www.basic. northwestern.edu/biotools/oligocalc.html).

SSR Validation, cross-amplification and diversity analysis

PCR amplification were performed in 20 μl final volume reactions, which contained 20 ng gDNA, 1 U Taq-DNA Polymerase (Inbio HighWay), 2 μl of buffer reaction TAS 10X, 2.5 mM MgCl₂, 0.5 μM of the reverse primer and forward primer, 0.4 mM of each dNTP, and Chromasolv® H₂O.

PCR reactions were performed with a touchdown program (Don et al. 1991). The protocol consisted of an initial denaturation step for 30 sec at 94 °C, followed by 10 cycles at 94 °C for 30 sec, 60 °C for 30 sec and 72 °C for 30 sec, followed by 30 cycles at 94 °C for 30 sec, 50 °C for 30 sec and a final extension step at 72 °C for 7 min. When necessary, annealing temperatures were adjusted for specific loci. Amplification products were visualized in 2 % agarose gels and in Fragment Analyzer. Amplification products close to the expected size were considered positive for both S. coryne and the other related species.

Additionally, basic diversity parameters were calculated for ten individuals belonging to selected populations (number of alleles (N_A), Major allele frequency (MAF), polymorphism information content (PIC), observed heterozygosity (H_o), expected heterozygosity (H_e) and Hardy-Weinberg Equilibrium (HWE). These analyses were performed by amplifying SSRs with a fluorescent dye-labeled forward primer, and then separating the amplified SSRs on Genetic Analyzer 3130xl (Applied Biosystems, Foster City, CA, USA). Allele assignments were made by size comparison with the standard allelic ladders, using the GeneMapper ID software provided by Applied Biosystems. The number of alleles, the expected heterozygosity (H_e), as well as the observed heterozygosity (H_o) were determined using Adegenet (Jombart 2008) in R environment. Hardy- Weinberg equilibrium (HWE) in each loci were tested using function hw.test from pegas package (Paradis 2010) in R, whereas the polymorphism information content (PIC) of each marker was calculated according to Botstein et al. (1980).

RESULTS

Evaluation and selection of the restriction enzyme pair for ddRADseq

In-silico simulations performed using the Carnegiea gigantea genome size (1400 Mbp; Copetti et al. 2017) and GC content information (36 %), as a proxy for enzyme restriction patterns in the non-model species S. coryne, showed notable differences between the enzyme pairs SphI-MboI and PstI-MboI (Figure 1, Supplementary Material Table SI). The first combination (SphI-MboI) is expected to generate 297,841.5 ± 712.78 type AB and BA fragments (i.e fragments digested with both enzymes, mean ± standard deviation (SD), 10 replicates), with 30,952.6 ± 109.48 of them of 300-400 bp (Figure 1b). On the other hand, the second enzyme pair combination (PstI-MboI) is expected to produce more fragments, of which 4,937,205.5 ± 3284.13 were type AB and BA, with 512,404.5 ± 537.85 of them of 300-400 bp (Figure 1a).

Figure 1
Histograms of in silico simulations (Numbers of fragments versus fragment size). (a) In silico digestion with PstI/MboI. (b) In silico digestion with Sph/MboI. The colored areas represent the portion of fragments of 300-400bp length. Blue areas specify only those fragments that also have the configuration AB+BA (i.e. were digested with different enzymes).

On turn, the in vitro digestion pattern was more homogeneous for the SphI-Mbol enzyme pair (Figure 2a), when comparing with PstI-MboI (Figure 2b). Taken the in silico and in vitro results together, and considering the observations of Aguirre et al. (2019), who suggested homogeneity and an optimum number of ~30,000 fragments as the most relevant features in enzyme selection, the pair SphI-MboI emerges as the most suitable combination, and thus used here for ddRADseq library preparation.

Figure 2
In vitro digestions profiles of S. coryne genomic DNA. Fragment Analyzer system runs. (a). Digestion with PstI/MboI and (b). Digestion with SphI/MboI.

NGS data processing, SSR identification and primer design

The high-quality sequencing reads that passed the quality filters corresponded to 98.36 % and 96.94 % of the reads for individual 1 and 2, respectively. However, sequencing depth was 10X higher in individual 1 than 2 (Table I). Sequence Read File (SRA) is available in PRJNA701953.

Table I
Generated reads with Miseq Illumina and demultiplexed reads before quality filters. R1: forward reads and R2: reverse reads.

Genotype	#generated reads	#demutiplexed reads
Individual 1 R1	2,037,765	2,012,346
Indiviudal 1 R2	2,037,765	1,996,607
Indiviudal 2 R1	176,810	172,847
Indiviudal 2 R2	176,810	169,968

The search for SSRs performed with MISA yielded 2,619 putative SSRs in the evaluated reads. The dominant repetitive motif was dinucleotide (53 %), followed by hexa (19 %), tetra (16 %), tri (7.4 %) and pentanucleotide (4.6 %) (Figure 3a). The dominant dinucleotide repetitive motif was TG (19.2 %) and the less abundant was GC (0.12 %). The most frequent trinucleotide, tetra, penta and hexanucleotide were AAT (11.01 %), TTTC (23.92 %), TTTCT (12.86 %), AAATGC (11%) respectively (Figure 3b, c, d, e, f).

Figure 3
(a) Number of total frequent motifs of SSRs. (b) Fequency of dinucleotids motifs (%) most abundant (c) trinucleotides (d)tetranucleotides, (e) pentanucleotides, and (f) hexanucleotides motifs found in S. coryne.

The analysis of only two individuals allowed us to detect 20 polymorphic loci (six between and 14 within samples), with the remaining 2,599 loci having the potential to be polymorphic in other individuals. Primers were successfully designed for 15 out of 20 regions. This means that 15 regions had sufficient flanking sequences to allow the design of appropriate unique primers (Table II).

Table II
Characteristics of the 15 polymorphic microsatellite loci developed for Stetsonia coryne.

LOCUS	Primer sequence	Primer size (bp)	Tm (ºC)	Repeat motif	Product size (bp)	GenBank code
STET1	F: GAATTTCCCCTTCACATCCA	20	59,73	(AG)11	111	MW600815
	R: TCTCATCTCAATCCCCATCC	20	59,82
STET2	F: TGCAAGAGAGTCATGATGGTG	20	59,85	(TC)11	105	MW600816
	R: GAGAGGGCAGGGACAAAAAG	21	61,12
STET3	F: CATGCTTGTTAAATTTTGCTGA	22	58,05	(TG)9	150	MW600817
	R: TACTTCCTTTAGTGTTCCTTACAATTC	27	57,73
STET4	F: CATGCCTCAAAGGGATATTTTT	22	59,36	(TTTC)5	150	MW600818
	R: TACCTTTGTATGTACTTGCACCTCTT	26	59,59
STET5	F: GATCCTCCGTTGACATCTTTACC	23	61,07	(TTC)13	150	MW600819
	R: GTTGTGCACTAAGCTAAAGAGAGG	24	58,85
STET6	F: CATGCCAATATGCCAATAAGG	21	60,18	(AG)17	150	MW600820
	R: GCTTCAAGGAGGATGAAAATG	21	58,78
STET7	F: GTTGCGTGGAATCTGTCTGA	20	59,84	(TC)14	192	MW600821
	R: GCAGTGCAGCACACTGAAAT	20	60,06
STET8	F: CAAGCAATTCGCGTGAAAAT	20	61,13	(TG)14	152	MW600822
	R: TCCTCAAATAATGCCACGAA	20	59,11
STET9	F: TGGTAGCATCACAAGCCAAG	20	59,86	(AAT)20	224	MW600823
	R: TTGATCTATGCAAATGTTCGATT	23	58,64
STET10	F: CATGCCTTTTGGCCTACCTA	20	60,09	(TG)10	140	MW600824
	R: CCCCCATAAAAAGAAAGCAA	20	59,05
STET11	F: TGAATCTTCTGCCGCTAACA	20	59,57	(CT)10	154	MW600825
	R: GTGCTTATTGCGACCTCCAT	20	60,1
STET12	F: ATCATTGTTCAGAGGGAGTG	20	55,01	(TTA)12	195	MW600826
	R: GCATTATGAAGGAGAGGACA	20	55,29
STET13	F: GCGTGCATGAGACAAAAACA	20	60,86	(AC)19	154	MW600827
	R: AAAGAGCAAAAGGGGAGAGG	20	59,82
STET14	F: TCAATTGATGAGCACAAACCA	21	60,1	(AG)13	192	MW600828
	R: TGGCAAATAGCAACTTGTATGA	22	58,38
STET15	F: TTGAATCATGCGTCATTAGTGT	22	58,17	(TG)9	184	MW600829
	R: CCCTCAGCCTTACTCCATGT	20	59,16

SSR validation, diversity analysis in S. coryne and evaluation of cross-amplification to close species

The 15 SSR loci derived from the GBS analysis were used for evaluating PCR amplification efficiency and polymorphism in ten individuals of S. coryne. Most of the primers amplified correctly (10 of the 15), whereas the rest showed nonspecific amplifications. One of them was monomorphic among the analyzed individuals.

The number of alleles per locus (N_A) ranged from 1 to 8 and the major allele frequency (MAF) varied from 0.273 to 1. The PIC value varied from 0.16 to 0.79. Ranges for observed (H_o) and expected heterozygosities (H_e) were 0.2- 0.64 and 0.18- 0.82, respectively, except for primer Stet6, which was monomorphic. Significant deviations from the HWE (p ≤ 0.05) were observed for only two of the ten markers evaluated for ten S. coryne individuals studied (Table III).

Table III
Microsatellite loci analyzed on 10 individuals from two evaluated populations of S. coryne.

Locus	N	NA	MAF	PIC	Ho	He	HWE p-value
Stet1	10	3	0.500	0.55	0.45	0.62	0.08
Stet2	10	6	0.682	0.49	0.36	0.51	0.018*
Stet3	10	3	0.850	0.25	0.3	0.27	1
Stet5	10	3	0.455	0.57	0.45	0.64	0.446
Stet6	10	1	1.000	0	0	0	1
Stet7	10	6	0.364	0.71	0.55	0.75	0.035*
Stet10	10	2	0.900	0.16	0.2	0.18	1
Stet11	10	6	0.364	0.72	0.64	0.76	0.051
Stet14	10	3	0.727	0.38	0.36	0.43	0.059
Stet15	10	8	0.273	0.79	0.64	0.82	0.037*

• N: Number of individuals, N_A: Number of alleles, MAF: minor allele frequency, PIC: Polymorphism information content, H_o: observed heterozygosity, H_e: expected heterozygosity and HWE: Hardy-Weinberg equilibrium, where: * is significant at p ≤ 0.05.

Ten of the 15 loci isolated from S. coryne were amplified in at least one of the other analyzed species, except for Opuntia species, which had no detectable amplifying loci. For the other species the amplification was successful to different degrees (Table IV). The Stet1, Stet2 and Stet7 loci were successfully amplified in P. microsperma, C. hankeanus, C. stenogonus, C. ahetiops and E. terscheckii, which were the phylogenetically species closest to S. coryne. Conversely, Stet8 and Stet12, Stet13, Stet14 and Stet15 failed to amplify in any studied species. The transferability success for P. microsperma was 60 %, whereas for the three species of the genus Cereus and of E. terscheckii, this parameter yielded values of 53.3 % and 46.6 %, respectively.

Table IV
Cross-amplification of 15 microsatellite markers isolated from S. coryne across seven species of Cactaceae.

Species	N	Stet1	Stet2	Stet3	Stet4	Stet5	Stet6	Stet7	Stet8	Stet9	Stet10	Stet11	Stet12	Stet13	Stet14	Stet15
Parodia microsperma	3	+	+	+	+	+	–	+	–	+	+	+	–	–	–	–
Cereus hankeanus	3	+	+	–	–	+	+	+	–	+	+	–	–	–	–	–
Cereus stenogonus	3	+	+	–	+	+	+	+	–	+	–	–	–	–	–	–
Cereus ahetiops	3	+	+	+	+	–	–	+	–	–	+	+	–	–	–	–
Echinopsis terscheckii	3	+	+	+	+	+	+	+	–	+	–	–	–	–	–	–
Opuntia tunilla	3	–	–	–	–	–	–	–	–	–	–	–	–	–	–	–
Opuntia sohorensis	3	–	–	–	–	–	–	–	–	–	–	–	–	–	–	–

DISCUSSION

Although cacti are an emblematic group in the arid and semi-arid regions of the American continent, to date the analysis of their genetic diversity has been limited due to the scarcity of available markers. Most of the studies carried out have been focused on species from the Northern Hemisphere (Cornejo-Romero et al. 2013, Fava et al. 2020), while the South American genetic studies are underrepresented (Larridon et al. 2018, Guerrero et al. 2019). In recent years, there has been an increase in the availability of markers and in population-genetic studies for species of this family in South America. However, these markers mainly belong to species present in different regions of Brazil (Fava et al. 2020). Genetic diversity studies and the design of molecular markers for species present in Argentina are still scarce.

NGS approaches have been used mainly in search of SNP polymorphic markers. This has allowed the generation of a large number of markers for both model and non-model species; however, when it comes to working with populations with large numbers of individuals, there is an economic impediment (Qin et al. 2017). For this reason, the application of these techniques in few individuals to generate other types of markers like SSR would be a great strategy, equally effective and more economical, that could be applied in poorly studied species that lack economic importance.

Today, there are a few reports of applications of NGS technology to non-model plant species, especially in the Cactaceae family, some of this in nuclear genome, (Bonatelli et al. 2015, Ossa et al. 2016, Fraga et al. 2020), and others in plastid genome (Chincoya et al. 2020, Hinojosa-Alvarez et al. 2020, Solórzano et al. 2019). A collateral application of NGS technique in plants is the cost-effective discovery of other genetic variants like polymorphic SSR loci, which still are the most commonly used markers for addressing genetic diversity in research with similar results to those of SNP markers (Bonatelli et al. 2015, Qin et al. 2017). The use of these technologies allowed us to discover markers for non-model species that still lacked available markers, or a reference genomes (Bonatelli et al. 2015).

A critical step to start the GBS techniques is the selection of the restriction enzymes for the analyses. In the present work, we evaluated two enzyme selection strategies. Both the in silico analysis and the in vitro analysis yielded different results. The first analysis showed that the PstI-MboI pair enzyme generated more fragments, although with the less homogeneous digestion profile. According to Aguirre et al. (2019) homogeneity is a desirable feature and this pair failed to meet this characteristic. The in vitro analysis showed the SphI-MboI enzyme pair generated the best profile. Possibly, these differences between simulated and real digestion profiles are due to the phylogenetical distance of the species and the information of the reference genome (that of Carnegiea gigantea) that was used, even though this species belongs to the same family. Although the size of the S. coryne genome is not available, the use of sim.DNAseq module of the SimRAD package allowed us to simulate the digestion profile of the genome by using information available from the genome of another species within the same family. Probably, the use of this tool would allow an early selection of potentially useful enzymes when de genome used belongs to a closer species.

In the present work, by applying the ddRADseq technique in only two S. coryne individuals, a total of 2619 SSRs were found. As reported for species of other families (Brassicaceae and Arecaceae), motifs rich in A/T was more dominant, while motifs rich in C /G was absent or scarce (Shi et al. 2014, Xiao et al. 2016). We were able to develop 15 SSRs, 10 of which were successfully amplified. The number of alleles found varied between 2 and 8 (except for one of the loci that was monomorphic among the samples analyzed). These values were similar to those reported by other works that use NGS techniques to search for SSR markers in the Cactaceae family. These values in those studies ranged from 2 to 10 (Bonatelli et al. 2015, Ossa et al. 2016, Fraga et al. 2020). Furthermore, these values were superior to those found in works using conventional SSR search techniques. In those cases, the number of alleles found varied between 2 and 7 (Hardesty et al. 2008, Fernandes et al. 2016). The large number of alleles found in the 10 samples analyzed, which is indicative of the high genetic variability in the S. coryne populations, seems to be consistent with their cross-breeding system.

One of the main characteristics of SSR markers is transferability, that is, they can be transferred between phylogenetically close species (Yan et al. 2017). In this study, we used NGS technics to develop markers specific to the assessed specie. Eight of these markers could be transferred to some very close species of this family. The transferability values were higher than those reported by Kalia et al. (2020) for the Bignoneaceae family (40.58%) and by Aiello et al. (2020) between two species of the family of the Apiaceae family (23%). Considering the recent divergence among cacti (Arakaki et al. 2011), the high percentage of transferability observed among the studied species (46.6-60 %) may be due to the degree of phylogenetic proximity between them. The obtained values are higher than those found by Bombonato et al. (2019), who determined 35 % transferability between species. However, the markers evaluated by them belong to more distant clade species. In this study, transferability was analyzed between species belonging to a more robust clade (core Cactoideae II), as proposed by Hernandez-Hernandez (2011). This may explain the highest percentage of transferability. Likewise, the enzymes used are sensitive to methylation, that is, we are sampling a greater number of coding regions, which are usually more conserved. This would also explain the higher transfer rate of the markers developed here. In contrast to this finding, as expected, the species of the genus Opuntia lacked transferability. Indeed, these species belong to a different subfamily and therefore the phylogenetic distance between them is greater.

The study of the variability and the genetic structure in cacti is still scarce. For this reason, the availability of markers that can be transferred to other species in the family is very low. The loci presented here are a valuable resource for the development of future research not only in the present species, but in other species of the Cactaceae family without available markers. Because of the lack of genetic-population studies conducted in South American cactus species (Cornejo-Romero et al. 2013, Fava et al. 2020), the markers developed in this study constitute a significant contribution in the advancement of the knowledge of the genetic variability existing in this plant family, especially South American cacti species and lineages not studied at present.

ACKNOWLEDGMENTS

The authors thank J. Sabio y Garcia for assistance with the English version of the manuscript, M. Galvan and D. Tosto for technical advice and M. Gonzalez, C. Mamani and P. Gorostiague for field assistance. The study was funded by FONCYT (Grant No. 1492/11), CIUNSa (Grant No.2355/0) and CONICET scholarship to A. Gutierrez.

REFERENCES

• AGUIRRE NC ET AL. 2019. Optimizing DDRADseq in non-model species: A Case Study in Eucalyptus dunnii Maiden. Agronomy 9: 1-21.
• ANDERSON EF. 2001. The Cactus Family. in press, T. (Ed.) Timber Press, Inc., Portland, Oregon.
• ARAKAKI M, CHRISTIN P-A, NYFFELER R, LENDEL A, EGGLI U, OGBURN RM, SPRIGGS E, MOORE MJ & EDWARDS EJ. 2011. Contemporaneous and recent radiations of the world’s major succulent plant lineages. Proc Natl Acad Sci 108: 8379-8384.
• ARENAS P & SCARPA GF. 1998. Ethnobotany of Stetsonia coryne (Cactaceae), the ‘Cardón’ of the Gran Chaco. Haseltonia 6: 42-51.
• ARENAS P & SCARPA GF. 2007. Edible wild plants of the Chorote Indians, Gran Chaco, Argentina. Bot J Linn Soc 153: 73-85.
• BARBARA T, PALMA-SILVA C, PAGGI GM, BERED F, FAY MF & LEXER C. 2007. Cross-species transfer of nuclear microsatellite markers: potential and limitations. Mol Ecol 16(18): 3759-3767.
• BEIER S, THIEL T, MÜNCH T, SCHOLZ U & MASCHER M. 2017. MISA-web: a web server for microsatellite prediction. Bioinformatics 33: 2583-2585.
• BOMBONATO JR, BONATELLI IAS, SILVA GAR, MORAES EM, ZAPPI DC, TAYLOR NP & FRANCO FF. 2019. Cross-genera SSR transferability in cacti revealed by a case study using Cereus (Cereeae, Cactaceae). Genet Mol Biol 42: 87-94.
• BONATELLI IAS, CARSTENS BC & MORAES EM. 2015. Using next generation RAD sequencing to isolate multispecies microsatellites for Pilosocereus (Cactaceae). PLoS ONE 10: e0142602.
• BORNET B & BRANCHARD M. 2001. Nonanchored Inter Simple Sequence Repeat (ISSR) Markers: Reproducible and Specific Tools for Genome Fingerprinting. Plant Mol Biol Report 19: 209-215.
• BOTSTEIN D, WHITE RL, SKOLNICK M & DAVIS4 RW. 1980. Construction of a Genetic Linkage Map in Man Using Restriction Fragment Length Polymorphisms. Am J Hum Genet 32: 314-331.
• CASAS A & BARBERA G. 2002. “Mesoamerican Domestication and Diffusion” In: Cacti. Biology and uses. Nobel P.S. (Ed) UNIVERSITY OF CALIFORNIA PRESS, Berkeley Los Angeles London.
• CATCHEN J, HOHENLOHE PA, BASSHAM S, AMORES A & CRESKO WA. 2013. Stacks: An analysis tool set for population genomics. Mol Ecol 22: 3124-3140.
• CHINCOYA DA, SANCHEZ-FLORES A, ESTRADA K, DÍAZ-VELÁSQUEZ CE, GONZÁLEZ-RODRÍGUEZ A, VACA-PANIAGUA F, DÁVILA P, ARIAS S & SOLÓRZANO S. 2020. Identification of High Molecular Variation Loci in Complete Chloroplast Genomes of Mammillaria (Cactaceae, Caryophyllales). Genes 11(7): 830.
• COPETTI D ET AL. 2017. Extensive gene tree discordance and hemiplasy shaped the genomes of North American columnar cacti. Proc Natl Acad Sci U. S. A. 114: 12003-12008.
• CORNEJO-ROMERO A, VARGAS-MENDOZA CARLOS F, VALVERDE P & RENDÓN-AGUILAR B. 2013. Estructura genética y filogeografía en cactáceas. Cactáceas y suculentas Mex 58: 4-28.
• DAVEY JW, HOHENLOHE PA, ETTER PD, BOONE JQ, CATCHEN JM & BLAXTER ML. 2011. Genome-wide genetic marker discovery and genotyping using next-generation sequencing. Nat Rev Genet 12: 499-510.
• DESCHAMPS S, LLACA V & MAY GD. 2012. Genotyping-by-sequencing in plants. Biology (Basel) 1: 460-483.
• DON RH, COX PT, WAINWRIGHT BJ, BAKER K & MATTICK JS. 1991. ‘Touchdown’ PCR to circumvent spurious priming during gene amplification. Nucleic Acids Res 19: 4008.
• EGAN AN, SCHLUETER J & SPOONER DM. 2012. Applications of next-generation sequencing in plant biology. Am J Bot 99: 175-185.
• ELSHIRE RJ, GLAUBITZ JC, SUN Q, POLAND JA, KAWAMOTO K, BUCKLER ES & MITCHELL SE. 2011. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE 6: e19379.
• FAVA WS, GOMES VGN, LORENZ AP & PAGGI GM. 2020. Cross-amplification of microsatellite loci in the cacti species from Brazilian Chaco. Mol Biol Rep 47: 1535-1542.
• FAVA WS, PAGGI GM, ZANELLA CM & LORENZ-LEMKE AP. 2016. Development and characterization of microsatellite markers for Echinopsis rhodotricha and cross-amplification in other species of Cactaceae. Biochem Syst Ecol 66: 19-23.
• FERNANDES VN DE A, DAS NEVES AF, MARTIN PG, MANGOLIN CA & MACHADO M DE FPS. 2016. Genetic structure and molecular divergence among samples of mandacaru (Cereus spp.; Cactaceae) as revealed by microsatellite markers. Biochem Syst Ecol 64: 38-45.
• FIGUEREDO CJ, NASSAR JM, GARCÍA-RIVAS AE & GONZÁLEZ-CARCACÍA JA. 2010. Population genetic diversity and structure of Pilosocereus tillianus (Cactaceae, Cereeae), a columnar cactus endemic to the Venezuelan Andes. J Arid Environ 74: 1392-1398.
• FRAGA DA, DE CARVALHO AF, SANTANA RS, MACHADO MC & LACORTE GA. 2020. Development of microsatellite markers for the threatened species Coleocephalocereus purpureus (Cactaceae) using next-generation sequencing. Mol Biol Rep 47: 1485-1489.
• GASPARRI NI & GRAU HR. 2009. Deforestation and fragmentation of Chaco dry forest in NW Argentina (1972-2007). For Ecol Manage 258: 913-921.
• GOETTSCH B ET AL. High proportion of cactus species threatened with extinction. Nat Plants 1: 15142.
• GUERRERO PC, MAJURE LC, CORNEJO-ROMERO A & HERNÁNDEZ-HERNÁNDEZ T. 2019. Phylogenetic Relationships and Evolutionary Trends in the Cactus Family. J Hered 110: 4-21.
• GUO W ET AL. 2007. A microsatellite-based, gene-rich linkage map reveals genome structure, function and evolution in Gossypium. Genetics 176: 527-541.
• HARDESTY BD, HUGHES SL, RODRIGUEZ VM & HAWKINS JA. 2008. Characterization of microsatellite loci for the endangered cactus Echinocactus grusonii, and their cross-species utilization. Mol Ecol Resour 8: 164-167.
• HERNANDEZ-HERNANDEZ T, HERNANDEZ HM, DE-NOVA JA, PUENTE R, EGUIARTE LE & MAGALLON S. 2011. Phylogenetic relationships and evolution of growth form in Cactaceae (Caryophyllales, Eudicotyledoneae). Am J Bot 98: 44-61.
• HINOJOSA-ALVAREZ S, ARIAS S, FERRAND S, PURUGGANAN MD, ROZAS J, ROSAS U & WEGIER A. 2020. The chloroplast genome of the pincushion cactus Mammilllaria haageana subsp. san-a ngelensis, a Mexican endangered species. Mitochondrial DNA Part B 5(3): 2038-2039.
• HODEL RGJ ET AL. 2016. A New Resource for the Development of SSR Markers: Millions of Loci from a Thousand Plant Transcriptomes. Appl Plant Sci 4: 1-6.
• INGLESE P, BASILE F & SCHIRRA M. 2002. “Cactus Pear Fruit Production”. In: Cacti. Biology and uses. Nobel PS (Ed) UNIVERSITY OF CALIFORNIA PRESS, Berkeley Los Angeles London.
• JOMBART T. 2008. “adegenet: a R package for the multivariate analysis of genetic markers.” Bioinformatics 24: 1403-1405.
• KALIA RK, CHHAJER S & PATHAK R. 2020. Cross genera transferability of microsatellite markers from other members of family Bignoniaceae to Tecomella undulata (Sm.) Seem Acta Physiologiae Plantarum 42(9): 1-14.
• KHAN G, GODOY MO, FRANCO FF, PEREZ MF, TAYLOR NP, ZAPPI DC, MACHADO MC & MORAES EM. 2018. Extreme population subdivision or cryptic speciation in the cactus Pilosocereus jauruensis? A taxonomic challenge posed by a naturally fragmented system. Syst Biodivers 16: 188-199.
• KIBBE WA. 2007. OligoCalc: An online oligonucleotide properties calculator. Nucleic Acids Res 35: 43-46.
• LANGE V ET AL. 2014. Cost-efficient high-throughput HLA typing by MiSeq amplicon sequencing. BMC Genomics 15: 1471-2164.
• LARRIDON I, VELTJEN E, SEMMOURI I, ASSELMAN P, GUERRERO PC, DUARTE M, WALTER HE, CISTERNAS MA & SAMAIN MS. 2018. Investigating taxon boundaries and extinction risk in endemic Chilean cacti (Copiapoa subsection Cinerei, Cactaceae) using chloroplast DNA sequences, microsatellite data and 3D mapping. Kew Bull 73: 55.
• LEPAIS O & WEIR JT. 2014. SimRAD: An R package for simulation-based prediction of the number of loci expected in RADseq and similar genotyping by sequencing approaches. Mol Ecol Resour 14: 1314-1321.
• LI YC, KOROL AB, FAHIMA T, BEILES A & NEVO E. 2002. Microsatellites: Genomic distribution, putative functions and mutational mechanisms: A review. Mol Ecol 11: 2453-2465.
• MASON AS. 2015. SSR Genotyping. In Plant Genotyping: Methods and Protocols, Methods in Molecular Biology. 1245: 77-89.
• MORAES EM, ABREU AG, ANDRADE SCS, SENE FM & SOLFERINI VN. 2005. Population genetic structure of two columnar cacti with a patchy distribution in eastern Brazil. Genetica 125: 311-323.
• NEFZAOUI A & SALEM H B. 2002. “ Forage, Fodder, and Animal Nutrition”. In: Cacti. Biology and uses. Nobel PS (Ed) UNIVERSITY OF CALIFORNIA PRESS, Berkeley Los Angeles London.
• NERD A, TEL-ZUR N & MIZRAHI Y. 2002. “Fruits of Vine and Columnar Cacti”. In: Cacti. Biology and uses. Nobel PS (Ed) UNIVERSITY OF CALIFORNIA PRESS, Berkeley Los Angeles London.
• NYFFELER R. 2002. Phylogenetic relationships in the cactus family (Cactaceae) based on evidence from trnK/matK and trnL-trnF sequences. Am J Bot 89(2): 312-326.
• ORTEGA-BAES P & GODÍNEZ-ALVAREZ H. 2006. Global diversity and conservation priorities in the cactaceae. Biodivers Conserv 15: 817-827.
• ORTEGA-BAES P, SÜHRING S, SAJAMA J, SOTOLA E, ALONSO-PEDANO M, BRAVO S & GODÍNEZ-ALVAREZ H. 2010. Diversity and Conservation in the Cactus Family. In Desert Plants: Biol Biotechnol, 1-503 p.
• OSSA CG, LARRIDON I, PERALTA G, ASSELMAN P & PÉREZ F. 2016. Development of microsatellite markers using next-generation sequencing for the columnar cactus Echinopsis chiloensis (Cactaceae). Mol Biol Rep 43: 1315-1320.
• PARADIS E. 2010. pegas: an R package for population genetics with an integrated–modular approach.Bioinformatics 26(3): 419-420.
• PEAKALL R & SMOUSE PE. 2006. GENALEX 6: Genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes 6: 288-295.
• PETERSON BK, WEBER JN, KAY EH, FISHER HS & HOEKSTRA HE. 2012. Double digest RADseq: An inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS ONE 7(5): e3713.
• PETERSON GW, DONG Y, HORBACH C & FU YB. 2014. Genotyping-by-sequencing for plant genetic diversity analysis: A lab guide for SNP genotyping. Diversity 6: 665-680.
• QIN H, YANG G, PROVAN J, LIU J & GAO L. 2017. Using MiddRAD-seq data to develop polymorphic microsatellite markers for an endangered yew species. Plant Divers 39: 294-299.
• QUIPILDOR VB, KITZBERGER T, ORTEGA-BAES P, QUIROGA MP & PREMOLI AC. 2018. Regional climate oscillations and local topography shape genetic polymorphisms and distribution of the giant columnar cactus Echinopsis terscheckii in drylands of the tropical Andes. J Biogeogr 45: 116-126.
• QUIPILDOR VB, MATHIASEN P & PREMOLI AC. 2017. Population Genetic Structure of the Giant Cactus Echinopsis terscheckii in Northwestern Argentina Is Shaped by Patterns of Vegetation Cover. J Hered 108: 469-478.
• SÁENZ-HERNANDEZ C, CORRALES-GARCIA JOEL & AQUINO-PÉREZ G. 2002. “Nopalitos, Mucilage, Fiber, and Cochineal”. 2002. In: Cacti. Biology and uses. Nobel PS (Ed) UNIVERSITY OF CALIFORNIA PRESS, Berkeley Los Angeles London.
• SHI J, HUANG S, ZHAN J, YU J, WANG X, HUA W, LIU S, LIU G & WANG H. 2014. Genome-wide microsatellite characterization and marker development in the sequenced brassica crop species. DNA Res 21: 53-68.
• SOLÓRZANO S, CHINCOYA DA, SANCHEZ-FLORES A, ESTRADA K, DÍAZ-VELÁSQUEZ CE, GONZÁLEZ-RODRÍGUEZ A, VACA-PANIAGUA F, DÁVILA P & ARIAS S. 2019. De novo assembly discovered novel structures in genome of plastids and revealed divergent inverted repeats in Mammillaria (Cactaceae, Caryophyllales). Plants 8(10): 392.
• SOSA-PIVATTO MS, FUNES G, FERRERAS AE & GURVICH DE. 2014. Seed mass, germination and seedling traits for some central Argentinian cacti. Seed Sci Res 24: 71-77.
• UNTERGASSER A, CUTCUTACHE I, KORESSAAR T, YE J, FAIRCLOTH BC, REMM M & ROZEN SG. 2012. Primer3-new capabilities and interfaces. Nucleic Acids Res 40: 1-12.
• VOLANTE JN, MOSCIARO MJ, GAVIER-PIZARRO GI & PARUELO JM. 2016. Agricultural expansion in the Semiarid Chaco: Poorly selective contagious advance. Land Use Policy 55: 154-165.
• WRIGHT CI, VAN-BUREN L, KRONER CI & KONING MMG. 2007. Herbal medicines as diuretics: A review of the scientific evidence. J Ethnopharmacol 114: 1-31.
• XIAO Y, XIA W, MA J, MASON AS, FAN H, SHI P, LEI X, MA Z & PENG M. 2016. Genome-wide identification and transferability of microsatellite markers between palmae species. Front Plant Sci 7: 1578.
• YAN Z, WU F, LUO K, ZHAO Y, YAN Q, ZHANG Y, WANG Y & ZHANG J. 2017. Cross-species transferability of EST-SSR markers developed from the transcriptome of Melilotus and their application to population genetics research. Scientific Reports 7(1): 1-11.
• ZALAPA JE, CUEVAS H, ZHU H, STEFFAN S, SENALIK D, ZELDIN E, MCCOWN B, HARBUT R & SIMON P. 2012. Using next-generation sequencing approaches to isolate simple sequence repeat (SSR) loci in the plant sciences. Am J Bot 99: 193-208.

SUPPLEMENTARY MATERIAL

Table SI.

Publication Dates

• Publication in this collection
  30 Aug 2021
• Date of issue
  2021

History

• Received
  13 Nov 2020
• Accepted
  8 Apr 2021