Assessment of Genetic Diversity and Phylogenetic Relationship of Local Coffee Populations in Southwestern Saudi Arabia using DNA Barcoding

1. Introduction:

Coffee is one of the most commercially significant crops, and the second most traded commodity after oil [1]. In addition to its high export value, coffee has also gained in cultural significance over the past few decades. Despite there being more than 125 reported species in the genus Coffea, only two species, Coffea arabica L. (also known as Arabica coffee) and C. canephora Pierre ex A. Froehner (known as Robusta coffee) are grown commercially [2]. The total annual global coffee production in 2020 was over 105.21 million tons, with 63.15 million tons of Arabica coffee accounting for 60.02% of the total, and 42.06 million tons of Robusta coffee comprising 39.97%. Coffee's genetic development is progressing at a sluggish pace despite its enormous economic relevance. The collection, characterization, and wise use of accessible germplasm material for any crop plant species contribute to its genetic development and long-term viability [3]. Therefore, enhancing diversity from both local and foreign sources is critical for the improvement of crops [4]. For historic reasons, the main issue with Arabica coffee has been its narrow genetic base that limits its adaptation to changing environments [2]. To get around this problem, breeders made use of wild coffee diploid species to introduce new genes into Arabica genotypes [2]. For instance, the leaf rust resistant Arabica cultivar Timor Hybrid got its resistance from its C. canephora parent; it was later used as a parent to develop several new rust resistant cultivars such as Catimor and Ruiru 11 [5]. For bean and liquor quality traits, the wild tetraploid Arabica genotypes from the species’ center of origin and the little-known ancient varieties from the Arabian Peninsula offer a wide gene pool to explore [6]. Despite the potential importance of coffee heirlooms in the Arabian Peninsula as a source of genetic diversity, there is limited information available on these genotypes. This information is essential for the development of new coffee varieties that can better adapt to changing environmental conditions, such as pest and disease pressure, abiotic stress, and changes in consumer preferences [7]. Furthermore, since over 60% of wild coffee species are in danger of extinction due to accelerated environmental change, gathering complete information, and characterizing this germplasm is of utmost importance [8]

Another issue facing the coffee industry as it struggles to cope with an over-supplied market, is adulteration. It has long been known that coffee is often adulterated with less expensive and readily available plant material [9]. Coffee adulteration has become a more serious issue for the industry in recent years due to the significant expansion in the variety of coffee recipes, stores, and ultimately consumers [10]. Therefore, developing molecular means like genetic barcodes to identify and authenticate the varieties can help mitigate the problem.

In Saudi Arabia and Yemen, C. arabica has been cultivated for at least four centuries on the terraced slopes and narrow valleys of the western mountains at different altitudes ranging mostly from 1200 to 2000m above sea level [11,12]. Most of what is grown now in southwestern Saudi Arabia are old cultivars that have been around for hundreds of years [12]. It is likely that these diverse populations are a result of successive introductions of genetic material from Eastern Ethiopia by Arab traders over centuries of uninterrupted exchange across the narrow strait of Bab El-Mandeb [13]. Therefore, it is safe to assume that the southwestern corner of the Arabian Peninsula contains the most genetic diversity of C. arabica outside the species' center of origin in the Ethiopian highlands [6]. Unfortunately, these genetic resources have not attracted enough attention from the scientific community, except for the 1989 FAO expedition to southern Yemen led by A. Eskes [14] and two recent studies by Tounekti et al. [12] and Montagnon et al. [6]. These three studies reported the existence of considerable diversity among coffee populations in the Arabian Peninsula. It is worth noticing that the present coffee populations have evolved over hundreds of years in a semi-arid environment [15] marked by recurring droughts, uneven distribution of rainfall, heat stress and high irradiance. Therefore, it is expected that these genotypes could be the source of interesting genes that confer stress tolerance [16].

