Genetic characterisation and population structure analysis of Anatolian figs (Ficus carica L.) by SSR markers

We used 310 fig accessions which were collected from different eco-geographical sites (Figure 1) in Anatolia and deposited at the National Fig Germplasm Repository, Fig Research Institute, Erbeyli-Aydın, Turkey. The names, accession numbers, locations and the gender of the figs studied here are presented in Supplementary Table 1.

Figure 1

DNA was extracted from leaf tissue as described by Lefort et al. (1998) and its concentration was estimated spectrophotometrically as described in detail by Akçay et al. (2014) and Burak et al. (2014).

Fourteen SSR markers, namely MFC1, MFC2, MFC3, MFC8 (Khadari et al., 2001), LMFC23, LMFC25, LMFC30 (Giraldo et al., 2005), FCUP068-1, FCUP038-6, FCUP008-2, FCUP070-2, FCUP027-4, FCUP066-7 (Bandelj et al., 2007) and FM4-70 (Zavodna et al., 2005) were used in this study (Supplementary Table 2). PCR amplifications were performed as previously described by Akçay et al. (2014). Briefly, 15–200 ng DNA, 0.5 mM dNTP, 5 pmol of labelled forward and reverse primers, 0.5 unit DNA polymerase (Promega) (containing 1.5 mM MgCl₂) and 1 ml 10X buffer were used and a total of 10 ml of PCR reaction mixture was prepared. PCR program: 1 cycle (94°C for 3 min), 35 cycles (94°C for 1 min, 53–58°C depending on the primers for 1 min, 72°C for 2 min), followed by 72°C for 10 min and 4°C forever (Supplementary Table 2). The amplification control of PCR products was checked by 2% agarose gel electrophoresis.

Capillary electrophoresis conditions were previously described by Akçay et al. (2014) and Burak et al. (2014) in detail. A Beckman CEQ fragment analysis software was used to determine the allele size of each SSR locus. The analyses were repeated at least twice to ensure that the results are reproducible. In each run, ‘Sarılop’ and ‘Kadota’ cultivars were included as controls.

For each locus, the number of alleles (n), allele frequency, expected (He) and observed heterozygosity (Ho) and the probability of identity (PI) values were calculated as previously described (Akçay et al., 2014; Burak et al., 2014). The software IDENTITY was used to detect identical accessions; the proportion of shared alleles was calculated using Microsat (version 1.5) (Minch et al., 1995) and a dendrogram was constructed with the unweighted pair-group method with arithmetic mean (UPGMA) method (Sneath and Sokal, 1973), using the software NTSYS-pc (Numerical Taxonomy and Multiware Analysis System) (version 2.0) (Rohlf, 1988).

Excluding the Beyaz Bukele accession which showed tri allelic cases for SSR loci, the Arlequin software (Excoffier et al., 2005) was used to estimate the population genetic parameters of 309 diploid fig accessions. A neighbour-joining tree constructed from Nei's genetic distances was used (1972). Genetix 4 (Belkhir et al., 1996–1998) was used to perform factorial correspondence analysis (FCA) and gene flow estimates. We used the BAPS (Bayesian Analysis of Population Structure) software (version 6.0) (http://www.helsinki.fi/bsg/soft-ware/BAPS) (Corander et al., 2008) to analyse individual elements of data from each province to distinguish population structures. The STRUCTURE software (Pritchard et al., 2000) was employed to analyse population structures of fig accessions. In these analyses, the same computing parameters were used with the exception of measuring the K level for K:1-8 for unknown reconstructed panmictic populations (RPPs) with 25 replications as reported by Pereira-Lorenzo et al. (2018). Structure Harvester (Earl and von Holdt, 2012) was also used for the estimation of the best K value supported by the data (Evanno et al., 2005). We used Bayesian model-based clustering methods to identify RPPs. The number of accessions strongly assigned to each RPPs was determined based on the qI (probability of membership) probabilities greater than 80%.

