RAD sequencing for the development of microsatellite markers for identification of Malaysian taro cultivars

Abstract

The progress in next-generation sequencing (NGS) has transformed the discovery and development of microsatellite markers, particularly for understudied species. In this study, one of the NGS platforms adopted was restriction site-associated DNA sequencing (RAD-seq) to identify and develop microsatellite markers for taro. A total of 8,490,604 clean reads were generated and assembled into 222,445 contigs with only 3,280 contigs found to contain a total of 3454 repeats. From this total, only 1,790 repeats were successfully designed as a flanking primer, and 14 microsatellites were selected at random for cultivar identification of Malaysian taro. The analysis of 14 novel simple sequence repeat (SSR) markers only showed that twelve SSR markers were polymorphic with the remaining two SSR markers, CesSSR9 and CesSSR11, monomorphic. The analysis of polymorphic SSR markers revealed the presence of 50 alleles ranging between two (CesSSR5) and seven (CesSSR13) having an average of 4.17 alleles per locus. The values of the polymorphic information content ranged from 0.05 (CesSSR5) to 0.59 (CesSSR7) with an average of 0.33 per locus. The unweighted pair group method with arithmetic mean (UPGMA) dendrogram showed that cultivar MINYAK was unable to be differentiated from cultivar PUTIH. As such, these newly developed microsatellite markers will provide a tool for molecular identification in taro and provide necessary information for future biodiversity conservation and gene bank management of taro.

Keywords:

• NGS sequencing
• RAD sequencing
• SSR marker
• Taro

Introduction

Taro is considered a root crop that is vegetatively propagated and grown in most tropical regions globally. It belongs to the monocotyledonous family of Araceae. Its flowers can blossom and set seed and it is considered an important socio-economic crop in the Pacific and Southeast Asia regions. Many studies have attempted to discover the origins of this crop, with New Guinea and northeast India being revealed as possible locations for the initial domestication of this crop [1]. Moreover, the enormous variety of the wild Colocasia species of the crop appears to extend from the north-eastern part of India to the Himalayas in southern China [2]. Even though an abundance of genetic resources exists, the numerous attempts to exploit and conserve the germplasm needed to sustain the production of this species remain unresolved. Therefore, precise and accurate molecular characterisation utilising molecular markers is needed for appropriate conservation management and future breeding programme for taro. Microsatellites, also known as simple sequence repeats (SSR), are one of the most broadly used DNA markers that has been employed over the last few decades in cultivar identification [3]. These markers are reproducible, highly polymorphic, codominant, and abundant throughout the eukaryotic genome [4,5], making them favourable DNA markers for a variety of genetic research studies.

Despite the more recent and remarkable marker system, namely single nucleotide polymorphism (SNP), microsatellites are still preferred in conducting various research including cultivar identification and authentication, population structure analyses, genetic drift and gene flow and pedigree or breeding estimation [6]. Likewise, the advancement in next-generation sequencing (NGS) has evolved into a mainstream approach in dealing with high-throughput in identifying and discovering a thousand microsatellites at reduced costs, even in non-model species [7,8].

Using the conventional approach, cloning and unrelenting Sanger sequencing were needed in order to discover even a small number of microsatellites. As such, the NGS approach is preferred. Also, one of the valuable NGS based approaches is RAD-seq, which targets the genomic sequences around the restriction sites before being assembled to generate contigs for microsatellite screening. The RAD-seq approach can produce efficient and effective results for multiple applications, namely, population differentiation, marker improvement, and phylogeography [9–11]. RAD-seq is also more advantageous compared to the more conventional hybridisation approach since the technique allows to sequence the same region across individual or diverse species.

Materials and methods

Plant material and DNA extraction for RAD sequencing

A sample of C. esculenta (cultivar WANGI) was obtained from the MARDI Taro Collection at MyGenebank MARDI. Following the manufacturer’s protocol, total genomic DNA was extracted and collected from fresh young leaves using the Qiagen Plant DNA extraction kit. The DNA integrity and concentration were next measured using 0.8% agarose gel and Fluoraskan Asken, respectively. The genomic DNA was then sent to MyTACG Bioscience Enterprise (Sri Petaling, Malaysia) for RAD sequencing.

RAD-seq sequencing and contigs assembly

