Genetic diversity and population structure of eddoe taro in China using genome-wide SNP markers

Plant varieties and DNA extraction

A collection of 234 taro accessions from 16 different provinces of China, including 16 multi-corm type and 218 multi-cormel type (Fig. 1, Table S1), were selected for SNP identification in the present study. The multi-corm and multi-cormel types of taro can be distinguished by the morphology of their underground corms and cormels (Fig. 2). These accessions were cultivated in separate rows with plant spacing of 50 cm and labeled with the accession name to avoid mechanical mixture on the experimental farm of the National Germplasm Repository for Aquatic Vegetable at Jiangxia District, Wuhan, China (30°20′, 114°14′; at an elevation of 32 m). Fresh healthy leaves from each taro accession were sampled, frozen in liquid nitrogen immediately and stored at −80 °C until DNA extraction. DNA was extracted according to the CTAB method. The concentration and quality of DNA samples was quantified with a NanoDrop-2000 spectrophotometer, and DNA samples were adjusted to 50 ng/µL.

SLAF-seq library preparation and sequencing

The SLAF library was perpared based on the procedures described by Sun et al. (2013) with some modifications. Genomic DNA of the accessions was digested with restriction enzymes HaeIII and EcoRV, with rice (Oryza sativa cv Nipponbare) DNA as a control to evaluate the efficiency of enzyme digestion. DNA fragments that were 364–414 bp (with indexes and adaptors) in length were chosen as SLAF tags to prepare paired-end (2 × 125 bp) sequencing with the Illumina High-Seq2500 system (Illumina, Inc.; San Diego, CA, US) at Beijing Biomarker Technologies Corporation.

Raw read processing, SNP calling and filtering