For clone differentiation, the GenAlEx v6.5 program (Peakall and Smouse, 2012) was used to identify multilocus genotypes (MLGs) in populations. Using the same program, number of different alleles (Na), effective alleles (Ne), observed heterozygosity (Ho), Nei's (1978) unbiased expected heterozygosity (uHe) and private alleles summary (PAS) values for each population were determined.

In addition, a histogram of pairwise distances generated using the software GenoType v1.2 (Meirmans and van Tienderen, 2004) was used to determine whether somatic mutations are present. A possible number of clones (representing clone number) and Simpson's diversity based on multilocus lineages (MLLs) calculations was conducted using the GenoDive v1.1 program (Meirmans and van Tienderen, 2004). In addition, an effective number of genotypes (accessions) (eff), genotypic diversity (div), eveness (eve) and Shannon–Wiener (shw) diversity index values were calculated using the GenoDive v1.1 program. In the analysis of MLGs, different mutational threshold or T values (T shows the maximum distance allowed to identify a clone between two individuals with the same ‘MLG (accession)’ value) were tested (e.g. from threshold = 0 to threshold = 10) to minimise potential scoring errors and mutational problems. The groupings of MLGs within MLLs were considered, and the accessions with similar mutational threshold values were considered to represent the clones.

In this study, a total of 310 fig accessions were analysed for 14 SSR loci and a total of 124 alleles were identified. The lowest number of alleles (2) was observed for the LMFC23 and the highest number of alleles (15) was found for the FCUP038-6 locus. The average allele number per locus was 7.75 and the mean He and Ho values were 0.652 and 0.685, respectively. The highest Ho values of 0.845 and 0.819 were observed for MFC1 and LMFC30, and the lowest values of 0.319 and 0.383 for MFC8 and LMFC23, respectively. FCUP038-6 with 15 alleles (PI: 0.093), FCUP070-2 with 11 alleles (PI: 0.094) and LMFC30 with 9 alleles (PI: 0.093) were the most informative loci. LMFC23 (PI: 0.600) with two alleles was the least informative locus (Table 1).

Various genetic parameters, such as Ho and He, polymorphic loci and the mean number of alleles per locus, estimated for six fig accession groups, are summarised in Table 2.

The highest gene flow (Nm) (31.23) was found between Central Anatolia and Marmara accession groups and the lowest gene flow (2.94) was found between Aegean and Black Sea accession groups (Table 3). Genetic similarity between the accessions was estimated using the coefficient for Nei's standard genetic distance (1972) (Supplementary Table 3), which showed relatively high genetic similarities. The tree constructed using neighbour-joining analysis was consistent with the findings from genetic distance analyses (Supplementary Figure 1).

Marmara accessions showed relatively high similarity to the geographically close Aegean and Central Anatolian accessions with high gene flow rates (10.16 and 31.23, respectively) between these accession groups (Supplementary Figure 1, Table 3 and Supplementary Table 3). Southeast Anatolian accessions showed relatively low similarity to the accessions from other regions, especially to Aegean and Central Anatolia groups, as evidenced by high genetic distance and low gene flow values. Overall, genetic distance and adjacent joining analyses revealed that genetic similarities between the groups were high.

The FCA (Figure 2) showed partial sub-structuring of fig accession groups. Mediterranean and Black Sea accessions were grouped separately and showed relatively little overlap with other accession groups. Southeast Anatolian accessions were similar to Mediterranean and Black Sea accessions and showed a partial overlap with these two groups (Figure 2). In contrast, Aegean and Marmara accessions showed a strong overlap with accessions from Central Anatolia, most likely due to potential gene flows occurring between these populations (Table 3).

Figure 2

