Intraspecific Genetic Variability of Brassica cretica Lam. (Brassicaceae) Using SSR Markers

1. Introduction

Resolving the evolution of the Brassica L. (Brassicaceae) crops and their wild relatives probably represents one of the most fascinating tasks in plant evolutionary biology [1,2,3,4]. Darwin [5] already highlighted the parallels between his theory of natural selection and the cultivation practices that shaped the diverse forms of Brassica oleracea L. crops. The domesticated forms of B. oleracea include six major botanical varietes [broccoli (var. italica), Brussels sprouts (var. gemmifera), cabbage (var. capitata), cauliflower (var. botrytis), kale (var. acephala), and kohlrabi (var. gongylodes)], and 12 additional cultivated crop types [6], representing a unique example of diversification through domestication in the plant kingdom. Global production of B. oleracea crops exceeded 100 million metric tons in 2023, indicating its agricultural and economic significance [7].

The generic delimitation of Brassica remains unresolved. Brassica oleracea crops and their wild relatives, however, constitute the monophyletic B. sect. Brassica within the polyphyletic genus Brassica [1]. The section includes about 13 species and several infraspecific taxa distributed across the Mediterranean Basin and the Atlantic coasts of Europe [8]. Several hypotheses on the progenitors of Brassica oleracea crops have been proposed. Brassica montana Pourr., B. rupestris Raf., B. cretica Lam., B. incana Ten., and wild B. oleracea have been nominated as single or multiple progenitor species [9,10,11]. Mabry et al. [6], using RNA-seq data of 224 accessions representing 14 B. oleracea crop types and nine potential wild progenitor species concluded that the Aegean endemic B. cretica is the closest living relative of cultivated B. oleracea, supporting an Eastern Mediterranean origin of cultivation. They also suggested that cultivated plants of this species can revert to a wild-like (feral) state with relative ease. A combination of those results with the evolutionary history of B. oleracea may contribute to a growing body of knowledge on B. oleracea crop domestication that will facilitate continued breeding efforts including adaptation to changing environmental conditions.

Brassica cretica is distributed in the Eastern Mediterranean (Albania, Greece, Türkiye, Lebanon, and Israel), with most of its populations concentrated in the Aegean Islands and the coastal regions of mainland Greece and Anatolia. The extensive phenotypic variability of B. cretica has led to the classification of multiple subspecies and varieties; however, there is no consensus regarding their taxonomic status. Gustafsson et al. [12] identify two subspecies: Brassica cretica subsp. cretica and subsp. nivea. In contrast, Snogerup et al. [11] recognize three subspecies: B. cretica subsp. aegaea subsp. cretica, and subsp. laconica. The disjunct distribution of Brassica cretica in the Middle East has been attributed either to early East Mediterranean trade [13] or to introductions into suitable natural habitats without cultivation [11]. Some populations within the natural range of the species in the Aegean may have also resulted from human-mediated dispersal. According to Mabry et al. [6], early forms of B. cretica may have played underappreciated roles in the domestication of B. oleracea crops. Therefore, resolving the domestication story of B. cretica is fundamental for fully understanding the evolutionary history of B. oleracea.

Brassica cretica is an extremely diverse wild relative of cultivated cabbage (B. oleracea), with genetic distances between its populations often exceeding those found between distinct species in the B. villosa/ B. rupestris/ B. macrocarpa group [14]. This diversity likely results from ‘isolation by distance’ processes, which have shaped significant morphological differences among local populations [11,15]. As it is a diploid (2n = 18) species, partially self-incompatible, hermaphroditic, and insect-pollinated, the geographical isolation of its populations, along with potential pollinator-driven selection based on flower color differences [16], may contribute to its high genetic diversity. Understanding this diversity is crucial for assessing the species’ adaptability and its potential use in breeding programs aiming at enhancing crop resilience [17].

Βreeders of cole crops (B. oleracea L.) have an interest in utilizing current and emerging PCR-based marker systems to differentiate elite germplasm. Among various systems available for genetic analysis in plants, molecular markers are more efficient, precise and reliable in discriminating closely related species or cultivars [18,19,20]. SSRs are co-dominant, highly polymorphic PCR-based markers and thus expected to be a very powerful tool in cultivar discrimination. Polymorphic SSR markers have been recognized as a reliable tool for assessing genetic diversity and population structure in Brassica species in many cases [21,22,23] as well as identifying salinity tolerant genotypes [24]. Furthermore, SSR markers identified in cabbage can be applied to ten other wild and cultivated species within the Brassicaceae family. This highlights their potential value in advancing Brassica oleracea research, including genetic mapping, marker-assisted selection (MAS), and comparative genome analyses [20].