Following the standard manufacturer’s protocol, MyTACG Bioscience Enterprise created a RAD library where total genomic DNA was digested with MseI before sequencing using an Illumina HiSeq 2000 sequencer. This was followed by the processing of raw data generated using Stacks [12], which was later used for de novo assembly using Velvet [13].

SSR identification and primer design

MISA, a MIcro-SAtellite identification tool, was used to screen the assembled contigs for microsatellite repeat motifs [14] in which mononucleotide repeats were omitted from the search criteria. The search criteria consisted of the following series of parameters 6, 5, 4, 3, with 3 repeats for the di-, tri-, tetra-, penta- and hexa-repeat types, respectively. The primer pairs were designed using the Primer3Plus software programme with the criterion of the amplicon size between 100 and 250 bp. All designed primer pairs were optimised before obtaining the best annealing temperature of the amplification process.

Molecular characterisation using newly developed microsatellite markers

A total of 14 novel SSR markers (Table 1) were chosen to characterise the thirteen individuals representing 9 different cultivars. The polymerase chain reaction (PCR) was carried out as proposed by Schuelke [15] by ligating the primers (either forward or reverse) and fluorescent dye (FAM, PET, NED or VIC) with a non-fluorescent labelled M13 sequence tail [M13 sequence: TGTAAAACGACGGCCAGT]. The mixture of the PCR reaction was prepared in order to achieve a 10 μL target volume, consisting of 1x buffer (Invitrogen, United States), 10 μmol/L of each forward and reverse primer, 5 μmol/L fluorescence-labelled M13 adaptor, 2 μmol/L of each deoxyribonucleside triphosphate (dNTP) (Invitrogen, United States), 0.1 μL of bovine serum albumin (BSA) as PCR enhancer and 1 U of Taq polymerase (Invitrogen, United States). Amplification was performed utilising a GeneAmp® PCR System 9700 (Applied Biosystems, United States). The PCR profile consisted of pre-denaturation for 2 min at 94 °C, followed by 35 cycles at 94 °C for 30 s, 41–65 °C for 45 s, and 72 °C for 45 s, and followed lastly by post-extension at 72 °C for 5 min.

Table 1. List of SSR markers used to characterise taro cultivars.

SSR ID	Primer sequence	Acc. number*
CesSSR1	F: TGAGGAGTTGTTGTGCCTC	LC484345
	R: AATAACGGCGCGAGTTTGG
CesSSR2	F: GAAATGATGTGAAGACAGCGG	LC484346
	R: ATTCCAGCCCAGGTGACAG
CesSSR3	F: GGGAGAAGGGCTCGGAAG	LC484347
	R: TACCCAACCCAAATACGGC
CesSSR4	F: AACGACTATAAGGAGATTTGATAGC	LC484348
	R: TGAGAAGGGTCAGCAGTTC
CesSSR5	F: CCCAGAAACTCTCAGCGAAG	LC484349
	R: AGCCCAATCAAGAGAACCC
CesSSR6	F: TGGGCTAAGTTCTTGGGTTAG	LC484350
	R: CCACATGTCTCGGTTGCAG
CesSSR7	F: ACGGGTCTTACTGTGTGGAC	LC484351
	R: TTCCTACGGCCGCTTTCTG
CesSSR8	F: CCACAATATTTGTTATGGTTAGTCCC	LC484352
	R: CGCTCTGACTTTGGGCTATG
CesSSR9	F: CGCGATCGAGCCATAGCC	LC484353
	R: TGTACTAGTGTCTGGTTACTATTGAG
CesSSR10	F: TCTTGAACGCGAAACATGAAAG	LC484354
	R: GCTCCTAAGTCCTAACCGTACC
CesSSR11	F: AGAGATTAGCAAGTTTGGATATGG	LC484355
	R: TCTAACTCAGATCTCTAATGCCAC
CesSSR12	F: GCTCTAGGTGACAAAGGATGC	LC484356
	R: ACCCTGCACAATGTATTTCAGC
CesSSR13	F: CTGCCTCAGTCATGCCAAC	LC484357
	R: TCTCACCAAAGGTATACTAAAGAAGG
CesSSR14	F: CCTTCATTTCACATAGCCTCTCAAG	LC484358
	R: TATAGCTCGAGGACACGGG

* GenBank accession number, NCBI/DDBJ.