Southeast Anatolian accessions, which showed relatively high genetic distance from the accessions of other regions, displayed a homogenous (monocoloured) BAPS structure (Supplementary Table 3 and Supplementary Figure 2). Central Anatolia, which showed overlap with other regions in the FCA, formed a homogenous structure according to the BAPS analysis. Marmara accessions, which showed the low genetic distance to other populations with a significant overlap in the FCA, displayed an admixture BAPS structure (Figure 2; Supplementary Figure 2).

In this study, we identified 22 homonymous (genetically different accessions known by the same name), 36 synonymous (genetically same accessions known by the different name) and 7 identical accession groups (Supplementary Tables 4–6).

A maximum K value at K = 2 corresponding to the two main RPPs (Figure 3; Supplementary Figure 3; Supplementary Tables 1 and 2) (RPP1 and RPP2) was identified using STRUCTURE analyses of diploid fig populations. Furthermore, RPP analyses showed that 87% (271) of accessions could be assigned to individual RPPs with at least 80% probability whereas 13% of accessions either could not be assigned or could be assigned only with a low probability value to a representative RPP (Figure 3).

RPP1 contained 81% of the Aegean population, followed by Central Anatolia (55%) and Marmara (51%) populations. In RPP1, 127 out of 165 accessions (77%) were observed with a probability of membership ratio qI (>80%), whereas a total of 38 accessions (23%) of different populations were identified as <80%. In RPP1, the highest number of accessions (126) was found in the Aegean population whereas the lowest number of accessions (4) was found in the Mediterranean population (Figure 3, Supplementary Table 1).

Figure 3

In RPP2, the populations with the highest numbers of accessions were Southeast Anatolia (100%), Mediterranean (93%) and Black Sea (79%). In RPP2, the qI (>80%) value of 72% was similar to that of RPP1. In RPP2 populations, the highest number of entries was found in the Mediterranean (39 accessions) and the lowest number of accessions (5 accessions) was found in Central Anatolia. In addition, the qI value was <80% for 40 accessions (28%), which included 4 accessions from the Black Sea, 3 accessions from Southeast Anatolia, 10 accessions from Marmara, 4 accessions from the Mediterranean, 2 accessions from Central Anatolia and 17 accessions from Aegean populations.

In both RPPs, RPP1 and RPP2, a total of 70 accessions were identified as admixed accessions (qI < 80%). The populations with the highest number of admixed accessions for both groups were Aegean (35 accessions) and Marmara (22 genotypes) (Figure 3; Supplementary Table 1).

When male–female distribution was examined in RPP analyses, female figs were distributed in two RPPs, whereas 40 out of 45 male (caprifig) figs originating from Aegean-Aydın and Aegean-İzmir, respectively, were found in RPP1, and the remaining 5 in RPP2.

In MLG analyses, a total of 96 different MLGs (accessions) and 213 unique accessions were identified. In line with the number of populations, the highest number of MLGs was 38 and these were found in the Aegean population; on the other hand, the lowest number of MLGs was 3 and these were found in the Central Anatolia population. The mean number of different alleles (Na) and effective alleles (Ne), which were 5.87 and 3.35, respectively, were similar in all 6 populations. Observed heterozygosity (Ho) values were 0.68 in Marmara, Black Sea and Southeast Anatolia populations. These Ho values were higher than the uHe values of the same populations. However, in the Aegean population, Ho and uHe values were identical. The lowest PAS value with 4 alleles at 2 different loci was found in the Black Sea and Southeast Anatolia populations, whereas the highest PAS value with 12 alleles at 5 different loci was found in Aegean and Mediterranean populations (Table 4).

In MLLs analyses performed using different threshold values, some minor differences were observed in the number of different MLLs. There was no difference in the number of MLLs between threshold = 2 and threshold = 3 values. Therefore, threshold = 2 is considered as the threshold value.

Based on clonal diversity values, all populations except Central Anatolia and one accession of Southeast Anatolia were found to contain unique accessions (Table 5).