The molecular phylogenies of coffee species have been established using variations in intergenic spacer sequences [17-19] and introns [20] of plastid DNA, as well as internal transcribed spacer (ITS) sequences of rDNA [21] and a combination of four plastid regions and ITS [22]. Chloroplast DNA (cpDNA) sequence variation is widely used in systematics and for making phylogenetic inferences at different taxonomic levels. Introns and intergenic spacers are known to exhibit high rates of mutation [23]. The trnT-trnL and trnL-trnF intergenic spacers, the trnL intron, and the atpB-rbcL intergenic region are useful for evolutionary studies at low taxonomic levels, and these regions have been extensively used in phylogenetic studies to analyze cytoplasmic polymorphism and to study the demographic history of several species [23,24,25,26,27].

Overall, further research is necessary to fully comprehend the diversity and potential of diploid and tetraploid coffee species and to utilize this information to develop new coffee varieties that can better meet the needs of farmers and consumers in the future. The present study aims to estimate the genetic diversity of local coffee populations in Saudi Arabia and to examine their genetic relatedness using chloroplast intergenic spacer markers.

2. Results

The successful amplification of all four intergenic spacer barcode sequences (atpB-rbcl, TrnT-TrnL, TrnL-trnF, TrnL) was achieved, resulting in a single band of the expected size (see supplementary Figure). The respective sequences for each barcode were submitted to the National Center for Biotechnology Information (NCBI) via Bankit submission, and accessions numbers for each barcode sequence are presented in Table 2. This study utilized the Coffea arabica (MN894552.1) chloroplast genome sequence, available on NCBI, as a reference genome sequence to detect polymorphisms among DNA sequences in 56 Coffea arabica accessions across all four barcoding primers.

The percentage of polymorphic sites for each sequence was determined by dividing the number of variable nucleotides by the length of the entire region and multiplying the result by 100. Table 3 presents the number of nucleotide sites (NNS), variable polymorphic sites (VPS), number of segregating sites (NSS), number of haplotypes (NH), nucleotide diversity (ND), and average number of nucleotide differences (ANND) for each barcode primer and the cumulative results for all four primers. The combined sequences showed the highest NSS, followed by the atpB-rbcl primer, while the TrnL primer had the lowest NNS. The TrnT-TrnL primer had the highest VPS (341), followed by atpB-rbcl, while the lowest (154) was recorded for TrnL-TrnF. The atpB-rbcl had the highest ND, followed by TrnT-TrnL (0.051), with TrnL-TrnF showing the lowest ND. Additionally, the atpB-rbcl had the highest ANND (185.54), while TrnL-TrnF exhibited the lowest value (25.23) for ANND.

NNS=Number of nucleotide sites, VPS= variable polymorphic sites, NSS= number of segregating sites, NH= number of haplotypes, ND= Nucleotide diversity, ANND= average number of nucleotide difference.

The nucleotide base composition of each barcode primer was determined and is presented in Table 4. The average nucleotide base composition of atpB-rbcl was recorded as 33.15% T(U), 16.60% C, 34.49% A, and 15.76% G. For TrnC, the composition was 26.7% T(U), 15.9% C, 37.6% A, and 19.8% G. TrnT-TrnL had a composition of 39.18% T(U), 13.88% C, 33.84% A, and 13.10% G (Table 4). The singleton variable sites (STVS) and parsimony information sites (PIS) for each chloroplast barcode are presented in Table 5. The TrnT-TrnL barcode recorded the highest number of STVS (338), followed by TrnL-TrnF (133), TrnL (52), and atpB-rbcl (1) had the lowest. Similarly, for PIS, TrnL had the highest number (155), followed by atpB-rbcl (137), while the minimum was recorded for TrnL-TrnF (45) across all sequences (Table 5). A phylogenetic analysis was constructed based on concatenated sequences of all four barcodes using the maximum likelihood method and Kimura 2-parameters model. The percentage of trees in which the associated taxa clustered together is shown next to the branches. This analysis involved 56 nucleotide sequences, and the final dataset comprised a total of 4114 positions. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The final phylogenetic tree divided the 56 accessions into five different groups.