In this study, 10 SSR molecular markers were used to analyze the genetic diversity and population structure of five natural populations of Brassica cretica and one B. oleracea botanical variety to reveal genetic diversity at B. cretica population levels, reveal the distribution of genetic variation within and between the five natural populations, and identify the closest relative among the five B. cretica populations to cultivated species. The genetic information obtained will help to exploit the most promising population in terms of genetic diversity and similarity to cultivated species and thus to explore its future implication in breeding programs.

2. Materials and Methods

2.1. Plant Material & DNA Extraction

This study includes five natural B. cretica populations (two B. cretica subsp. aegaea populations originated from Evia and Ymittos, two B. cretica subsp. cretica populations originated from Akrokorinthos and Crete, and one B. cretica subsp. laconica population from Leonidio) from different regions of Greece and one commercial botanical variety of B. oleracea (population F) (Table 1; Figure 1). Ripe seeds from each population were collected on their natural habitats (Table 1; Figure 1) during physiological ripening period. Seedling preparation was conducted at the greenhouse facilities of the Laboratory of Vegetable Production at Agricultural University of Athens. Specifically, 20 ripe seeds per population were sown on rockwool sheets (AO Plug, Grodan, Roermond, the Netherlands), covered with a thin layer of vermiculite. During germination, seedlings were exposed to ambient light conditions and temperatures ranging from 25 °C to 27 °C during the day, and 17 °C to 19 °C during the night. Additionally, the seedlings were irrigated regularly with a nutrient solution containing 16.4 mM N, 6.6 mM K⁺, 4.4 mM Ca²⁺, 2.8 mM Mg²⁺, 2.75 mM SO₄^2-, 1.3 mM H₂PO₄^-, 20 μΜ Fe (EDDHA), 10 μΜ Mn²⁺, 6 μΜ Zn²⁺, 0.8 μΜ Cu²⁺, 30 μΜ B, 0.7 μΜ Μο, resulting in an electrical conduc-tivity and pH of 2.4 dS/m and 5.6, respectively. One hundred and twenty (120) newly developed seedlings (20 per population) were analyzed in total.

Total genomic DNA was extracted from fresh plant tissue which was stored in deep freeze (–80°C) before DNA isolation. DNA was extracted according to classical CTAB protocol [25] with minor technical modifications. DNA concentrations and purity ratios were measured with Nanodrop ND-1000 Spectrophotometer (NanoDrop Technologies, Inc., Wilmington, DE, USA). Then, all individual DNA samples were diluted at final concentration 50 ng/µl with molecular grade water and stored at freeze −20 °C for PCR analysis.

2.2. Molecular Analysis; PCR Amplification & SSR Markers

The present study used ten SSR markers (Table 2) to analyze genetic variation namely, Ol10B11, Ol10B01, Ol09A01, Ol10F11, Ni4-B10, sORA26, BN12A, Na10-F06, nga111, MB4. Seven SSR primers were studied by [26]; Ol10F11, Ni4-B10, sORA26, BN12A, Na10-F06, nga111, MB4 and five markers described by [23]: Ol10B11, Ol10B01, Ol09A01, Ol10F11, sORA26.

PCR parameters were identical to previous studies with some technical modifications. Selection of those primers was based after an initial screening of a larger number of 14 primers during preparative experiments standardizing PCR conditions according to previously used protocols [23,26].

PCR reactions were performed in 96-well plates Nippon Genetics thermal cycler (FastGene Ultra Cycler Gradient) where different PCRs could be performed at the same time in the same gradient machine using different SSR loci with different Tm (Table 2). The conditions applied for PCRs agreed with corresponding literature [23,26]. The annealing temperature (Ta) was kept 2–3 °C below the melting temperature (Tm) of that primer sequence.

For SSR genetic analysis, PCR amplifications were performed in a total volume of 20 µL for each molecular marker (Eurofins Genomics Custom DNA oligos) separately containing 50 ng/µl of DNA template, 1X reaction buffer (5X KAPA Taq reaction buffer with no magnesium), 2 mM of MgCl2 (Kapa Biosystems) in the case of Ol10B11, Ol10B01and Na10-F06, 2,5 mM of MgCl2 and 0,25 μM of each dNTP (Enzyquest dNTP mix, 10 mM each), 0.4 μM each of the reverse and forward primer, except Na10-F06 and nga111 where we added 0.6 μM of each (Eurofins Genomics Custom DNA oligos) and 0.5 unit KAPA Taq DNA polymerase (Kapa Biosystems, 5u/μl). The first amplification step consisted of initial denaturation at 95 °C for 5 min followed by 35 cycles of denaturation at 94 °C for 35 s, primer annealing at 50 °C - 59 °C (varying with SSR primer pair) for 35 s, primer extension at 72 °C for 35 s and a final extension step at 72 °C for 7 min.