After filtering out adapter reads and low-quality reads, the sequencing quality was then evaluated by calculating the guanine-cytosine (GC) content and the Q30 value (Table S2). The clean reads were uploaded to NCBI SRA database (http://www.ncbi.nlm.nih.gov/sra/) under the BioProject ID PRJNA644023. Sequences of sample with the most tags were used as the reference. BWA 0.7.17 (Li & Durbin, 2009) was used to align reads to the reference. GATK 3.7.0 (McKenna et al., 2010) and SAMtools 1.9 (Li et al., 2009) were used for SNP calling. And the intersection of SNPs called by GATK and SAMtools were considered as potential high-quality SNPs. Samples with a coverage less than 10 folds were removed, and SNPs with >0.2 missing data rate and minor allele frequency (MAF) <0.05 were removed by VCFtools 0.1.16 (Danecek et al., 2011). After that, 132,869 high-quality SNPs (Table S3) were identified, which were used to perform the following analyses.

Population structure, principal component analysis and maximum likelihood (ML) tree

Before running ADMIXTURE, SNP data set was screened by Hardy-Weinberg equilibrium (HWE), linkage disequilibrium (LD) and outlier loci. SNPs with p value of HWE less than 0.0001 and LD decay with r² great than 0.1 were filtered out through PLINK v1.90b4 (Purcell et al., 2007), and outlier loci were detected by BayeScan v2.1 (Fischer et al., 2011) and Lositan (Antao et al., 2008) with default settings. Population structure was estimated by ADMIXTURE v1.3.0 (Alexander, Novembre & Lange, 2009) with ten independent simulations for each K ranging from 1 to 10. A 10-fold cross-validation (CV) procedure was performed in order to calculate the CV errors for each K, and the value of K that exhibited a lowest CV error was selected as the optimal number of populations. Additionally, all the accessions were assigned to a corresponding group based on their ancestry fractions. The R package ‘pophelper’ that implements CLUMPP was used to plot the population structure (Francis, 2017).

To complement the results from ADMIXTURE, the DAPC function of the “adegenet” R package (Jombart, 2008) was used to determine the optimal number of clusters and assign individual accessions to the different clusters with the 132,869 high-quality SNPs. The optimal number of clusters from DAPC were determined by the K value corresponding to the lowest Bayesian Information Criterion (BIC). Thereafter DAPC clustering was performed on the clusters identified using the 40 principal components. Principal component analysis (PCA) was performed with TASSEL v.5 (Bradbury et al., 2007). A maximum likelihood (ML) phylogenetic tree was produced using RAxML-NG v.1.0.0 (Kozlov et al., 2019) under the GTRGAMMA model with 200 bootstrap replicates, and the best ML tree was displayed through iTOL v4 (https://itol.embl.de/) (Letunic & Bork, 2019).

Genetic diversity and development of core germplasm

Nei’s genetic diversity index (H), polymorphic information content (PIC), MAF, and observed heterozygosity (Ho) of the population and inferred groups were calculated with the popgen function in the R package “snpReady” (Granato et al., 2018). The analysis of molecular variance (AMOVA) was assessed using the R package “poppr” (Kamvar, Tabima & Grünwald, 2014; Kamvar, Brooks & Grünwald, 2015), and the function ‘randtest’ was performed with 999 replicates for permutation test. Pairwise levels of differentiation and nucleotide diversity (π) were estimated using the Populations program in Stacks pipeline (Rochette & Catchen, 2017), and one thousand bootstraps for all loci were set to build 95% confidence intervals for pairwise FST. The R package Core Hunter 3.0 (De Beukelaer, Davenport & Fack, 2018) was employed to develop a core germplasm set with the biparental genomic matrix that contained the total SNPs. The genotypic data of the core germplasm accessions were screened out with VCFtools for further analysis.

SLAF-Seq genotyping

A total of 234 taro accessions, including 218 of multi-cormel type and 16 of multi-corm type from 16 provinces of China (Table S1), were surveyed for SNP identification in the present study. Through SLAF-seq, 1774.09 million of clean reads were obtained from this experiment, and number of SLAF tags ranged from 231,112 to 600,960 with an average sequencing depth of 18.04 (Table S2). The average of Q30 reached 93.50% (Table S2) and the average of GC was 45.11% (Table S2) among 234 samples. SNPs were developed for 234 samples based on the polymorphic SLAF tags. After filtering out the invalid SNPs, 132,869 high-quality SNP markers (Table S3) were retained for further analysis.

Population structure analysis

In all, 25505 SNPs passed the filters with HWE of p value less than 0.0001 and LD with r² threshold of 0.1. BayeScan detected 13 outlier SNPs which were the subset of the 7,582 ones that identified by Lositan. After screening by these filters, 17,923 SNPs were set to ADMIXTURE for population structure analysis. The estimated ancestry fractions of the 234 accessions for different values of K ranged from 1 to 10 were determine d by ADMIXTURE (Fig. 3). When K = 2 ∼ 5, the SNP panel set cannot separate the two morphological types (TMC and TNC) clearly. From K = 6 on, the genetic substructure between TMC and TNC progressively takes shape. The cross-validation error (CV) decreased to the lowest value at K = 8 (Fig. 4), which indicated that the entire population could be categorized into eight inferred groups. According to membership probability values, six, 12, 54,12, 84, 49, seven and 10 accessions were assigned to group I, II, III, IV, V, VI, VII and VIII, respectively (Fig. 3, Table 1). Accessions in group III, IV, V, VI and VIII are all multi-cormel type, while accessions of multi-corm type dominate in groups I and II, and only one of seven accessions in group VII belongs to the multi-corm type. Grouping of the 234 taro accessions was not according to their collected regions (Table 1).

The Bayesian information criterion (BIC) from discriminant analysis of principal components (DAPC) was performed to determine the appropriate number of clusters. The BIC showed rapid decline from 1 to 8 (Fig. S2), indicating that the accessions can be grouped into eight main clusters, consistently with the results obtained using ADMIXTURE (Fig. 5A, Fig. S1 and Table S1). To gain deeper insight into the accessions grouping and the pattern of variation, principal component analysis (PCA) was performed to measure the variation in taro collection based on the 132,869 SNP markers. In our study, the first two principal coordinates explained 51.95% of the total variations. It revealed that the eight groups inferred by ADMIXTURE were clearly differentiated in the top two eigenvectors, which accounted for 35.52% and 16.43% of the total variation, respectively (Fig. 5B, Fig. S2). To determine the genetic relationship among taro accessions, a maximum likelihood (ML) phylogenetic tree was constructed based on the SNPs unveiled. As Fig. 6 exhibited, the ML tree can be classified into eight clusters, which showed high consistency with the assignments made by the ADMIXTURE program.

Genetic diversity analysis

To evaluate the genetic diversity of the eight inferred groups, genetic parameters including nucleotide diversity (π), Nei’s genetic diversity index (H), observed heterozigosity, PIC, and MAF were estimated separately. The value of π and Nei’s genetic diversity index indicated that accessions in group V exhibited the highest genetic variation, whereas the accessions in group IV and VIII showed the lowest genetic variation (Table 2). Consistently, the highest Nei’s genetic diversity index was present in group V (0.474), and the lowest was found in groups VIII (0.353). Furthermore, The mean PIC values ranged from 0.276 to 0.356, with the highest polymorphism in group V, while the lowest value in group VIII; and mean MAF values revealed similar tendency for the eight groups inferred, which ranged from 0.278 in group VIII to 0.425 in group V and VI.

Population differentiation

Pairwise F_st and analysis of molecular variance (AMOVA) was performed to measure genetic differentiation among the eight inferred groups. The results revealed the highest differentiation between group IV and VII, with an F_st of 0.298. On the other hand, the lowest Fst of 0.100 were found between group VI and II, and group VI and VII, which indicated that they were most closely related (Table 3). The hierarchical AMOVA for all retained SNPs indicated that 77.94% of the total variance was attributed to genetic differentiation among the eight inferred groups, while 22.06% of the diversity was attributed to the genetic differentiation within inferred groups, suggesting a rich genetic diversity level among groups, while 20.46% variance was found among two types and less variance among surveyed 16 regions (11.69%) (Table 4). This is in agreement with the results of population structure, where the inferred groups of taro accessions were not clustered accordingly to their collected regions.

Development of the core germplasm set

A core germplasm set comprising 41 accessions that included eight multi-corm accessions and 33 multi-cormel accessions, was screened out using Core Hunter (Table 1, Table S1). The genetic diversity of this core germplasm set was estimated to determine the extent to which it might effectively represent the genetic diversity of the entire collection. Comparisons of all genetic parameters revealed that Nei’s genetic diversity index, π, Ho, PIC, and MAF values for the core germplasm set were equal to or greater than those for the entire collection (Table 5). Specifically, mean Nei’s genetic diversity, PIC and MAF values of the entire collection were 0.374, 0.299 and 0.279, respectively. However, these values increased to 0.388, 0.307 and 0.295 for the core germplasm set, respectively. Also, the mean π and Ho values of the core germplasm set exhibited no significant difference to that of entire collection. The AMOVA result showed that there was no significant difference between the developed core collection and the rest of the collection, and that the ‘within collection’ differences contributed 99.71% of the total genetic variation, indicating that the core collection represented the entire collection well (Table 6).