The genotypic diversity (div) value, also known as expected heterozygosity, was similarly found to be high in all populations (mean 0.991). The lowest eveness (eve) value, which shows distribution status of accessions within a population, was 0.702 in the Aegean population, whereas the highest eveness value, which was the same as the number of effective number of genotypes (accessions) (eff), was 1.000 in the Central Anatolia population. This indicates that all accessions of the Central Anatolia population have the same frequency. The higher Shannon–Wiener (shw) value (2.035) which was observed in the Aegean population as compared with other populations, was consistent with the high-diversity features of this population (Table 5).

Among accessions, 54 clones/multilocus accessions were found for threshold = 2. Especially in Aegean (33 clones) and Mediterranean (9 clones) populations, accession-clone matches were found to occur mostly within the same population. Detailed information (accession name, accession no, geographical region and province) on clones determined based on threshold = 2 value is given in Supplementary Table 7.

It has been observed that the clone matching groups determined in the MLL analysis are partially similar to the groups determined in the RPP analysis.

Based on the PI values, the most informative locus was FCUP038-6 (PI: 0.093; 15 alleles), and the least informative locus was LMFC23 (PI: 0.600; 2 alleles). The FCUP038-6 locus has also been identified as the locus which demonstrates the highest number of alleles in previous investigations (Bandelj et al., 2007; Caliskan et al., 2012; Ferrara et al., 2016; Caliskan et al., 2018). The number of alleles per locus ranged from 2 for LMFC23 to 15 for FCUP038-6 and the allele sizes ranged from 121bp for MFC3 to 261bp for LMFC30. Similarly, the LMFC23 locus displayed low polymorphism in other studies (Giraldo et al. 2005; Giraldo et al., 2008; Ferrara et al., 2016). The MFC and LMF group loci (Khadari et al., 2003; Saddoud et al., 2008; Achtak et al., 2009; Aradhya et al., 2010; Do Val et al., 2013; Caliskan et al., 2018; Ganopoulos et al., 2015; Teoman et al., 2017) and the FM4-70 locus (Zavodna et al., 2005; Caliskan et al., 2012; Ikten et al., 2018) have also been used by other workers for genetic identification of figs and, overall, our results are in accordance with these previous studies. In a study by Teoman et al. (2017), 24 LMF loci were studied on 45 caprifigs (F. carica var. caprificus) and 2 female figs from the Marmara and Aegean regions of Turkey, and the LMF-30 locus was found to have the highest number of alleles per locus (9 alleles). Similar results were also obtained by Achtak et al. (2009) and Aradhya et al. (2010).

Estimates of genetic similarity were obtained from SSR markers data. The genetic distance matrix varied between 0.055 and 0.232 suggesting that fig accessions analysed had relatively high genetic diversity. The lowest similarity value was found between the accessions from Southeast and Central Anatolia whereas the maximum similarity was found between Marmara and Aegean accessions. Additionally, high similarities between accessions from Marmara and other regions were found (Supplementary Table 3). The Marmara region is known for having the highest rate of human commute due to its historical past (the former capital of Turkey throughout history), and this genetic similarity of Marmara accessions and other regions can be due to the comparatively high level of genetic material exchange with these regions throughout history.

Based on SSR analyses, 310 accessions were classified into six accession groups (dendrogram not shown). Similarly, Caliskan et al. (2018) divided 90 caprifig accessions from the Eastern Mediterranean into five groups. In our study, reference accessions from geographically distant Aegean region clustered separately from eastern accessions. The most significant genetic differentiation value (0.078) was found between Aegean and Black Sea fig accession groups whereas the lowest F_ST value (0.007) was identified between Marmara and Central Anatolian fig populations. Similarly, significant genetic differentiation values were found within Anatolian fig populations, with the highest and lowest F_ST values being 0.182 and 0.007, respectively (Caliskan et al., 2018). The results of gene flow analyses revealed that the accessions from geographically distant regions had lower levels of gene flow, whereas the accessions originated from closer geographical locations had higher levels of gene flow.