Three technical independent PCR reactions were performed for each SSR marker. The amplified fragments were resolved on 2.2% agarose gel running in 0.5X TBE electrophoresis buffer on device Mupid One Electrophoresis System using 50 bp DNA step ladder (Sigma Aldrich) as size standard reference. DNA bands were visualized under UV light using the gel illumination system Mini Bis Pro (DNRBio-Imaging Systems) connected with software Gel capture.

2.3. Data statistical Analysis

Utilizing binary data generated by SSR primers, the software GenAlex vs. 6.5 was used for the analysis of the results [27].

2.3.1. Variation Within Population

The effective numbers of alleles (Ne), observed heterozygosity (Ho), expected heterozygosity (He), total Expected Heterozygosity (Ht), and Shannon’s Information Index (I) were determined using GenAlex software.

2.3.2. Variation Among Populations

To estimate among-populations genetic variation, Fst (Fixation Index) was calculated to measure the genetic differentiation between the populations based on allele frequencies. For allelic richness (AR) and private allelic richness (pAR) calculations software HP-RARE was used [28].Means and significant values over the loci and populations were obtained by bootstrapping to assess the variability and accuracy of the estimates. The hierarchical distribution of genetic variation (AMOVA) among and within the populations for the ten markers was assessed with GenAlEx software [27], allowing for a partition of the total genetic variation into among-group and within-group components. The effective size of each population (Ne) can be used in the calculation of the Shannon Indices (Shannon’s diversity index) (I) to quantify the richness and evenness of genetic diversity, as well as in the estimation of the effective number of migrants (Nm), which provides insights into gene flow between populations. Nei’s unbiased genetic distance and significance among populations were also estimated, providing a measure of genetic divergence based on allele frequencies. Moreover, genetic distance for each pair of populations was calculated to investigate the relatedness between populations. FreeNA software [29] was used to estimate Nei’s unbiased genetic distance and significance among populations [30], incorporating methods for correcting null alleles. Cavalli-Sforza and Edward’s chord genetic distances (DC) [31], were also estimated for each pair of populations, utilizing the computed methods ENA (DCENA) and INA (DCINA) [29].

The dendrogram among populations was designed with UPGMA (Unweighted Pair Group Method with Arithmetic Mean) based on Nei’s genetic distance using R software. The most likely number of clusters (K) was defined using DELTA Evanno’s method 9 [32] which evaluates the most likely number of genetic clusters based on a model of population structure, using.

3. Results

3.1. Genetic Variation Within Populations

A total of 219 alleles were obtained from the ten SSR loci analysed (Table 3). SSRs, Ol10B11, Ol10B01and Ol09A01 detected >25 alleles. The (total) number of alleles per locus (Na) varied from 16 (minimum Na10-F06) to 32 (maximum; Ol09A01). Respectively, effective number of alleles (Ne) ranged from minimum 2.413 (Na10-F06) to maximum 4,107 (Ol09A01), with a mean value of 3.050 per locus. The number of individuals per population was twenty (20), and the effective number of alleles was 21.9. Expected heterozygosity (He) within population for all SSR loci studied was 0.65.

3.2. Genetic Variation Between and Within Populations

The average number of alleles (Na) per population ranged from 3.200 (population C) to 4.400 (population A), with a mean of 3.650 alleles per individual population. The effective number of alleles (Ne) ranged between 2.747 in population E and 3.519 in the population A, with a mean of effective 3.050 alleles per individual population. Population A, obviously, indicates a higher genetic stability. The expected heterozygosity (He) varied between 0.617 (minimum) (population F) and 0.678 (maximum) (population B). Observed heterozygosity (Ho) was generally noticed in low levels, ranging from 0.060 (population C) to 0.220 (population F). The overall average allelic richness (AR) was 3.26; the highest (3.69) was noticed in population A, whereas the lowest (3.05) was in populations C & E. Private allelic richness (pAR) varied between 0.09 (population D) and 0.86 (population B) with a mean value 0.62 (Table 4).

The genetic differentiation between populations was high (F=0.812; p>0.05). AMOVA test proved that differentiation within-population was 80% and among populations was 20% (Figure 2).

