Exploring Genetic Diversity in Romanian Lycium Species: A Comparative Study of Wild and Cultivated Genotypes

1. Introduction

Goji berry plants have been used in Chinese medicine and cuisine for thousands of years [1] and presently, the berries are acknowledged as one of the most recognized 'superfruits' of the 21st century [2,3,4], being considered a functional food [1,5]. The goji berry has attracted significant attention in Western countries due to its nutritional profile, especially for its abundant vitamins and antioxidants. Its ORAC (Oxygen Radical Absorbance Capacity) values, which lie between 25,000 and 30,000, surpass those of other nutritionally beneficial fruits like pomegranates and blueberries, indicating its superior antioxidant capacity [6]. Its medicinal uses go from improving visual acuity [5,7,8], infertility [5], abdominal pain [5], dry cough, fatigue, and headache [5], immune system support, cancer prevention [7,8], antidiabetic activity [7,8] to increased longevity [8,9,10] and enhanced fertility [10,11,12,13].

In China, out of the existing nine Lycium taxa [14], only four are traditionally utilized, with L. barbarum and L. chinense being the main species traded worldwide [14,15]. For over two millennia, various pharmacological substances have been collectively referred to as "gou qi" [4,5,16], although it is likely that a diverse range of genetic materials is hidden behind this generalized term. The genus Lycium, distributed in temperate and subtropical regions worldwide [17,18,19], actually comprises 436 species names [20]. Out of these, only 92 are acknowledged as accepted species. 241 of these names are considered synonyms, while the remaining 103 names are categorized as unplaced [20]. Yao and al. (2018), cites 97 Lycium species, out of which 35 species and 2 varieties with recorded uses of food and/or medicine worldwide [14]. The Plants of the World Online records 101 officially recognized Lycium species, indicating their indigenous presence in 71 countries, across 130 regions [21]. The database also indicate a list of 91 countries and 73 regions where Lycium species have been newly introduced, including Romania as one of these recent additions [21].

The flora of Romania recognizes Lycium halimifolium L. as a native species [22], however, a manuscript from 1867, referenced by Marian (2008), documents the traditional usage of Lycium vulgare Dun. in Romania, also mentioning its identification as L. barbarum in the Transylvania region [23]. Although it is scientifically acknowledged that L. halimifolium is typically considered a synonym of L. barbarum [24,25,26], widespread public belief, likely shaped by online media and local traditions, treats L. halimifolium and L. barbarum as separate species, attributing different culinary and toxicological properties to them [27,28]. In fact, both L. halimifolium and L. vulgare are today synonym names of L. barbarum [25,29]. Evidence of the long-standing presence of wild goji berry species in Romanian regions is found in the variety of local names used historically. These include names like ‘licie’, ‘lician’, ‘lițian’, ‘răchișoară’, ‘răchițică’, and ‘zaharică’ in Bucovina, ‘licină’ in Transylvania, and ‘cătină’ and ‘cătină de garduri’ in Moldova region [23]. The plant is characterized as a shrub with red berries, which not only serves as a natural fence in the countryside, but it is also recognized for its medicinal properties. Traditionally, it has been used to treat conditions related to fear and anxiety, demonstrating psychological effects, and for epilepsy and spasms, indicating neurological benefits [30]. In few Romanian regions it is considered an invasive plant, as in Oltenia, Danube riverbanks and Dobrogea [31,32].

Due to goji berry increased fame, the market demand grew exponentially in the last two decades [33,34]. China dominates goji berry production, particularly in the Northwest regions, like Ningxia and Xinjiang, two main exporting regions [9,35,36]. In contrast, production in North America and Europe is limited due to a lack of traditional use and the challenge of adapting Chinese varieties to these regions [37,38,39,40]. Ideal U.S. cultivation spots include Northern California and parts of the Northwest, while in Europe, Italy and Romania are leading in production [41,42]. Italy boasted the largest European goji farm in 2013, with 35 hectares and 60,000 plants [6]. Romania has emerged as a significant producer of goji berries [41], focusing also on plant material for cultivation [43], with a market that is showing a rising trend [42]. Especially in the difficult context of climate change constraints, goji berry planting material adapted to local conditions is required by European farmers. Therefore, initiatives on identifying new cultivars and crop production processes were launched, together with breeding programs that considered Lycium as an important target [44]. In Romania, by 2023, seven varieties of Lycium spp. were registered on the Official catalogue of cultivated plant varieties: 'Erma’, ‘Transilvania’, ‘Kirubi', '‘Kronstadt‘, ‘Bucur’, ‘Sara’ and ‘Anto’, belonging to both L. barbarum and L. chinense [45].

