Sixty finger millet accessions were used in this study. A total of 55 accessions were landraces originally collected from four regional states of Ethiopia (i.e., 15 from Amhara, 15 from Oromia, 15 from Tigray, and 10 from Southern Nations, Nationalities, and People's region) whereas the remaining five were improved cultivars. The landrace accessions and cultivars were obtained from the Ethiopian Biodiversity Institute (EBI) and the Bako Agricultural Research Center (BARC), respectively (Supplementary Table 1). Seeds of each accession were planted in a greenhouse at the Swedish University of Agricultural Sciences (SLU), Alnarp, Sweden. Leaf tissue was collected from young seedlings at about 3 weeks after planting, for DNA extraction. Each accession was represented by 15 individual plants. However, the samples of 29 accessions were collected individually (i.e., each accession was represented by 15 separate samples) whereas the samples of 31 accessions were collected in the pool (i.e., each accession was represented by a single pool of leaf tissue from 15 plants). Hereafter, the two groups of samples will be referred to as “indiv-accessions” and “pool-accessions,” respectively, for the sake of simplicity. For sampling, we used a deep-well plate containing two glass beads in each well. Immediately after sampling, the plate was sealed with a perforated lid, placed in a freeze drier for 4 days, and stored at −80°C until DNA extraction. Frozen samples were homogenized using a Retsch MM400 shaker (Haan, Germany) at 300 Hz for 30 s. DNA was extracted from the homogenized tissue using QIAcube HT (QIAGEN GmbH, Hilden, Germany) following a modified CTAB (Cetyltrimethylammonium bromide) procedure as described previously (Bekele et al., 2007; Tesfaye et al., 2007). A NanoDrop ND-1,000 spectrophotometer (Saveen Werner, Sweden) was used for estimating the quality and quantity of the DNA.

A total of 1,956 finger millet EST sequences downloaded from the National Center for Biotechnology Information (NCBI) were analyzed for the presence of SSRs within the EST sequences using a web-based software WebSat (Martins et al., 2009). A total of 100 sequences that contained SSRs with 2–6 repeat motifs were identified among the downloaded ESTs. After excluding duplicates, overlapping and very short sequences, and sequences with more than one SSRs, 50 EST sequences suitable for designing primers were maintained. Primer pairs were designed for the 50 ESTs targeting the SSRs using a web-based Primer3 primer-designing program version 4.1.0 (Rozen and Skaletsky, 2000). Of note, 10 samples from diverse finger millet accessions were used to evaluate the primer pairs for the quality of their amplified products. This led to the selection of 13 primer pairs that specifically amplified their target SSR loci. These primer pairs were used to amplify the target EST-SSRs (for PCR conditions, see the “SSR-PCR amplification and capillary electrophoresis” section). The amplified products of the 13 loci (five samples each) were purified using the GeneJet PCR Purification Kit (Thermo Fisher Scientific, Vilnius, Lithuania) followed by the Sanger sequencing to confirm that their sequence matches the original EST sequences. The sequencing was performed at Eurofins (www.eurofins.com) using a mixture of 2 μl of 10 μM of a sequencing primer and 15 μl of a purified PCR product. Each amplified product was sequenced using both forward and reverse primers. The comparative analysis of the EST sequences containing these target EST-SSRs through the Basic Local Alignment Search Tool (BLAST) was also performed to evaluate the transferability of the repeat motif across the Poaceae species and to determine the position of the SSRs within the corresponding genes.

A total of 10 genomic SSRs reported to be highly polymorphic in the previous research (Dida et al., 2008) were screened for good amplification, polymorphism, specificity, and suitability for multiplexing using the same set of DNA samples used for testing the EST-SSRs, and three SSR loci were selected. Hence, the 13 EST-SSRs developed in this study and the 3 genomic SSRs from the study by Dida et al. (2008) were used to analyze the 60 finger millet accessions. The PCRs were performed in a volume of 25 μl containing 25 ng genomic DNA, 0.3 μM of each primer, 2 mM MgCl₂, 0.3 mM dNTPs, 1 U Taq polymerase, and 1 × PCR buffer (10 mM Tris–HCl, pH 8.3, and 50 mM KCl). S1000™ Thermal Cycler (Bio-Rad, Hercules, CA, USA) was used for the amplification of the target loci using the following temperature profile: initial denaturation at 95°C for 5 min followed by 35 cycles of 30-s denaturation at 95°C, 30-s primer annealing at optimized temperature for each primer pair (ranging from 55°C to 59°C), and 1 min primer extension at 72°C, with a final primer extension at 72°C for 10 min.

The forward primers were 5′-labeled with 6-FAM™ or HEX™ fluorescent dyes. To prevent the non-template addition by Taq polymerase to the PCR products, the reverse primers were PIG-tailed with GCTTCT according to the study by Ballard et al. (2002). Prior to capillary electrophoresis, PCR amplification was confirmed by running 5 μl of the PCR products on 1.5% agarose gel containing GelRed and visualized using BioDoc-It™ Imaging System (Upland, CA, USA). The multiplexing of the PCR products was performed as described in the study by Geleta et al. (2012). The capillary electrophoresis of the PCR products was performed using an Applied Biosystems 3,500 Genetic Analyzer (Thermo Fisher Scientific, Waltham, MA, USA) at the Department of Plant Breeding, Swedish University of Agricultural Sciences, Sweden.

After capillary electrophoresis, GeneMarker version 2.7.0 (SoftGenetics, LLC, State College, Pennsylvania, USA) software with default settings was used for peak identification at recommended threshold intensity. The determination of the fragment size was based on the GS600 size standard. Each peak was treated as an allele at a codominant locus, and the genotype of each individual/pool at each locus was recorded. Polymorphic information content (PIC) (Botstein et al., 1980) for each marker across all accessions was determined using PowerMarker version 3.25 (Liu and Muse, 2005). Arlequin version 3.5.2.2 (Excoffier and Lischer, 2010) was used for the analysis of molecular variance (AMOVA). The pairwise F_ST output in an XML format was used to generate a graph through the application of a series of R scripts within Rcmd (i.e., a console version of the R statistical package) triggered through the Rcmd command button added to Arlequin version 3.5.2.2 toolbar. The Numerical Taxonomy System (NTSYS) statistical program version 2.1 (Exeter Software, New York, USA) was used to calculate Nei's standard genetic distance (GD) as described by Rohlf (2002). These GD data were used for the unweighted pair group method with arithmetic mean (UPGMA) and neighbor-joining (NJ)-based cluster analyses using the MEGA7 program (Kumar et al., 2016). The Nei's standard GD-based principal coordinate analysis (PCoA) was performed to further display the genetic relationship between the finger millet accessions using GeneAlEx 6.41 (Peakall and Smouse, 2006). The allele data set obtained from the polymorphic loci was used to identify the genetic populations for 29 indiv-accessions using STRUCTURE software version 2.3.4 (Pritchard et al., 2000). The admixture model was adopted with 100,000 burn-in periods and 200,000 Markov chain Monte Carlo (MCMC) chain iterations. To find the optimum number of clusters, a K-value was set from 2 to 10 with 10 independent runs. STRUCTURESELECTOR (Li and Liu, 2018) was used to visualize the threshold K-value (i.e., number of clusters) based on different approaches, the median of medians (MedMedK), median means (MedMeaK), maximum of medians (MaxMedK), and maximum of means (MaxMeaK) (Puechmaille, 2016). After the optimum K-value was determined, the CLUMPACK beta version (Kopelman et al., 2015) was used to display the graphical representation of the population structure.

Thirteen new EST-SSR markers were developed in this study. The list of these EST-SSRs is given in Supplementary Table 2 along with information such as repeat motif, SSR position, the position of forward primer (start), the position of reverse primer (end), expected product size, forward and reverse primer sequences, and the full corresponding EST sequences. Six of these markers are trinucleotide repeat whereas the remaining seven are dinucleotide repeat SSRs (Supplementary Table 2). A comparative analysis of the EST sequences containing the 13 EST-SSRs against annotated genomes of Poaceae species available at https://phytozome.jgi.doe.gov/pz/portal.html (Phytozome version 12.1) showed that, on the one hand, the SSRs of Elco-33, Elco-37, Elco-39, and Elco-41 are located within the coding sequences. On the other hand, the SSRs of Elco-27 and Elco-40 are located in the 5′UTR whereas those of Elco-42, Elco-43, and Elco-48 are located in the 3′UTR (Supplementary Table 2). For example, the SSR of Elco-27 is located 13 nucleotides upstream of the start codon of the aquaporin gene. However, there was no sufficient sequence information to determine the SSR location within the corresponding genes for Elco-14, Elco-35, Elco-45, and Elco-47. Among the 13 EST-SSRs, 6 (i.e., Elco-14, Elco-35, Elco-41, Elco-43, Elco-45, and Elco-47) were monomorphic across the 60 accessions studied. The remaining seven were polymorphic and used for the analyses of population genetics together with the three polymorphic genomic SSRs (Supplementary Table 2). The three genomic SSRs (i.e., UPEG24, UPEG27, and UPEG66) were selected from those used in the study by Dida et al. (2007).

In total, 26 alleles were recorded across the 10 polymorphic loci and 60 accessions of finger millet. The observed number of alleles (Na) per locus ranged from three (in Elco-37 and Elco-42) to nine (Elco-39) with an average value of 5.62 (Table 1). The polymorphic SSR loci accurately followed their repeat motif patterns. For example, locus Elco-27 is a dinucleotide repeat SSR and had 131, 147, and 153 bp alleles (Figure 1) whereas the trinucleotide repeat locus Elco-33 had 158, 161, and 172 bp alleles. The allele pattern of the three samples in Figures 1A–C is an example that demonstrates finger millet is a tetraploid species.

The total gene diversity (H_T) of each locus across all populations ranged from 0.004 in locus Elco-48 to 0.19 in locus UPEG-66 with an average value of 0.09. Among the 29 indiv-accessions, Elco-48 was polymorphic only in two accessions (242622 and 243623), which were originally collected from western Tigray. Similarly, two accessions that were polymorphic only for Elco-48 (accessions 215908 and 215929) were identified among the pooled accessions. These accessions were originally collected from Agew Awi and Gojam, respectively. The within-population gene diversity (H_S) of each locus varied from 0.004 to 0.16 with an average value of 0.08. The observed heterozygosity (H_O) for each locus varied from 0.30 (Elco-48) to 0.81 (UPEG-24). The estimates of population differentiation (G_ST) at each locus varied from 0 (Elco-48) to 0.27 (Elco-37) with a mean of 0.12. Similarly, the lowest (0.045) and highest (0.71) PIC values were recorded for Elco-48 and UPEG-66, respectively, with a mean of 0.51 (Table 1).

The mean gene diversity estimates of each accession across the 10 polymorphic loci showed narrow variation ranging from 0.44 (accession 216054) to 0.68 (accession 237443) (Table 1). Similarly, a narrow range of variation was recorded at the level of regional states ranging from 0.54 (Oromia) to 0.59 (Amhara and Tigray). Interestingly, the average gene diversity of the landrace accessions (0.57) was similar to that of the cultivars (0.58).

The AMOVA based on allele frequency of the 60 accessions (combined data set of indiv-accessions and pool-accessions) was conducted to determine the genetic differentiation at accession and higher hierarchical levels (Table 2). AMOVA that was conducted by grouping the accessions according to their geographic regions of origin and without grouping revealed a highly significant variation among accessions and regions as well as within accessions (Table 2; P < 0.0001). Of the total genetic variation, 18.8% accounted for the differentiation among the accessions (F_ST = 0.18) whereas 81.2% was varied within the accessions. Furthermore, when the accessions were pooled to their respective regions, variation among groups (regions) and among accessions within groups were 6.9 and 13.2%, respectively.

The Nei's pairwise GD among the 29 indiv-accessions ranged from 0 indicating high genetic similarity (215908 vs. 243640; 245088 vs. 215932; 238317 vs. 237443, and Bako-09 vs. 215943) to slightly over 0.5 showing low genetic similarity (Bareda vs. 215929, Bareda vs. 216046, and Boneya vs. 216046). Accession 216046 from Oromia is the most genetically distinct accession with a mean GD of 0.32 from the other 28 accessions. The highest GD between accessions within a region was 0.42 for Amhara (215929 vs. 215943), 0.38 for Oromia (216046 vs. 237969), and 0.31 for Tigray (238321 vs. 242623). The highest GD within the cultivar group (0.41) was recorded for Boneya vs. Addis-01 cultivar (Table 3). The mean GD between accessions within Amhara, Oromia, Tigray, and cultivar groups were 0.13, 0.10, 0.14, and 0.15, respectively. Among regions, the mean GD of 0.15, 0.18, and 0.20 were recorded for Amhara vs. Oromia, Oromia vs. Tigray, and Amhara vs. Tigray, respectively. Against the cultivars, the mean GDs of 0.27, 0.25, and 0.14 were recorded for Amhara, Oromia, and Tigray regions. The pairwise F_ST values demonstrated in Figure 2 are in agreement with the overall significant differentiation (F_ST = 0.18; P < 0.0001) among all accessions studied.

The Nei's standard GD-based NJ dendrogram revealed six clusters with each cluster comprising accessions from different groups (i.e., geographical regions and cultivar) except Cluster-4, which comprised only two accessions from Oromia. Cluster-1 was the largest comprising eight accessions composed of three accessions from Amhara, one accession from Tigray, and four cultivars. Cluster-2 was composed of four accessions from Tigray, two from Oromia, and one from Amhara. Cluster-3 comprised one accession from Oromia and a cultivar. Cluster-5 comprised two accessions from Tigray and one from Oromia, whereas four accessions from Amhara, two from Oromia, and one from Tigray were placed in Cluster-6 (Figure 3).

The PCoA was used to visualize the differentiation among the 29 indiv-accessions (Figure 4). The analysis revealed that the first and second coordinates accounted for 39.8% and 18.9% of the total variation, respectively. Hence, the two coordinates together explained 58.7% of the total variation. Along the first coordinate, a group comprising eight accessions (i.e., three accessions from Amhara, one from Tigray, and four from the improved cultivar group) was separated (highlighted in light yellow) from the major group (highlighted in light green) whereas an accession from Tigray (242621) was an outlier. Similarly, three accessions were separately grouped along the second coordinate. An outlier accession (216041) was also observed along the second coordinate (Figure 4). Overall, the PCoA revealed the absence of a clear clustering pattern of the accessions according to their geographic regions of origin, similar to the results obtained through the cluster analysis (Figure 3).

The analysis of the population structure of the 29 indiv-accessions was conducted based on the data from the 10 polymorphic loci using STRUCTURE software. For this analysis, the approach of Puechmaille (2016) was used, which revealed that the optimum K = 3 (i.e., MedMeaK, MaxMeaK, MedMedK, and MaxMedK = 3), indicating that the individual genotypes of the 29 indiv-accessions most likely came from three genetic populations (Figure 5). The STRUCTURE output at K = 3 displayed partial membership of all indiv-accessions to more than one cluster (Figure 6), suggesting a strong population admixture.

In general, compared with other markers, the EST-SSR markers are appropriate and provide useful genetic information in diversity studies due to their high transferability across species and association with differentially expressed genes (Gupta et al., 2003; Saha et al., 2004; Varshney et al., 2005). In this study, most of the newly developed EST-SSR markers showed high transferability across different crops and they are linked to the coding regions of different genes. For example, the EST-SSR of Elco-27 is linked (only 13 nucleotides upstream of the start codon, data not shown) to the aquaporin gene in maize and rice. The aquaporin gene is responsible for water stress-induced chilling tolerance (Li et al., 2000). The EST-SSR of Elco-42 is linked to a gene coding for a chemocyanin-like protein that plays a positive role in a response to salinity stress and stripe rust in wheat (Feng et al., 2013). In addition to marking aquaporin coding genes and chemocyanin-like protein genes, the newly developed EST-SSR markers are also important for species identification, gene identification, and conservation of genetic material. In this study, 6 of the 13 newly developed EST-SSRs were monomorphic, and those that are monomorphic within E. coracana can probably serve as species-specific markers, differentiating it from other Eleusine species.

Among the seven new polymorphic EST-SSRs, three of them were located in the CDS (coding sequences) whereas four were located in either 3′UTR or 5′UTR. Interestingly, those in the CDS were more informative (mean PIC = 0.40) than the others (mean PIC = 0.24). The three SSRs in the CDS were trinucleotide repeat SSRs and hence do not cause frameshift mutations, which might have contributed to their higher polymorphism. Among the new EST-SSRs, Elco-27 and Elco-39 were the most informative with PIC values of 0.53 and 0.65, respectively. Hence, they need to be prioritized for use in population genetics studies. The level of H_O obtained in this study was similar to that of Bharathi (2011) and Babu et al. (2018) who reported values ranging from 0.20 to 0.85 and from 0 to 0.88, respectively. Given that finger millet is a predominately self-pollinating species, such a relatively high level of heterozygosity is not expected. The possible explanation for high heterozygosity is the fact that it is an amphidiploid, and the SSR loci used are most likely found in both A and B genomes. Hence, although the heterozygosity within a genome can be very low due to self-fertilization, H_O could be high due to the dominance of different alleles in the A and B genomes. Among the three UPEG SSR loci, UPEG24 was mapped to chromosome 3B (Dida et al., 2007). If the locus is specific to the B genome, the high heterozygosity at this locus may suggest that the locus is under selection that favors heterozygosity. This locus was associated with productive tiller number whereas UGEP66 was associated with grain yield and thousand-grain weight (Lule et al., 2018). Hence, breeders can use these markers in their marker-aided breeding programs to increase the productivity of finger millet. It is interesting to investigate if any of the alleles identified in this study at these loci are associated with desirable characteristics of these traits.

The level of genetic variation in populations can be estimated by different statistical parameters, such as H_O, H_T, H_S, percentage of polymorphic loci (PPL), and Nei's gene diversity (H_E). The PIC of an SSR locus indicates the extent of its usefulness in the determination of population genetic diversity. Similar to other population genetics parameters, the PIC of a locus is dependent on the number and frequency of alleles at that locus. The higher the PIC value of a locus, the more informative it is. In this study, a broad range of PIC values (0.05–0.71) was obtained. Interestingly, the average PIC of the seven EST-SSRs (0.31) was lower than the mean PIC of the three previously published SSRs (0.69) (Dida et al., 2008), which is expected as EST-SSRs are more conserved than genomic SSRs in general. The mean PIC values of these three primers in this study (0.69) are lower than previously reported (0.88) (Lule et al., 2018). The higher values reported in the study by Lule et al. (2018) can be explained by the more diverse germplasm used which led to the recording of a larger number of alleles (sum = 61) than those recorded in this study (sum = 20). The gene diversity estimates of the accessions varied from 0.44 (accession 216054 from Oromia) to 0.68 (accession 237443 from Amhara) with an overall mean of 0.58. To identify the local genetic diversity hot spots for finger millet, further study needs to be conducted in the localities where accessions with the above-average gene diversity estimates were obtained. Interestingly, the average gene diversity of the landrace accessions and that of the cultivars were similar, implying that breeding activities did not have a significant effect on the genetic diversity of the crop.

The AMOVA revealed significant genetic variations both among and within the finger millet accessions. A higher percentage of variation (81.23%) was detected within accessions compared with that of among accessions (18.77%), which is in line with the fact that finger millet is a predominantly self-pollinating species. Babu et al. (2014) reported 73% within-population and 27% among-population variation in the study conducted on 190 accessions of finger millet using 75 genic SSR markers. Similarly, Lule et al. (2018) found 69.52% within-population and 30.48% among-population variation in 138 finger millet accessions using SSR markers. Furthermore, research by Pandian et al. (2018) on 83 finger millet genotypes using 43 genic SSR markers revealed 77% within-population and 23% among-population genetic variation. The significantly higher within-population variation than the among-population variation in finger millet could be partly explained by the fact that farmers select finger millet genotypes based on different criteria, such as good performance in grain yield and tolerance to various biotic and abiotic stresses. Only the farmer-selected individuals contribute seeds to the next generation, and hence, the selection of genetically diverse genotypes leads to high genetic diversity within populations. The amphidiploid nature of finger millet is also another contributor to genetic variation within populations. The lower genetic differentiation among populations compared with the variation within populations, which is in agreement with other published studies (Babu et al., 2014; Lule et al., 2018; Pandian et al., 2018), could be explained by strong gene flow through the market channel-based seed exchange and the use of the same improved cultivars among neighboring regions in Ethiopia.

The NJ cluster analysis of the 29 indiv-accessions revealed a weak clustering pattern of the accessions according to the geographic regions of origin of the landrace accessions, indicating that the genetic makeup of the accessions does not have a strong correlation with their geographic origin. The clustering together of accessions from different regions, as well as the placement of accessions from geographically close areas in different clusters in this study, agreed with the results of previous studies. For instance, Pandian et al. (2018) reported 3 clusters for 83 accessions of finger millet using 43 genic SSR markers, Lule et al. (2018) reported 3 clusters for 138 accessions evaluated using 20 SSR markers, and Ramakrishnan et al. (2016) reported 3 clusters for 128 Indian finger millet accessions evaluated using 87 genomic SSR markers. Given that the first two principal coordinates of the PCoA biplot (Figure 5) explained only 58.7% of the total variation, a discrepancy between the clustering patterns of the accessions displayed in the cluster analysis and PCoA is expected. However, the PCoA analysis also revealed a weak clustering of the accessions according to their geographic origin, in agreement with the result of the cluster analysis. The results suggest a strong gene flow between the geographic regions. The fact that most of the cultivars were clustered closely together suggests that they might have been developed based on genetically similar germplasm and/or selected for similar traits during the breeding process. However, Boneya and Adis-01 were described as the best cultivar in terms of grain yield and stability, tolerance to disease, and also other agronomic performance (Negash et al., 2019).

The population structure analysis distinguished three genetic groups (K = 3), which were not based on the geographical origin of the accessions. Some individual genotypes from different accessions showed highly similar genetic profiles and hence have a close relationship. Previous studies, using SSR markers on different finger millet accessions, also produced three genetic clusters having week grouping based on their geographic origin (Dida et al., 2007; Ramakrishnan et al., 2016; Lule et al., 2018; Pandian et al., 2018). All accessions and individual samples analyzed had alleles from the three genetic populations indicating strong gene flow that led to poor genetic differentiation among the accessions. Although there are some exceptions, there is a good general agreement between different analyses in terms of grouping the accessions. For example, the seven accessions in Cluster-6 of the NJ tree (Figure 3) were clustered closely together (highlighted in light yellow) in the PCoA biplot (Figure 4). Similarly, most accessions in Cluster-1 of the NJ tree (Figure 3) were mainly represented by the combinations of purple and blue clusters in the STRUCTURE output (Figure 6B), whereas most of the accessions in Cluster-6 of the NJ tree (Figure 3) were dominated by orange cluster in the STRUCTURE output (Figure 6B). Overall, given the result of the cluster and genetic structure analyses, there are indications of the complex genetic composition of finger millet landrace accessions, and this needs extensive study using different molecular tools, as in other crops (Thurber et al., 2013; Qiu et al., 2014).