The DELTA method of Evanno [32] was employed to determine the optimal number of clusters (K) (Figure 2), which revealed the presence of four distinct gene pools (K=4) (Figure 3) in the Structure analysis performed using Clumpak software [33].

The STRUCTURE analysis results revealed (optimum) four gene pools in Greece (Figure 3). Group I included B. cretica subsp. aegaea populations from Evia and Ymittos. At group II, populations of B. cretica subsp. cretica from Crete and northern Peloponnese were included. Groups III and IV were formed by B. cretica subsp. laconica and B. oleracea, respectively (Figure 4).

To further examine relationships among populations the INA and ENA methods were used to calculate Pairwise populations Fst Values. (Table 5). The greatest genetic distance was observed between populations C and F (Fst= 0.236) indicating a high genetic divergence, while the smallest genetic distance between populations C and D (Fst= 0.106), indicating their close genetic relatedness.

A UPGMA dendrogram (Figure 4), based on Nei’s (unbiased) genetic distance matrix [30], was constructed using R software and the Ape and Phangorn packages [34,35], to illustrate the genetic relationships among the six studied populations. Populations of B. cretica subsp. aegaea and B. cretica subsp. cretica form two distinct clades. The B. oleracea cultivar forms a clade that is sister to the clade comprising these two subspecies. Notably, B. cretica subsp. laconica appears genetically distinct from the other B. cretica subspecies, forming an entirely separate clade. These findings are largely consistent with the results of the STRUCTURE analysis.

Figure 4. UPGMA tree using Nei’s genetic distance. Bootstrap values are presented at the branch intersections.

Figure 4. UPGMA tree using Nei’s genetic distance. Bootstrap values are presented at the branch intersections.

4. Discussion

Ιn the Mediterranean basin, the fragmentation of species into genetically isolated populations has significant implications for diversification, adaptation, and speciation. The rich Aegean flora is a product of the region’s complex geological history, which has shaped a fragmented insular landscape. Many plant species exist in small, isolated populations with limited migration [36,37,38], contributing to the characteristic mosaic variation of Aegean plant taxa [11,39,40]. Herein, ten SSR markers were used to assess the genetic diversity and population structure of five Brassica cretica populations and one B. oleracea botanical variety. This study aimed to evaluate genetic diversity within and genetic differentiation among B. cretica populations, examine the distribution of genetic variation, and identify the closest relative to the cultivated plants. We found high genetic diversity, with the number of alleles per locus (Na) ranging from 18 (sORA26 and BN12A) to 32 (Ol09A01), with a mean value of 20. Notably, the observed Na was higher than that reported by [26] (Na=5.4) for several shared loci (Ol10-F11, Ni4-B10, sORA26, BN12A, Na10-F06, nga111 and MB4). Furthermore, our study recorded a mean number of private alleles of 0.62, significantly higher than the 0.141 reported in the previous study, indicating a unique gene pool to our study populations.

The average number of alleles per population (Na) in our study was 3.65, which is higher than the 1.9 reported for seven B. cretica populations from Crete [26] and greater than the values found in 25 B. oleracea accessions from Ireland using similar molecular markers [41]. This higher allele number may be attributed to the isolation by distance among our populations compared to the cultivated species studied by El-Esawi et al. [41]. Additionally, the wider geographical distribution of our sampled populations across long-established biogeographical barriers, compared to the more geographically restricted Cretan populations examined by [26], may further account for the observed increase in genetic diversity.

The populations of all B. cretica subspecies showed similar heterozygosity (He) values. The average expected heterozygosity (He) across the studied populations was 0.647, significantly higher than the values reported in previous studies. Edh et al. [26] and Maggioni et al. [42] found a He of 0.211, while Watson-Jones et al. [43] reported a He of 0.26 in their analyses of nine wild Brassica populations, respectively. Allelic richness was 3.26 and private allelic richness was pAR (0.62) higher than the value (0.141) reported by Edh et al. [26]. This higher number of heterozygosity and private allelic richness can be partly attributed to outcrossing breeding systems or the methodology used.

The mean Fst value for our studied populations was 0.81, indicating a high level of differentiation among them. Similar results were reported by Edh et al. [26] using nuclear and chloroplast markers, which revealed exceptionally high levels of population diversity (overall F_ST = 0.628 and 1.000, respectively). In addition, although Watson-Jones [43] reported lower Fst values (0.2257), they also detected considerable differentiation among populations. These results suggest that Brassica species, particularly B. cretica populations, are fragmented across various habitats and have adapted to varying environmental conditions, including soil type, temperature, and moisture levels. Such environmental differences may drive local adaptation and contribute to genetic divergence between populations [26].