Genetic diversity is the cornerstone of biodiversity, offering significant value in numerous areas, both scientific and practical [46]. It plays a crucial role in comprehending the origins and evolutionary processes of species, conserving genetic resources, investigating valuable germplasm resources and use the genetic materials displaying valuable traits [47,48]. Having a deeper understanding of native goji berry genetic resources is important both for preserving local biodiversity and for the breeding sector [36,47,49]. Since the rise of goji berry market demand, extensive molecular studies have been undertaken, to determine valuable genes in both cultivated and native or spontaneous plants, as a resource for valuable traits in the future breeds [1,18,36,49,50,51].

The process of crop plant breeding focuses on identifying and cultivating new plant varieties that exhibit beneficial characteristics like higher yields, better disease resistance, and enhanced nutritional value [52]. The high-throughput technologies, including genomics, transcriptomics, and metabolomics, have opened a new phase in crop breeding, enhancing the efficiency and precision of this process [52,53]. The last two decades have seen a significant growth in both the volume and quality of publicly available plant genomes [53]. This progress is largely due to advancements in sequencing technologies such as PacBio Sequel, Illumina HiSeq, Oxford Nanopore PromethION, and Illumina NovaSeq, alongside bioinformatic tools like the chromosome conformation capture technique (Hi-C) and BioNano optical maps. These developments have greatly enhanced the efficiency of genome sequencing, assembly, and annotation [53,54,55].

In the Solanaceae family, which includes around 90 genera and 3,000 species, 170 full genomes of 46 species have been reported [54]. These genomic studies primarily targeted horticultural crops and their wild relatives, such as the cultivated tomato (Solanum lycopersicum) and its wild counterpart (Solanum pimpinellifolium), potato (Solanum tuberosum), hot pepper (Capsicum annuum), and eggplant (Solanum melongena). They also extend to tobacco (Nicotiana tabacum), decorative flowers like Petunia inflata and Petunia axillaris, and herbs with medicinal properties, including Datura stramonium and Lycium barbarum [52,56], or its invasive relative, L. ferocissimum [57]. L. barbarum genome contains 12 chromosomes [34] (2n = 2x = 24) and is 1.8 Gb in size, with a level of heterozygosity of approximately 1% [56]. The sequenced and annotated genome ASM1917538v2 [58] was obtained by sequencing a haploid plant developed from pollen culture, using PacBio Sequel technology [56]. The annotation allowed identification of 47,740 genes and 34,339 protein-coding sequences. The availability of another annotated genome, of L. ferocissimum, of 1.2 Gb size, 40,291 genes and 30,549 protein-coding genes [57] will ease even more the characterization of the future goji berry sequenced genomes, allowing identification of new genes of interest.

Exploring the genetic diversity in Romanian wild and cultivated genotypes of Lycium species could bring the breeding sector one step closer to rapidly developing new varieties, adapted for specific ecopedological conditions and market requirement, indicating potential genetic markers for required traits. Generally, the goji berry in Romania has three major biotic threats, the powdery mildew, the goji berry gall mite and the stink bugs [59,60], being much more suitable for organic production than in China [40,61]. The escalating threat of extreme weather events caused by climate change is set to pose an increasingly serious challenge to goji berry production, major threats being extreme drought and insolation [62,63]. The Romanian homologated varieties were developed based on Chinese varieties germplasm, due to their high fruit quality traits [64,65], without using the local germplasm.

The current study marks the initial phase of a broader project focused on genes related to resistance to abiotic and biotic stress, focusing here on cuticle thickness. This aspect is critically important due to cuticle significant influence on resistance to drought, insolation, and pests [66,67,68], as the cuticle is essential for plant survival in periods of severe drought [69] and pests and disease attack [70]. BODYGUARD genes in plants typically encode for proteins that are directly involved in defending the plant against various stresses, including the biosynthesis of cutin in Arabidopsis [69,71,72], mainly the α/β hydrolase protein [69,71,73]. Analysing the genetic diversity of cultivated and wild goji plants genes has the final aim of providing information required by the goji berry breeders supporting the development of goji berry production in Romania.