Figure 1. Evolutionary analysis by Maximum Likelihood method using four barcodes with 1000 bootstraps constructed in MEGA 10.0 using the concatenated sequence of atpB-rbcl, TrnLc, TrnT-TrnL and TrnL-TrnF.

Figure 1. Evolutionary analysis by Maximum Likelihood method using four barcodes with 1000 bootstraps constructed in MEGA 10.0 using the concatenated sequence of atpB-rbcl, TrnLc, TrnT-TrnL and TrnL-TrnF.

The phylogenetic analysis grouped the coffee accessions into four different groups. The first group contained six accessions (KSA42, KSA29, KSA2R, KSA41, KSA43, KSA11R, and KSA38) that were mostly from the Rayda district of Assir region. The second group contained seven accessions (KSA51R, KSA3R, KSA27, KSA60, KSA4R, KSA7R, and KSA1R), all from Jazan Region except KSA60 was from Assir. The third group contained 12 accessions (KSA45R, KSA13R, KSA39, KSA25, KSA35, KSA59, KSA52, KSA36, KSA24, KSA22, KSA37 and KSA46), all collected from Jazan Region except KSA59 from the north of Assir Region. The fourth and largest group contained 43 accessions that can be further subdivided into four subgroups. The first subgroup (IVa) was a diverse one and contained 12 accessions originating from the three regions. Subgroup IVb contained three accessions (KSA33, KSA28 and KSA5R) all from Jazan Region. Subgroup IVc contained seven accessions, six from Jazan and one from Al-Baha. Subgroup IVd contained 11 accessions, eight from Jazan region, two from Assir and one from Al-Baha.

3. Discussion:

The genetic diversity present in crop wild or primitive relatives (CWR) plays a crucial role in the effectiveness of crop improvement programs. These wild or unknown genotypes exit in diverse habitats, many of which are currently facing significant threats due to habitat degradation and climate change [8,14]. It is estimated that approximately 60% of wild coffee species are at risk of extinction worldwide; similarly, underutilized old varieties are disappearing from the orchards. This emphasizes the urgent need to conserve these species through both in situ and ex situ measures, which can help preserve their genetic diversity for future utilization.

While morphological descriptors are commonly used to characterize different coffee species, molecular markers are considered more efficient in distinguishing closely related species and cultivars [23], and more precise and reliable than morphological and biochemical markers [24]. Furthermore, several studies have demonstrated that specific regions of the chloroplast genome can serve as DNA barcodes for a wide variety of plant species [25,26]. Selection of suitable plastid genomes offer sufficient genetic information for distinguishing between genotypes. Additionally, when choosing suitable DNA barcoding loci, the variable regions should be given a primary consideration [27]. Therefore, the objective of this study was to evaluate the evolutionary and phylogenetic relationships among fifty-six Arabica coffee species by utilizing four DNA barcoding markers (atpB-rbcl, TrnL-TrnF, TrnL, and TrnT-TrnL). This research aimed to investigate the potential of four DNA barcode loci (specifically, atpB-rbcL, TrnL, TrnL-TrnF, and trnL-trnL from the chloroplast region) for the identification and provision of phylogenetic information on local Arabica coffee species. All four regions were successfully amplified using universal primers, yielding clear and reliable results. However, earlier studies have indicated that there were cases of partial amplification from the respective barcode loci's universal primers [28,29].