Following amplification, the PCR product was then multiplexed up to four primers with a different combination of fluorescent dyes. This was followed by mixing the products with Hi-Di formamide and GeneScan 500 LIZ (standard molecular weight ladder) prior to being resolved using an ABI 3130xL Genetic Analyser (Applied Biosystems, United States). GeneMapper Version 5 (Thermo Fisher Scientific, United States) was used to score the size of the allele. The produced electropherograms were accordingly scored and analysed, as described by Arif et al. [16].

To calculate the number of alleles, the Major allele Frequency, Gene Diversity, the extent of heterozygosity, and Polymorphism Information Content (PIC) of each microsatellite marker, the PowerMarker software was used. The same software was used to calculate the shared allele based genetic-distance [17]. To visualise the dendrogram, MEGA7 was used founded on the Unweighted Pair Group Method with Arithmetic Mean (UPGMA), constructed by applying the produced matrices genetic distance of the shared microsatellite alleles [18].

Results

RAD-seq sequencing and microsatellite analysis

In this study, a total of 8.49 million clean reads were produced. Using Velvet, these clean reads were then assembled into 222,445 contigs. A total of 3,280 contigs were found to contain a total of 3,454 repeats (Table 2). From the 3,454 repeats, 2,760 (79.9%) were dinucleotide repeat motifs, 631 (18.3%) trinucleotide repeat motifs, 40 (1.2%) tetranucleotide repeat motifs, 15 (0.4%) pentanucleotide repeat motifs and 8 (0.2%) hexanucleotide repeat motifs (Table 3). The most abundant dinucleotide repeat motifs were CT having 707 repeats, and the least number of dinucleotide repeat motifs was for CG. A total of six types of dinucleotide repeat motifs (CT, AG, AT, AC, GT and CG) were found in this study. For the trinucleotide repeat motifs, a total of 20 types of trinucleotide repeat motifs were found, where the AAG repeat motif was the most abundant trinucleotide repeat motif found having a total of 103 repeats (16.3%), and the least repeated motif among the trinucleotides was ACG with only one repeat identified. The details of the di- and trinucleotide repeat motifs are summarised in Table 4. Out of all the repeats found, only 1790 were suscessfully designed with their flanking primer. Further, fourteen microsatellites were randomly selected from 1790 microsatellites for characterisation of the Malaysian cultivated taro.

Table 2. Statistical summary of RAD sequencing data for microsatellites.

Clean bases	2,479,455,351
Clean reads	8,490,604
Total Contig	222,445
Contigs containing repeats	3280
Total repeats found	3454
Total successfully design Primer	1790
No. of selected SSR for cultivar identification	14

Table 3. The abundance of repeat motifs type found in this study.

Repeat types	No. of repeats	% No. of repeats
Dinucleotide	2,760	79.9
Trinucleotide	631	18.3
Tetranucleotide	40	1.2
Pentanucleotide	15	0.4
Hexanucleotide	8	0.2
Total	3,454	100.00

Table 4. Abundance of the motif.

Repeat motif	Frequency	% Frequency
CT	707	20.5
AG	701	20.3
AT	612	17.7
AC	414	12.0
GT	311	9.0
CG	15	0.4
AAG	103	3.0
CTT	99	2.9
AGC	51	1.5
GCT	49	1.4
CCT	42	1.2
AAT	41	1.2
AGG	41	1.2
ATT	39	1.1
AAC	37	1.1
GAT	35	1.0
ACC	21	0.6
GTT	21	0.6
ATC	20	0.6
GGT	10	0.3
CGG	8	0.2
CCG	5	0.1
ACT	3	0.1
CGT	3	0.1
AGT	2	0.1
ACG	1	0.0
Others	63	1.8
Total	3454	100.0

Molecular characterisation of Malaysian cultivated Taro

The analysis of 14 novel SSR markers in this study found that there were only twelve polymorphic SSR markers, with the remaining two SSR markers, CesSSR9 and CesSSR11, being monomorphic. Accordingly, these two markers were discarded from the analysis. The analysis of the polymorphic SSR markers revealed the presence of 50 alleles that ranged between two (CesSSR5) and seven (CesSSR13) with an average of 4.17 alleles per locus, whereas, the number of polymorphic information content ranged between 0.05 (CesSSR5) and 0.59 (CesSSR7) with an average of 0.33. Similarly, the value of gene diversity and heterozygosity ranged between 0.06 (CesSSR5) and 0.00 (CesSSR4 and CesSSR5) to 0.63 (CesSSR6) and 0.68 (CesSSR6), respectively. The details of SSR characterisation are summarised in Table 5.