2. Results

2.1. NGS Data Analysis

2.1.1. Sequencing Data Quality Control

The genomes of eight Romanian L. barbarum varieties, out of which five were cultivated varieties, part of a population obtained from Chinese seeds [74] and three were spontaneous plants growing in the wild in three different Romanian counties [31,75] were sequenced using NGS technology. The quality distribution of sequencing was analysed across the entire length of all sequences to identify any locations with abnormally low sequencing quality, that could indicate the inclusion of incorrect bases at higher-than-normal rates. The proportion of clean data relative to raw data, referred to as the effective rate, was higher than 99.52% for all reads. For 1,669,720,889 bp reference genome, the mapping rate of each sample ranged from 96.66% to 99.36%. Referring to the reference genome (without Ns), the average depths were between 10.01 X and 9.29 X, and 1 X coverages ranged from 77.43% to 97.41%, the results falling within the acceptable normal range and could be utilized in variation detection and genetic analysis.

2.1.2. SNP Detection, Distribution, and Mutation Frequency

SNP (Single Nucleotide Polymorphisms) variations were observed in all eight genotypes, but the quantity and genomic distribution of these variations differed across the genotypes. A total of 108,290,958 SNPs were identified within the eight genotypes, with an average ranging from 14,079,300.6 SNPs/genome for the cultivated specimens and 12,631,485 SNPs/genome for the wild specimens. The Lb2 genome exhibited the largest quantity of SNPs, totalling 15,983,773.

In the eight genotypes, the transitions – point mutations that change one purine nucleotide to another or one pyrimidine to another – were more frequent, with an average count of 8,559,837.5. This was higher than the number of transversions, which are mutations that switch a purine for a pyrimidine or vice versa, averaging at 4,976,532.25. This resulted in an average ts/tv ratio (transitions to transversions) of 1.72, which was relatively consistent across the two categories of genotypes, cultivated and wild. However, there was a difference between cultivated and wild-grown plants: the average ts/tv ratio was higher in cultivated genotypes at 1.736, compared to 1.688 in wild-grown plants.

The genotype Lb7w exhibited the highest heterozygosity rate, measured as 5.037‰ (per thousand), whereas the lowest rate was found in genotype Lb4, at 3.245‰. On average, wild plants showed a higher heterozygosity rate, averaging at 4.851‰, compared to cultivated plants, which had an average heterozygosity rate of 3.607‰.

Figure 1 illustrates the distribution of the six types of SNP mutations. Genotypes Lb1 and Lb2 showed the highest number of SNPs across all six SNP types, while genotype Lb8w had the lowest count. Among these six types of SNP mutations, for every genotype, the most common was the C:G > T:A mutation, followed by T:A > C:G. Conversely, the least frequent type of SNP was the C:G > G:C mutation.

2.1.3. Insertion/Deletion Detection and Distribution

InDel (Insertions and Deletions) variations were detected in all the studied genotypes, amounting to a total of 11,225,960 InDels, averaging 1,403,245 InDels per genome. The genotype Lb2 had the highest number of InDels at 1,773,325, whereas Lb8w had the fewest, totaling 1,197,605. For cultivated specimens, the average InDels per genome was 1,495,506.4, compared to 1,249,476 for wild specimens. When examined individually, the overall count of insertions, which is 5,166,205, is less than the total total number of deletions, totaling 6,043,098. However, the mean count of insertions in cultivated specimens, at 684,984.8, exceeds the average insertion count in wild specimens, which is 580,427. In the case of deletions, the average number of deletions in cultivated specimens, amounting to 808,410.6, is greater than the average deletion count in wild specimens, recorded at 667,015.

The InDel heterozygosity rate, expressed in per mille (‰) and calculated as the ratio of InDels to the total genomic bases, is 0.381‰. This value is lower on average for cultivated specimens, at 0.335‰, compared to the higher average rate for wild specimens, which is 0.458‰.

InDels distribution within the genome (Figure 2) was the highest for 1 bp insertion/deletion and decreased with the increase in sequence length.