Similarly, other studies [30,31] have indicated that additional barcode primers, including matK, rbcL, and trnL-trnF, have demonstrated successful amplification within coffee species, which aligns with the findings of the present study. Despite the abundance of available data on DNA barcoding of angiosperms, there is currently limited information regarding specific barcodes that can guarantee an accurate species identification in all cases [32]. Often, a barcode that performs effectively for one group of plants may prove inadequate for another group, especially in the case of recently diverged species [33]. The current study successfully identified all fifty-six accessions as Coffea arabica, showcasing the effectiveness of the universal DNA barcode primers. Likewise, multiple studies have extensively documented the reliability of matK and rbcL, either individually or in combination, as DNA barcodes that can be used with confidence across various plant species [34]. Additional reports have recommended the utilization of rbcL as a valuable DNA barcode locus, primarily due to its relatively compact length of 500 bp, robust universal primer, high success rate in PCR amplification, and excellent sequencing quality [35]. However, other DNA barcodes, such as TrnL-TrnF and the TrnL spacer, have also been suggested as reliable alternative barcodes for identification of species [36]. The extent of sequence variation among the species or terminals under analysis is a crucial factor in determining the effectiveness of any barcoding locus [34]. Other reports have recommended the utilization of rbcL as a valuable DNA barcode locus due to its relatively compact length of 500 bp, robust universal primer, high success rate in PCR amplification, and excellent sequencing quality [25], along with other DNA barcodes such as TrnL-TrnF and the TrnL spacer [36]. The degree of sequence variation present among the species or terminals being analyzed is a determining factor in the efficacy of any barcoding loci [34].

The number of singleton variable sites was found to be higher in TrnLc, TrnL, and TrnL compared to atpB-rbcl. Similarly, TrnL and atpB-rbcl had more parsimony information sites than the rbcL barcode spacer region. These findings are consistent with a previous study conducted by [37], which reported that TrnL-trnF and matK barcodes exhibited greater variability than rbcl in Indian coffee arabica genotypes. The present study also found similar results for PIC among the four barcodes analyzed. Similarly, previous research has indicated that TrnL-TrnF and matK loci exhibit greater sequence polymorphism than rbcL, as suggested by [38,39], and [40]. The current study's results support these findings. Hence, the present study found that all four barcode sequences, which were evaluated as candidate barcode markers, met the DNA barcoding criteria outlined by [33]. Specifically, these markers exhibited sufficient sequence variability to enable effective discrimination among Saudi coffee genotypes. Similarly, [41] previously suggested that indels involving multiple residues may serve as useful diagnostic markers for species discrimination and phylogenetic analyses [42].

The phylogenetic analysis grouped the Saudi Coffea arabica genotypes into four groups with a clear influence of geographic origin suggesting the genotypes of each region share one or more common ancestor (Figure 2). For instance, accessions KSA11R, KSA41, KSA42 and KSA43 from the isolated Rayda district of Assir region were grouped in clusters I and II. The accessions representing very old trees (KSA36, KSA44, KSA46, KSA47) segregated in the middle of the phylogenic tree in groups III and IVa. Similar results were reported by [37] where the grouping using single and multilocus barcode primers was strongly influenced by the geographic origin of the genotypes. A molecular analysis of coffee arabica genotypes from Saudi Arabia using SRAP markers grouped them into five distinct groups based mostly on their origin [23]. The accessions collected from Jazan region primarily clustered in groups II and IV, whereas those from Al-Baha and Assir regions constituted a different group. Similar surveys of genetic diversity among coffee populations in northern Yemen [6] and southern Yemen [14] found that each district (valley) have its own cultivars. Another study using genotyping by sequencing (GBS) showed that genetic closeness correlated with geographic proximity [43].The current study provides further evidence to support this finding. It was also suggested chloroplast sequences provide high level effective and accurate information for plant evolution and high level preservation [44]. For future studies on this economically significant crop, we recommend using resequencing and genome-wide association studies (GWAS) to discover additional polymorphic markers associated with important agro-morphological traits. These markers would be beneficial for a range of investigations in Coffea crop. Ultimately, the polymorphic markers established and confirmed in this research hold potential as a valuable genomic asset for molecular breeding, identification, and biogeography studies of Arabica coffee.

4. Materials and Methods