The observed heterozygosity values, which were higher than expected heterozygosity for all groups except for those from Central Anatolia, indicated a high level of cross-pollination in different F. carica populations. This finding agrees with that of Aradhya et al. (2010), who studied fig samples from Europe, Asia and North America. According to the BAPS structure analysis, fig accessions from Southeast Anatolia displayed a homogeneous structure which is consistent with the FCA analysis. In addition, Southeast Anatolian accessions showed the highest genetic distance among the remaining regions studied (Supplementary Table 3). Central Anatolian accessions showed a homogeneous structure. Based on genetic distance and the rate of gene flow analyses, genetic similarity between accessions from this region and those of Marmara indicated gene flow from Central Anatolia to the Marmara region. Marmara and Aegean accessions with an admixture (heterogeneous) structure had some overlap with other accessions (Figure 2; Supplementary Figure 2), suggesting that gene flows to these regions from other regions had been also occurred.

The RPP analyses confirmed the results of structuring by other population analyses (FCA, BAPS structure, gene flow, etc.). RPP1 included 81% of the Aegean population, half of the Marmara population and 50% of the Central Anatolia population. As a result, these three populations overlapped in FCA, showing high gene flow and low genetic distance to one another.

RPP2 was the most diverse RPP and consisted of accessions from all populations. In particular, the Mediterranean and Southeast Anatolia, which showed high similarity to each other, were included in this group (Figure 3; Supplementary Tables 1 and 2).

In MLG analysis carried out in our research, 96 MLGs that showed correspondence to 31% of the total population were distinguished. In similar MLG analyses performed on different plants, this rate was 78% in Halophila ovalis population (Xu et al., 2019), 48% in olive population (Barazani et al., 2014) and 14% in wild population of Ziziphus celata (Rhamnaceae) (Gitzendanner et al., 2012). The rate of MLGs in a population may vary significantly depending on the total population number, the number of SSR markers used and discrimination power (Hamadeh et al., 2018). Accordingly, about 56% of 96 MLGs determined in MLG analysis were distinguished as clone in MLL analysis (threshold = 2) (Table 5) and proved that the discrimination power of SSR markers used in the clonal analysis is sufficient.

As it is known, heterozygosity values (unbiased expected heterozygosity (uHe) and observed heterozygosity) determined in MLG analysis provide information about genetic diversity and kinship relations between MLGs (Harris and DeGiorgio, 2017). In MLG analysis performed on 1,747 Pueraria montana (Lour.) Merr. var. lobata (also known as kudzu) plants belonging to 87 different locations in the United States, the average number of samples for 87 locations was 20, and based on the unbiased expected heterozygosity (uHe) value (0.398) and observed heterozygosity value (Ho) (0.444), it was clear that kudzu plant had high genetic diversity and it was explained that this situation was caused by high clonal reproduction (Bentley and Maurıcıo, 2016). In addition, in MLG analysis based on SSR markers in 429 genotypes belonging to 22 different populations of the perennial tree Platanus orientalis living in the Mediterranean region, unbiased expected heterozygosity (uHe) was found between 0.267 and 0.607, and observed heterozygosity (Ho) was found between 0.207 and 0.564, and despite the high number of analysed population, low heterozygosity values observed at the clonal level were attributed to geographic isolation and low gene flow (Rinaldi et al., 2019).

In MLG analysis of our study, it was interestingly revealed that expected heterozygosity (uHe) and observed heterozygosity (Ho) values were quite close to each other in all six populations, and values of both heterozygosities were observed to be 0.620 and higher. This situation reveals the possible genetic variation among MLGs in each population (Harris and DeGiorgio, 2017).

As explained in detail later, in this study, 36 synonymous, 22 homonymous and 7 identical fig accessions were identified (Supplementary Tables 4–6).