In the analysis of the eight genotypes, the densities of SNPs and InDels across each chromosome appeared to be relatively similar. However, a visible reduction in InDel density was observed (Figure 3). Additionally, distinct diferences were evident between the genomes of cultivated and wild plants. Generally, the ratio of SNPs to InDels was around 10. This ratio was slightly lower in cultivated plants, ranging from 9.29 to 9.52, and somewhat higher in wild plants, with values ranging from 10.03 to 10.16.

For both SNPs and InDels, it is visible that the density of variation is higher at the end of the all 12 chromosomes. In addition: (1) all genomes of cultivated specimens have a high variation density, for both SNPs and InDels, almost in the middle

2.1.4. Sequence analyses of BODYGUARD genes in Romanian goji berry genoms

In both Lycium species available genomes, three BODYGUARD (BDG) genes have been identified (Table 1). However, these genes are situated on different chromosomes in each species, on chromosomes 4, 8, and 9 in Lycium barbarum, and on chromosomes 1, 3, and 9 on Lycium ferocissimum.

Analyzing the eight studied genomes using Genome Workbench, a distinct divergence was noted between the genomes of cultivated and wild plants.

Regarding the BDG gene situated at LOC132634709 (Table 2) on Chromosome 4, 22 SNPs were identified within its coding region. Out of these, 14 are synonymous mutations, meaning they do not change the amino acid sequence. With the exception of the SNP at the position 1312, situated within a codon that encodes for lysine/arginine (basic amino acids), the rest of the SNPs are situated within codons that encode for either nonpolar (10) or polar amino acids (11), and none of these SNPs changed the polarity or charge of the encoded amino acid. A key finding was the apparent distinction in most SNPs between cultivated and wild varieties of the plant, differences being observed as homozygous versus heterozygous SNPs.

In the BDG gene at LOC132607278 (Table 3) on Chromosome 8, analysis revealed 28 SNPs within its coding sequence, including 4 synonymous mutations. Notably, SNPs at positions 287-288 and 392-294 each impact a single codon, changing AAA to CGA (Lysine to Arginine) and AAA to GAC (Lysine to Aspartic Acid), respectively. A distinct pattern was observed between cultivated and wild genotypes in 17 specific SNP locations: 174, 209, 288, 392-394, 617, 648, 688, 694, 799, 891, 1030, 1049, 1076, 1150, 1151, 1266, and 1304, predominantly presenting as homozygous versus heterozygous SNPs. Apart from the SNPs at positions 146 and 1274, which alter the charge of the encoded amino acids (Lysine to Glutamic Acid and Glutamic Acid to Lysine, respectively) from basic to acidic and vice versa, the other SNPs do not cause any changes in either the polarity or charge of the encoded amino acids.

The BDG gene at LOC132609965 (Table 4) on chromosome 9 features 29 SNPs within its coding region, with 18 of these being silent mutations. Two SNPs at positions 555-557 lead to the formation of four different codons: TCT, GCT, TCA, and TCG, which correspond to the amino acids Serine, Alanine, Serine, and Serine, respectively. In a similar manner, the two SNPs at positions 1671-1672 result in the codons GGA, AGA, and GAA, which encode for the amino acids Glycine, Arginine, and Glutamic Acid, respectively. The SNPs at positions 555-557 are unique in that they alter the polarity of the encoded amino acid from Serine (as found in the reference genome) to Alanine (as observed in the cultivated genotypes). The other SNPs, however, do not cause any changes in the polarity or charge of the amino acids they encode. Approximately half of these SNPs show sequence-level differences between the cultivated and wild type genotypes.

3. Discussion

The sequencing and annotation of the goji berry genome in 2023 [56] represents a crucial resource for future resequencing projects. The advancements in Next Generation Sequencing/Whole Genome Sequencing (NGS/WGS) [55] are poised to generate a wealth of data, which will be instrumental in developing new goji berry varieties. The examination of variations in density at the genomic level, particularly for SNPs and InDels, highlights specific genome areas that warrant further investigation, to identify potentially beneficial genes from wild genotypes that could be integrated into new varieties. The observed increased polymorphism near the ends of chromosomes can be attributed to the higher frequency of recombination in these areas [76].