4.1. Plant Material

A survey was carried out at several sites in the Sarawat mountain range, running parallel to the Red Sea from the southeast to the northwest through the three administrative regions of Jazan, Assir, and Al-Baha. The survey covered a strip of terraced mountains located between latitudes 17˚N and 20˚N, the most northern location where coffee is commercially grown in the world. The coffee gardens included in the survey were found at altitudes ranging from 900 to 2000 m a.s.l. In total, we collected young leaves from 56 accessions, from Jebel Fayfa (Fayfa district), Eddayer, Maadi (Haroub district), Jebel Al-Gahr (Al-Rayth district), Rayda valley (Assouda district in Assir region), Mahayel Assir district, Al-Majarda district and Jebel Shada (Al-Mekhwah district of Al-Baha region) (Table 1). We tagged and sampled 1-3 trees representing each tree population. Each accession was given a code starting with the acronym “KSA” (e.g., KSA-1), but, for the sake of simplicity, we dropped the acronym in the figures. The letter “R” was added to the code of accessions 1-19, 45, and 51 to indicate that they were sourced from a small, local coffee germplasm collection established in the Fayfa district.

4.2. DNA extraction

Plant material, consisting of young leaves from various C. arabica accessions, was collected from representative trees in each population and transported to the lab in a cooler. The leaves were sanitized by immersing them in a 5% sodium hypochlorite solution for 1-2 minutes and then rinsing them with sterile distilled water. The material was then ground in liquid nitrogen and stored in an -80ºC freezer. DNA was extracted from 100 mg of mixed powder using an innuPREP Plant DNA Kit (Analytik Jena), following the manufacturer's protocol. DNA quality and concentration were determined using a Nanodrop ND-1000 spectrophotometer (Saveen Werner, Sweden).

4.3. Chloroplastic DNA amplification and sequencing

Four chloroplast DNA regions were considered (Table 2). PCR was performed in a 25 μl volume containing 2 μl of template DNA, 10 μl of 1X innuMix Standard PCR, and 1 μM of each primer (Table 2). The Gene Amp PCR System 9700 was used with the following program: initial denaturation at 94°C for 5 min, 35 cycles of denaturation at 94°C for 1 min, annealing at 49-52°C for 60-75s, and elongation at 72°C for 60-75s, followed by a final polymerization at 72°C for 10 min (Table 2). To check the effectiveness of PCR, positive control using sterile water was included in all amplifications. The PCR products were checked by electrophoresis on 1% agarose gel in TAE buffer, and DNA was visualized under UV light after staining with ethidium bromide.

Table 2. General information about the PCR primers used in this study.

Sr#	Sequence 5′-3′	Target	PCR condition	Source
1	CATTACAAATGCGATGCTCT	trnT-trnL	Hybridation : 50°C/1 min Elongation : 72°C/1min	[45]
	TCTACCGATTTCGCCATATC
2	CGAAATCGGTAGACGCTACG	TrnL	Hybridation : 49°C/1.15 min Elongation : 72°C/1.15 min	[45]
	GGGGATAGAGGGACTTGAAC
3	GGTTCAAGTCCCTCTATCCC	TrnL-trnF	Hybridation :52°C/1 min Elongation : 72°C/1min	[45]
	ATTTGAACTGGTGACACGAG
4	GAAGTAGTAGGATTGATTCTC	atpB-rbcL	Hybridation : 50°C/1 min Elongation : 72°C/1min	[46]
	TACAGTTGTCCATGTACCAG

Table 2. General information about the PCR primers used in this study.

Sr#	Sequence 5′-3′	Target	PCR condition	Source
1	CATTACAAATGCGATGCTCT	trnT-trnL	Hybridation : 50°C/1 min Elongation : 72°C/1min	[45]
	TCTACCGATTTCGCCATATC