Table 5. Characterisation of novel microsatellite markers.

Marker	Major allele frequency	Allele number	Gene diversity	Heterozygosity	PIC
CesSSR1	0.64	6	0.54	0.24	0.51
CesSSR2	0.82	5	0.32	0.05	0.31
CesSSR3	0.93	4	0.13	0.13	0.12
CesSSR4	0.89	3	0.19	0.00	0.18
CesSSR5	0.97	2	0.06	0.00	0.05
CesSSR6	0.50	5	0.63	0.68	0.56
CesSSR7	0.57	5	0.62	0.18	0.59
CesSSR8	0.86	3	0.25	0.08	0.23
CesSSR10	0.75	3	0.40	0.03	0.36
CesSSR12	0.75	4	0.41	0.16	0.38
CesSSR13	0.67	7	0.50	0.37	0.45
CesSSR14	0.89	3	0.19	0.16	0.18
Mean (polymorphic only)	0.77	4.17	0.35	0.17	0.33

The UPGMA dendrogram (Figure 1) showed that only cultivar MINYAK was unable to separate from cultivar PUTIH. Moreover, the dendrogram showed no variation within cultivar PUTIH. Also shown was an intra-variety variation of cultivars WANGI A1 and WANGI B4. The calculated pairwise genetic distance based on shared alleles showed that the value within the cultivar Putih was 0.00, and the value between cultivar PUTIH and cultivar MINYAK was 0.00. The highest pairwise genetic distance was observed between cultivar TELUR and cultivar WANGI A1 with the value of 0.75 (refer to Table 6).

Figure 1. UPGMA dendrogram generated based on shared allele genetic distance of studied cultivar.

Table 6. Pairwise genetic distance of studied cultivars based on 12 novel microsatellite marker.

OTU	BADAK_A7	CINA_A16	MAWAR_A11	MAWAR_B2	MERAH_60	MINYAK_C17	PUTIH_20	PUTIH_58	PUTIH_G_MUSANG	TELUR	UDANG_C17	WANGI _A1	WANGI_B4
BADAK_A7	0.00	0.07	0.25	0.25	0.12	0.06	0.06	0.06	0.06	0.70	0.61	0.40	0.52
CINA_A16	0.07	0.00	0.27	0.27	0.10	0.02	0.02	0.02	0.02	0.72	0.62	0.43	0.52
MAWAR_A11	0.25	0.27	0.00	0.00	0.29	0.25	0.25	0.25	0.25	0.71	0.67	0.50	0.52
MAWAR_B2	0.25	0.27	0.00	0.00	0.29	0.25	0.25	0.25	0.25	0.71	0.67	0.50	0.52
MERAH_60	0.12	0.10	0.29	0.29	0.00	0.10	0.10	0.10	0.10	0.71	0.53	0.39	0.51
MINYAK_C17	0.06	0.02	0.25	0.25	0.10	0.00	0.00	0.00	0.00	0.71	0.63	0.43	0.52
PUTIH_20	0.06	0.02	0.25	0.25	0.10	0.00	0.00	0.00	0.00	0.71	0.63	0.43	0.52
PUTIH_58	0.06	0.02	0.25	0.25	0.10	0.00	0.00	0.00	0.00	0.71	0.63	0.43	0.52
PUTIH_G_MUSANG	0.06	0.02	0.25	0.25	0.10	0.00	0.00	0.00	0.00	0.71	0.63	0.43	0.52
TELUR	0.70	0.72	0.71	0.71	0.71	0.71	0.71	0.71	0.71	0.00	0.62	0.75	0.70
UDANG_C17	0.61	0.62	0.67	0.67	0.53	0.63	0.63	0.63	0.63	0.62	0.00	0.50	0.62
WANGI_A1	0.40	0.43	0.50	0.50	0.39	0.43	0.43	0.43	0.43	0.75	0.50	0.00	0.55
WANGI_B4	0.52	0.52	0.52	0.52	0.51	0.52	0.52	0.52	0.52	0.70	0.62	0.55	0.00