In Romania, seven new goji berry varieties have already been created and registered, originating from Chinese genomic material, without genetic elements from native plants. Present research is dedicated to enriching the diverse gene pool found in wild germplasm, potentially enhancing the unique characteristics of Romanian goji berries. By examining the morphological and phenological traits of wild goji berries and correlating them with genetic data, characteristics like early or late flowering, high drought tolerance, strong resistance to low temperatures (Figure 4), as well as features like thicker cuticles and leaves, could become valuable assets in breeding programs.

In earlier research, the morpho-anatomical features of leaves and flowers of both wild and cultivated goji berries in the Bucharest region were analysed. The study aimed to identify key traits of interest to both goji berry breeders and taxonomists [31]. An other study also involved mapping the spontaneous genetic resources found across Romania [75]. Notable morphological distinctions were observed in leaf shape, orientation, and width of Romanian L. barbarum, results that are similar with findings reported in Republic of Moldova, in a similar study between cultivated and wild goji berry [77]. The leaf anatomical characteristics are particularly significant in relation to biotic and abiotic stress factors, the wild plants having the leaves covered with a thick cuticle, prominently developed vascular bundles, and sheaths surrounding the vascular bundles within the mesophyll. Additionally, the palisade cells in these plants were observed to be considerably larger than in cultivated plants [31]. These findings motivate further investigation into genes putatively linked to these phenotypic differences. The formation of the plant cuticle involves several proteins that play crucial roles in the biosynthesis and regulation of cutin and waxes, as BDG, CER, KCS, VLCFAs, GPAT, LACS, ABC, SHN/WIN, LTPs, CD1 [69,71,78,79].

The sequence analyses of the BDG genes in Romanian goji berry genomes revealed several differences among the three genes. The genes located on chromosomes 4 and 9 are encoding probable lysophospholipase BODYGUARD 3 proteins, whereas the gene located on chromosome 8 encodes a probable lysophospholipase BODYGUARD 4 protein [56]. The BDG 4 gene from chromosome 8 is shorter than the BDG 3 genes, from chromosomes 4 and 9. Although located on different chromosomes, two of the genes presented SNPs affecting the same amino acid, as in 12, 65, 84, 235, 254, 410, 426, 467, and 473 positions. Even if some SNPs are located within conserved regions, many of them are silent (Figure 5).

4. Materials and Methods

4.1. Plant Material

In this study, five selected genotypes from Romanian goji berry cultivated and three spontaneous-growing genotypes were examined. The five genotypes are integral to an extensive breeding program for goji berries that commenced in 2014 at the Experimental Field of the Faculty of Horticulture, at University of Agronomic Sciences and Veterinary Medicine in Bucharest [64,65,74]. The initial biological samples were derived from the seeds of Lycium barbarum L., including five distinct biotypes: Lb1-Lb5 [74]. The native plant samples were chosen from robust and well-established populations in the counties of Bucharest (Lb6w), Ilfov (Lb7w), and Călărași (Lb8w). Specifically in Bucharest, specimens were gathered from the shores of Morii Lake (44.453424, 26.013337), a natural area on the periphery of the western segment of the Romanian capital. This location was also selected for a comparative morpho-anatomical study of the leaves and flowers of both wild and cultivated goji berry plants [31]. The plants encountered in Ilfov county are believed to have originated from cultivated specimens within a military base, subsequently becoming naturalized in the area. (44.447382, 26.019239). The specimens from Călărași were found proliferating along a roadside on Lehliu city (44.434389, 26.858775).

4.2. DNA Extraction

Genomic DNA from fresh goji berry leaves was isolated using the InnuPure C16 automated system (Analytik Jena GmbH, Jena, Germany), which employs magnetic particle separation technology for the fully automated extraction and purification of DNA. This process took place at the Research Center for Studies of Food Quality and Agricultural Products at the University of Agronomic Sciences and Veterinary Medicine in Bucharest, Romania. For genomic DNA extraction was used the InnuPREP Plant DNA I Kit-IPC16 (Analytik Jena GmbH, Jena, Germany), adhering to the protocols provided by the manufacturer. Initial processing involved breaking down the plant material externally, with the sample being mashed into a fine powder under liquid nitrogen and then homogenized with a SLS lysis solution (with CTAB as detergent), proteinase K, and an RNase A solution. Following this external lysis step, the automatic DNA extraction continued in the InnuPure C16 automated system, as per the manufacturer's guidelines. DNA quantification was performed with the NanoDrop™ 1000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE) [80].