2	CGAAATCGGTAGACGCTACG	TrnL	Hybridation : 49°C/1.15 min Elongation : 72°C/1.15 min	[45]
	GGGGATAGAGGGACTTGAAC
3	GGTTCAAGTCCCTCTATCCC	TrnL-trnF	Hybridation :52°C/1 min Elongation : 72°C/1min	[45]
	ATTTGAACTGGTGACACGAG
4	GAAGTAGTAGGATTGATTCTC	atpB-rbcL	Hybridation : 50°C/1 min Elongation : 72°C/1min	[46]
	TACAGTTGTCCATGTACCAG

The amplified products were purified using the GFX PCR kit (GE Healthcare). Sequencing reactions were carried out by Congenic using Sanger technology, separately for each strand to obtain independent forward and reverse sequences. The forward and reverse fragments were aligned, and additional reactions were conducted in case of any discrepancies.

4.4. Sequence Analysis:

The scanner software-2 was utilized to determine the quality of the sequences. The four barcode samples of each coffee arabica genotype were manually curated and aligned using the contig assembly program in Bio Edit 7.0 software to ensure high-quality sequences. Nucleotide sequences obtained from the 57 accessions were initially aligned using CLUSTAL W [47] and analyzed with MEGA program version X. The number of individuals, number of nucleotide sites, variable polymorphic sites, number of segregating sites, number of haplotypes, nucleotide diversity, and average number of nucleotide differences of each barcode marker and consensus sequence were measured using DNAsp (v6) [48]. The quantification of insertion events in the sequence was determined by the number of variable sites where the addition of one or more nucleotides signals polymorphism. Likewise, the number of deletions was determined by the variable sites where polymorphism arises due to the removal of one or more nucleotides. The identification of the number of transitions in the sequences was based on the number of variable sites where polymorphism occurred due to the exchange between two purines (A and G) or two pyrimidines (C and T). On the other hand, the number of transversions was determined by the variable sites where polymorphism resulted from the replacement of a purine with a pyrimidine (Table 4). To determine the number of mutation events that have occurred in a sequence, the sum of variable sites and the number of distinct mutations observed at the same nucleotide site across different samples are combined. This quantification considers both different types of polymorphisms and multiple occurrences of mutations within the sequence. Various parameters were estimated for each sequence region to differentiate them, based on the number of monomorphic or polymorphic sites, the number of parsimony informative sites (PIC), nucleotide diversity (π), haplotype diversity (Hd), and the total number of mutations [49], singleton variable site (STVC) [50] (Table 5).

4.5. Evolutionary analysis by Maximum Likelihood method:

The evolutionary history was inferred by using the Maximum Likelihood method and Kimura 2-parameter model [51]. The tree with the highest log likelihood (-22360.57) is shown. The percentage of trees in which the associated taxa clustered together is shown next to the branches. Initial tree(s) for the heuristic search were obtained automatically by applying Neighbor-Join and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood (MCL) approach, and then selecting the topology with superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. The proportion of sites where at least 1 unambiguous base is present in at least 1 sequence for each descendent clade is shown next to each internal node in the tree. This analysis involved 57 nucleotide sequences. There was a total of 5381 positions in the final dataset. Evolutionary analyses were conducted in MEGA X [52].

5. Conclusion:

To summarize, this study utilized a DNA barcoding approach to investigate the molecular relationships among fifty-six Arabica coffee samples collected from the southern region of Saudi Arabia. The three-barcode region, namely TrnT-TrnL, TrnL-TrnF, and TrnL, exhibited higher sequence variability compared to the atpB-rbcl barcode region and effectively differentiated the local coffee genotypes by the presence of unique variable sites (singletons and parsimony). Moreover, the combination of DNA sequences from these barcode loci, analyzed through phylogenetic analysis using the maximum likelihood method, grouped similar coffee genotypes, providing improved resolution, and understanding of the studied genotypes. These findings will contribute to future research programs in the identification and conservation of Arabica coffee using DNA barcoding markers.