Discussion

Conventional vs. NGS based microsatellites development

Enhancements in the NGS platform have provided a speedy and financially smart method for microsatellite discovery and progression using paired-end RAD-seq which has subsequently created longer assembly sequences for improving the rates of microsatellite identification. The method has also proficiently and effectively produced a vast series of polymorphic microsatellite markers which can be used for a broad range of research for non-model or understudied species from a hundred to thousands of microsatellites in a responsive and cost-effective way [19–21].

In contrast to the more conventional hybridisation approach for the development of microsatellite markers, the RAD-seq approach, along with de novo sequencing is more efficient and effective considering both time and money. Moreover, the conventional hybridisation approach requires between two and four weeks to complete the technical work, acquiring a restricted number of microsatellite markers. On the other hand, the embraced approach only requires between one and two weeks to obtain a vast number of microsatellite markers [22].

Customary, microsatellite development methods tend to create a sequence from one specimen and fail to consider the existence of repetitive components and insertion-deletion mutations that flank the candidate locus, which can impact the success of both genotyping and precision. Additionally, the conventional hybridisation approach for microsatellite identification and development typically only creates a sequence from a single individual and fails to consider the presence of repetitive components and the insertion-deletion of mutations that flank the target locus thus influencing genotyping success and precision.

On the other hand, there are a vast number of studies on microsatellite development using NGS [23–25]. In a study by Yang et al. [25], they detected the loci of 650 microsatellites resulting from 4.5 million raw RAD reads, with only 285 (43.84%) primer pairs that were designed successfully. However, in the present study, 3,285 microsatellites were found in 222,445 assembled sequences, with 1,790 primer pairs that were designed successfully. Regarding the conventional selective hybridisation approach, according to Badri et al. [26], they only succeeded in obtaining 350 white colonies, with only 52 shown to have unique SSR regions for sesame. Paliwal et al. [27] also applied the conventional approach to creating microsatellite markers, obtaining 356 white colonies, and 114 (32.02%) colonies with SSR repeats. Therefore, this demonstrates the RAD-seq approach for microsatellite development to be more useful and efficient regarding time, cost and labour.

Genetic conservation and cultivar identification

The conservation of taro and its genetic assets requires identification, characterisation and evaluation of the genotypes for the efficient and practical application in breeding activities and acceptable conservation management programme. The methodologies and techniques such as biochemical, morphological, and molecular-based identification of different plant species have been applied worldwide, in which molecular markers have proven to be the best, if not the preferred approach, as the marker system towards recognising the features explicit to plants, which can be recorded with both confidence and precision.

Conventionally, the characterisation of taro depends on the morphological characterisation of its leaves, flowers, fruit, seeds and any other traits. However, even though the utilisation of morphological characterisation for the cultivar identification remains significant, the sole utilisation of morphological characteristics often results in incorrect and conflicting outcomes. These outcomes are attributable to ecological and phenological factors, often limited by the number of different characteristics. As a result, molecular markers offer an increasingly reliable option to that of morphological markers. As in this study, we are able to detect intra-variety variation between WANGI A1 and WANGI B4 by using microsatellite markers, which morphologically we assume both are WANGI cultivars.

Furthermore, despite the expanding ubiquity of various SNP genetic studies, microsatellites remain one of the most informative and flexible markers accessible to various genetic domains. High throughput microsatellite markers are extensively used in plant genomic research given their exceptionally polymorphic, multi-allelic nature, simple logical system and transferability over genotypes [7, 28]. As a result, microsatellite markers are seen as a valued tool for the identification of taro cultivars, providing crucial, if not essential data for sustaining future biodiversity preservation and gene bank management.

Author contributions

S.A.R.: designed the whole project, analysed the data and wrote the manuscript. M.N.G.: sampling and morphological characterization and identification. N.H.E.N.A. and A.M.A.M.: performed the experiment. S.N.I.: scored the alleles and reviewed the manuscript. All authors have read and approved the manuscript.

Acknowledgement

The authors would like to acknowledge anyone who indirectly involved in this project.

Disclosure statement

No potential conflict of interest was reported by the authors.

Data availability

All data that support the findings reported in this study are available from the corresponding author upon reasonable request.

Additional information

Funding

This research was funded by Malaysian government under NKEA EPP14 project.