In general, the ecotypic variation due to somaclonal variations of fig cultivars appears to be high. Ecotypes here are defined both within and between geographical regions and provinces. Especially in Aegean, Marmara and Mediterranean ecotypes, different accessions were found homonymous (e.g. ‘Datça 1’ (145) – ‘Datça 3’ (147) – ‘Datça 4’ (148) – ‘Datça 5’ (157)), whereas similar accessions (‘Bektaşi’ (40) – ‘Mor Güz’ (95) – ‘Mor İncir’ (107), etc.) were found synonymous (Supplementary Tables 4 and 5).

Clones may be different in MLGs for reasons such as allele scoring error and somatic mutations. Somatic mutations, especially in perennial plants, can lead to differences between individuals with the same clonal genotype (Coughlan et al., 2017). Clonal variations were mostly detected in different varieties that do not form ecotypes within and between geographical regions/provinces. Clonal diversity detected among non-ecotype varieties is thought to be caused either by variations found in each geographical region/province accessions (‘Mor Armudi’ (172-Mediterranean, Adana), ‘Kırmızı İncir’ (174-Mediterranean, Adana), ‘Kırmızı İncir’ (187-Mediterranean, Adana)) or by the emergence of new clones, especially in case of other Aegean and/or Marmara varieties (‘Hamri’ (182-Mediterranean, Adana), ‘Siyah’ (215-Marmara, Balıkesir), ‘Siyah İncir’ (216-Marmara, Balıkesir)) grown in other ecologies (Black Sea, Mediterranean, etc.) (Supplementary Table 7).

In recent years, RPP analyses have been used by various researchers in figs for population structure groupings (Belttar et al., 2017; Teoman et al., 2017; Ikten et al., 2018).

In our study, overall RPP-based groupings showed correlations with genetic similarity relationships of the 310 accessions examined. Regardless of the geographical region/province, almost all 36 synonymous, 7 identical accession cases were included in the same RPP (RPP1 or RPP2) (Supplementary Tables 1, 5 and 6).

In addition, the finding that synonymous, identical accession are grouped with the same RPPs is an indication of their unique structures. Clonal cases except for a few clone groups (unnamed (306) – Sarı Bardak (227); Yediveren (249) – Dumanı Kara (218) and Mor İncir (200) – Kilis inciri (194)) have been found to be in the same RPPs groups (Supplementary Tables 5–7).

Common fig is a gynodioecious species where female and hermaphroditic trees are found together within a population. Male fig trees produce syconia containing both male and female flowers. In contrast, the syconia produced by female fig trees contain only female flowers. Because only male trees can produce pollen, the common fig is considered to be functionally dioecious (Boudchicha et al., 2018).

A total of 45 male accessions from Aegean-Aydın (37) and Aegean-İzmir (8) were included in this study. So far, the number of studies conducted on the structure distribution of male and female accessions in fig accessions has been limited. Teoman et al. (2017) reported that the structure analysis of 47 fig genotypes produced 2 RPPs without any clear male–female separation.

In contrast, in our study, caprifigs showed a distinct structure and mostly located in RPP1, although only five male accessions from the Aegean-Aydin male group (‘Siyah İlek 48’ (48), ‘Kizilay 1’ (58), ‘ Şeytan 1’ (64), ‘Taşlık’ (73), ‘Kızılburun’ (74)) were also found in RRP2.

One of the reasons for these five male accessions to differ from the others could be due to their different genetic structures. Another reason why these five male accessions differed from others would be related to plant genetics as well as gender phenotype. An orthologue of RAN1 gene loci was reported to be associated with sex determination in fig (Mori et al., 2017). Mori et al. (2017) reported that mutational variations of the two SNP regions in RAN1, which is associated with sex determination in figs, caused heterozygous accessions in the 18 male fig accessions. In the current study, five accessions previously reported to be female (Ikten and Yilmaz, 2019) were identified as male based on RAN1-associated cleaved amplified polymorphic sequences (CAPS) marker screenings. These heterozygous genome variations may be another reason for the fact that 5 of the 45 caprifigs accessions were found in different RPPs.