4.3. Sequencing

Whole-genome sequencing (WGS) was conducted using the next-generation sequencing (NGS) technology of an Illumina platform by Novogene Co., Ltd. (Cambridge, UK). The original image data from Illumina's high-throughput sequencing were converted into sequenced reads (raw data) through the CASAVA base recognition process (Base Calling) at Novogene Co., Ltd. These raw data were saved in FASTQ (.fq) format files [81], which included the sequencing reads along with their respective base quality scores.

4.4. Computational Data Processing, and Sequencing Analysis

The BWA software was utilized to align the effective sequencing data with the reference sequence, using the parameters: mem -t 4 -k 32 -M [82]. The alignment outcomes were used to calculate the mapping rate and coverage. To ensure the quality of subsequent analyses, reads of low quality or those containing adaptors were filtered out to achieve clean reads. The Phred quality score (Qphred), which is indicative of the base quality, was determined by the formula Qphred = -10log10(e), with ‘e’ being the error rate of the sequencing.

The reads were aligned with the reference genome of the goji berry, GCF_019175385.1, downloaded from the NCBI database [56]. This process produced sequence alignment format files, which were subsequently transformed into binary sequence alignment format (*.bam) files. These were then processed to generate a variant file containing SNP (Single Nucleotide Polymorphism) data. The mapping rates for the samples indicate the degree of resemblance between each sample and the reference genome. Additionally, depth and coverage serve as metrics for the consistency and extent of correspondence to the reference genome, as conducted by Novogene Co., Ltd.

4.5. SNP Detection and Annotation

SNP detection was conducted using SAMtools with the specified parameters 'mpileup -m 2 -F 0.002 -d 1000' [82], facilitated by Novogene Co., Ltd. To minimize the likelihood of errors in SNP identification, the data underwent a two-step filtration process: a SNP was only considered if it was supported by over four reads, and its mapping quality had to exhibit a root mean square value exceeding 20, based on the supporting reads' mapping qualities. The overall heterozygosity rate of SNPs across the genome (het rate, denoted in permille ‰) was determined by the number of heterozygous SNPs over the total count of genomic bases. SNPs were assorted into six mutation classifications: T:A > C:G, T:A > G:C, C:G > T:A, C:G > A:T, T:A > A:T, and C:G > G:C. Take, for instance, mutations from T:A to C:G, which entail alterations from T to C and A to G. A T-to-C mutation on one strand of the DNA double helix will correspond to an A-to-G mutation at the identical position on the opposite strand. As a result, mutations of T > C and A > G were grouped together in one category.

4.6. Insertion/Deletion (InDel) Detection and Annotation

An InDel was identified as either an insertion or a deletion of a DNA sequence that is 50 base pairs (bp) in length or shorter. Detection of InDels was carried out using SAMTOOLS with the parameters set to ‘mpileup -m 2 -F 0.002 -d 1000’ [82], annotated with ANNOVAR software [83], by Novogene Co. The criteria for filtering InDels to enhance detection accuracy were consistent with those applied during SNP detection. The length distribution of InDels was examined as a proportion of the entire genome.

4.7. Sequence Analysis of the BDG genes

Next-generation sequencing BAM files containing the nucleotide sequence data for the eight goji berry genotypes were uploaded onto the NCBI Workbench software and aligned to the reference genome [84]. For each variety, the differences in nucleotide sequence were noted. Amino acid sequences of the three probable BDG proteins were aligned using the MultAlin software [85].

5. Conclusions

The present study re-sequenced the whole genome for eight L. barbarum genotypes, cultivated and wild-type, and analysed the variability of three BDG genes, involved in the cuticle biosynthesis, at the coding sequence level. NGS sequencing revealed clear differences between cultivated and wild type genotypes, not only in the whole genome, but also among the BDG genes. Future studies will be done to confirm the role of the DBG genes in cuticle biosynthesis, and furthermore, their implication in the resistance to biotic/abiotic stress. In addition, the data generated by the whole genome resequencing of these genotypes will allow for the analysis of additional useful genes that could be introgressed from the wild type into future varieties in the goji berry breeding programs in Romania.