Genetic distance analysis and UPGMA dendrograms consistently reveal genetic divergence among Brassica samples, classifying them into four major genetic groups. This grouping is further supported by population structure analysis. Three key conclusions emerge from these findings. Brassica cretica subsp. cretica and subsp. aegaea are genetically distinct but monophyletic, partly supporting the intraspecific classification proposed by Snogerup et al. [11]. While this classification is widely accepted by sources like POWO [44], Euro+Med [45], and the Flora of Greece web (2025), https://portal.cybertaxonomy.org/flora-greece/intro, among others, the classification into two subspecies (i.e., B. cretica subsp. cretica and subsp. nivea) proposed by Gustafsson et al. [12]) is still followed in some studies, e.g., [46,47]. Our results indicate that the population from southern Peloponnese (Leonidio), classified as B. cretica subsp. laconica, significantly deviates from the other two subspecies of B. cretica. The flavonoid composition of B. cretica subsp. laconica also contrasts to that of the other two subspecies [48], suggesting a distinct differentiation in southern Peloponnesian populations. Their taxonomic re-evaluation could provide insights into the evolutionary history of B. oleracea/B. cretica complex, advancing our understanding on the domestication of B. oleracea crops.

The STRUCTURE analysis reveals a high degree of genetic similarity between the populations from northern Peloponnese (Akrokorinthos Castle) and eastern Crete. According to Snogerup et al. [11], populations from both regions belong to Brassica cretica subsp. cretica. Notably, the Akrokorinthos Castle population, located at the acropolis of ancient Corinth with a continuous human presence for millennia, was designated as the type locality for this subspecies [11]. These authors also suggested that the disjunct distribution of B. cretica subsp. cretica in Crete and northern Peloponnese likely reflects remnants of a once more continuous range rather than human introduction. This conclusion was supported by observed chromosomal rearrangements and reduced interfertility between populations. However, our analyses indicate high genetic similarity between the Akrokorinthos and the Cretan populations and reduced genetic variability in the Akrokorinthos population (Ho=0.060), suggesting that an early human introduction to northern Peloponnese cannot be ruled out. Chromosomal rearrangements and reduced interfertility among populations may have resulted from a founder effect, leading to reduced genetic diversity and the fixation of chromosomal mutations, as well as genetic drift, selection, or local adaptation. Interestingly, the Akrokorinthos population has rather low genetic similarity to B. oleracea cultivar, possibly indicating an introduction to suitable habitats without cultivation as suggested by Snogerup et al. [11].

The Brassica oleracea cultivar clusters closely with B. cretica subsp. cretica, particularly with the Cretan populations, which exhibit the highest genetic similarity. This observation agrees with previous findings by Mabry et al. [6], who identified B. cretica as the closest extant relative of domesticated B. oleracea, lending further support to the hypothesis of an eastern Mediterranean origin for cultivated cabbages and kales. Historical evidence also suggests an eastern Mediterranean origin for B. oleracea. Descriptions of kales, a form of B. oleracea, first appear in the works of the Greek scholar Theophrastus (370–285 BC), whereas earlier texts from 800–600 BC make no mention of such plants. This temporal context implies that domestication likely took place in this region during that period [49].

Using SSR markers we were able to detect the genetic diversity in five B. cretica populations. Our findings indicate the need for additional molecular and phenotypic studies to clarify the taxonomic classification of B. cretica in the Aegean region. These investigations should include a broader range of populations from different regions of Greece and the eastern Mediterranean. Our results also revealed high genetic diversity and allelic richness in B. cretica populations, making them a valuable genetic resource for breeding with Brassica crops.

Climate change threatens both the survival of native Brassica populations and the sustainability of Brassica crop production. Rising temperatures disrupt flowering times, pollination, and seed viability, while extreme heat stress reduces germination and biomass production [49]. Drought and water stress further hinder seed germination, early development, and overall yield, with wild populations facing habitat loss due to desertification [50]. Extreme weather events further damage plants, displace seeds, and accelerate soil erosion, while habitat fragmentation and genetic erosion reduce adaptive potential. Climate-induced pollinator declines and shifts in flowering periods disrupt reproduction [51]. To mitigate these threats, conservation efforts must focus on preserving genetic diversity, protecting natural habitats, and implementing sustainable agricultural practices. To our knowledge, very few accessions of B. cretica currently exist in genetic banks in Greece. Therefore, establishing ex situ conservation measures for Greek populations, including seed collection and vegetative propagation, should be a high priority.