There are around 1200 native cultivars, 3000 synonyms [8], and currently there are 5 published *Olea europaea* assemblies comprising 3 different cultivars, *Olea europaea* cv. 'Farga' (Oe6/Oe9) [17,18], *Olea europaea* cv. 'Picual' (Oleur0.6.1) [13] and *Olea europaea* cv. 'Arbequina' (Oe\_Rao) [19], and a purported wild variety *Olea europaea ssp. sylvestris* (Oe451) [20]. Due to the prohibitive cost and effort involved in whole-genome sequencing, it was once an endeavour only available to consortia and solely focused on model organisms. However, the rapidly evolving sequencing technologies from companies such as Illumina, Pacific Biosystems, and Oxford Nanopore have consistently surpassed Moors law [21], allowing even small labs to afford to sequence their own species of interest. This affordable access to rapid high-throughput sequencing has been revolutionary for many areas of biology, with applications in de novo genome assembly, genotyping, gene–trait associations, metagenomics, transcriptomics, and epigenetics [22,23]. Given the economic and cultural importance of olive, it is not surprising that we have seen so many olive genome-sequencing projects in a short space of time. Each of the available olive genome assemblies uses different sequencing technologies, or a combination of methods. It is necessary when choosing a sequencing approach to consider the application of the genome, as the different types of errors, error rates, and biases that can come from a particular methodology will affect the overall quality and completeness of the assembly [24]. Furthermore, the different choices of assembly tool and pipeline will impact on the contiguity, accuracy, and handling of repeat regions in highly polyploid species [25].

Genetic profiling in many crops, as well as the analysis of genetic variation within and between their populations, has been achieved using cheap and effective biochemical and DNA markers such as random amplified fragment polymorphic DNA (RAPDs) [26,27], amplified fragment length polymorphisms (AFLPs) [28,29], SSRs [10,15,30,31], and singlenucleotide polymorphisms (SNPs) [32–34]. A sequence assembly of any quality is not required for many projects such as genetic profiling, establishing genetic relatedness, QTL mapping, or to perform a GWAS, although access to a complete and correctly annotated genome assembly provides the best account of individual genome variation and provides more information, increasing the potential resolution when using methods capable of

recovering a higher density of variants, such as genotype by sequencing (GBS) [35], or whole-genome sequencing (WGS).

Both GBS and WGS could be used to call SNPs but differ massively in terms of missing data and their cost effectiveness. The main difference is that GBS is a reducedrepresentation approach to sequencing, that while quick and cost-effective, results in much more missing data due to its DNA fragmentation step. GBS uses restriction-site-specific digestion enzymes to fragment DNA samples. The DNA fragments then have unique barcode sequences ligated to the ends of the DNA fragments before fragment size selection is performed. The primary advantage of GBS is that by assigning sample-specific barcoded adapters, it is possible perform multiplexed sequencing in a single Illumina flow-cell lane for a large number of samples [35], making it much more cost effective than WGS. The number of SNPs that can be identified from within a WGS dataset can be significantly more than with GBS; however, this level of resolution is not always necessary in geneticlinkage-based research [36]. GBS has already been used in olives to generate genetic maps [20,37,38], study the diversity, and perform association analysis [39]. While the GBS library preparation itself is relatively simple, demultiplexing of the raw data is required to process the samples. This step can add extra difficulty for any researcher not already familiar with data processing.

It is likely that as the cost of sequencing and data processing continues to fall, we will see even more cultivar genomes assembled using differing technologies and assembly pipelines. It then becomes important to be able to assess the quality and functionality of a genome in order to choose the right assembly for a project's goal. Several tools and methods already exist to estimate genome completeness and contiguity [40]; however, the potential impact and bias that may arise from a reference genome's genetic background is not fully understood, nor to what extent this issue may affect different types of analysis.

In the context of a typical population genetics study, we wanted to understand if the genetic background would have a significant impact on the interpretation of group clustering and genetic relationships. Given the difference in data production and cost, we asked how comparable GBS datasets are to WGS datasets for this mode of research. Furthermore, we assessed the practicality of combining WGS and GBS datasets (Table 1), as this would allow for unrelated sequencing projects to access and combine publicly available data using two very different genotyping methods.

**Table 1.** List of *Olea europaea* cultivars sequenced using GBS and available WGS data.

#### **2. Results and Discussion**

*2.1. Genome Assembly Comparison*

In this study, we used five publicly available *Olea europaea* subsp. *europaea* assemblies (Table 2) comprising three different cultivars, 'Farga' (Oe6/Oe9) [17,18], 'Picual' (Oleur0.6.1) [13], and 'Arbequina' (Oe\_Rao) [19], and a purported wild accession, *Olea europaea* subps. *sylvestris* (Oe451) [20]. Our full list of GBS samples and WGS samples is listed in Table 1. We removed 12 GBS samples from the initial 36-sample dataset that either failed to pass quality controls or for which no WGS data were available. The analysis was carried out on the remaining 24. The selection process is described in more detail in Section 4.5. These samples likely failed during the library preparation, although, there could be many reasons behind this problem. One of the most common is related to the quality of the DNA. Low DNA quality (e.g., due to impurities in the DNA extraction) reduces the efficiency of the restriction enzymes and leads to a partial digestion and reduction in the fragment population. Nevertheless, the low yield for the samples that failed (e.g., Barnea) indicates a bias in the amount of template used in the pool. Because the libraries were performed by an external service, it is difficult to assess where the problem was, but our guess is that it was related to the amplification step of some samples during the library preparation [35]. This would be an issue for a typical study of olive cultivar relatedness and would therefore need to be repeated. Fortunately, in this case the remaining samples were sufficient to compare GBS performance against WGS and the impact of assembly bias.

The assembly of GCA\_002742605.1 (Oe451) used a whole-genome shotgun sequencing approach with the Illumina HiSeq 2000 platform to sequence *Olea europaea* var. *sylvestris* (wild olive). The assembly was performed using SOAPdenovo to generate a genome coverage of 220.0× [20]. For Oleur0.6.1, *Olea europaea* cultivar 'Picual' was selected as the genetic background [13]. The assembly process integrated Illumina HiSeq2500 and PacBio RSII sequencing to improve gap filling. *Olea europaea* cultivar 'Farga' based GCA\_900603015.1/Oe6 [17] was assembled from Illumina HiSeq2500 with a genome coverage of 380×. Oe9 is an updated version of Oe6 [18], where a genetic map was used to anchor scaffolds to chromosomes. The most recently published is the *Olea europaea* cultivar 'Arbequina' (GWHAOPM00000000/Oe\_Rao; NGDC) [19], which used Oxford Nanopore long-read sequencing and Hi-C data to construct chromosomes. The quality and

completeness of each assembly was assessed before read mapping, SNP calling, and further analysis. Our comparison of genome assemblies was based on four metrics, (1) contiguity stats, (2) gene space completeness, (3) a k-mer completeness assessment using Merqury, and (4) the LTR Assembly Index (LTR) that evaluates the completeness of the genome using LTR retrotransposon elements (Table 2).

**Table 2.** Genome assembly and contiguity statistics.

Assembly stats generated using custom script. Gene space completeness assessed using BUSCO v5 with eudicot\_db10 dataset. Completeness, repeat regions, and sequencing data incorporation evaluated with Merqury v1.3 and Illumina short-read-derived k-mers. LTR Assembly Index (LAI) was estimated using the LTR\_Retriever tool v2.9.0 with the default parameters.

> With respect to (1), assembly contiguity statistics were generated using a custom script (see Section 4.4). First, we note that total assembly size varied between genetic backgrounds; Oleur0.6.1, at 1.68 Gb long, was by far the largest assembly, and over 360 Mbs longer than both Oe6 and Oe9 (Table 2). Indeed, it was over 500 Mbps longer than Oe451 or the 'Arbequina' genome (Oe\_Rao). Jiménez-Ruiz et al. proposed that the bigger size of the 'Picual'-variety genome could be explained by the presence of a large number of duplicated DNA fragments coming from a very recent partial genome duplication event or artificially introduced repeat regions from assembly issues [18]. This comparison is interesting to observe because we are comparing the shotgun approach (Oe451, Oe6 and Oe9, Oleur0.6.1), where the genome is fragmented into reads of 250–800 bp [20], to sequencing with Pacbio first-generation RSII long reads (Oleur0.6.1), and the Oxford Nanopore third-generation long-read approach where reads are commonly 10–30 kb [41]. Shotgun approaches are very cost effective, but the process of sequence reassembly is more complicated than with long reads. The short reads make it difficult to correctly assemble repeat regions, particularly tandem repeats. This can also cause issues when using a genetic map to anchor

to chromosomes, as the shorter read may not span a large enough section of the sequence to contain genetic markers. The trade-off is that on a per read basis, short reads are currently still more accurate and cheaper than long reads [41]. However, when only using short reads for an assembly, greater coverage is required, and this increases the number of errors in the dataset. Error correction and the filtering of low-quality reads is therefore an important step in the genome. We noted that Oe451, according to the Merqury assessment, had the highest QV score even with a 220× coverage, indicating careful management of the raw data during assembly. It is important to remember that each assembly will have its own unique set of errors introduced by sequencing or assembly issues, and thus any assembly should not be considered as a definitive sequence but rather, as stated in [42], only a working hypothesis.

Next, we checked (2) gene space completeness using Benchmarking Universal Single-Copy Orthologs (BUSCO) [43].  Highly fragmented assemblies may contain greater percentages of missing or fragmented genes. The number of duplicated orthologs may be used as an indicator of possible errors that can occur when an assembly tool mistakenly assembles, or fails to assemble, reads into repeat regions. Nevertheless, polyploidy events and partial duplications may also lead to an increase in the number of duplicated genes detected by BUSCO, so often a deeper analysis needs to be carried out in order to distinguish a biological cause from a technical problem. In the publications for each of the olive tree genomes, BUSCO was used to evaluate assembly quality; however, the results were not comparable in terms of gene set database or BUSCO software version used. Thus, we ran BUSCO v.5 on all assemblies, using the eudicot\_db10 gene set, which contains 2326 eudicot-specific genes. The results are different from those in the publications of the genomes used; however, this can be explained by the fact that the gene dataset's specificity can vary greatly depending on which is used, and older versions of the gene datasets may be less complete. Using the eudicot\_db10 dataset, we obtained an identical completeness score of 96.6% for Oleur0.6.1, Oe6, and Oe9, which indicates that a high proportion of the core gene space was captured by these assemblies [44].

In addition to a recent genome duplication event in 'Picual', there are signs of at least two older whole-genome duplications in Olive [20], so it is unsurprising to see high levels of duplicated genes here. However, the 'Picual' genome has 51.5% duplicated genes, more than twice any other assembly. This could be an indication of a technical error such as a failure in the consensus calling due to high levels of heterozygosity, or in this case, there is also evidence of a biological origin for a recent genome duplication. BUSCO gene sets cannot have taken into account very recently uncovered duplications, and so this high duplication percentage in Oleur0.6.1 is difficult to evaluate. The number of missing BUSCO genes was low in the 'Picual' and 'Farga' assemblies (1.4–1.5%), as was the percentage of fragmented genes (1.9–2%), indicating an overall high level of completeness. Oe451 scored the lowest using this dataset, with a completeness score of only 85.9, an indication that a large portion of the genome is still unassembled despite its pseudo chromosomes; this may be an issue for use in particular studies, such as synteny, or the discovery of candidate genes by QTLs, but may not be an issue for SNP calling as we had intended if the GBS sites are not in the unassembled regions. Oe\_Rao scored 93.4 for completeness and had a lower percentage of duplicated genes compared to the 'Farga' and 'Picual' assemblies.

We assessed each of the assemblies using (3) Merqury, a k-mer-based method able to evaluate the quality and completeness of a de novo assembly without the need for a reference genome. Merqury uses a similar method to KAT in which high-quality sequencing reads are decomposed into k-mers datasets, then the k-mer sets are compared to the genome assembly. Merqury summarises the quality assessment with two values: a completeness score that measures the completeness of the assembly based in the k-mer populations of the assembled and the unassembled reads, and a phred-scaled consensus quality (QV) score that measures the error produced during the haploid sequence consensus calling of the assembly. Additionally, Merqury's copy number spectra plots (Supplemental Figure S1) allow for a visual inspection of unassembled reads and artificial duplications [45]. The evaluation of Oe\_Rao could not be carried out correctly as Merqury requires short reads and the Oxford nanopore reads used in the assembly are not compatible with the process, and no Illumina reads were publicly available at the time of this analysis. The two 'Farga' assemblies and the 'Picual' assembly had an identical BUSCO gene completeness score, yet, using Merqury, the overall completeness of the 'Farga' and 'Picual' assemblies was vastly different. This indicates a high amount of reads which were never used in the final 'Farga' assemblies (Oe6/Oe9). Indeed, in the Supplemental Figure S1 plots Oe6 and Oe9, we can see that this is the case. This result may be partially explained by the co-sequencing of the fungus genome *Aureobasisium pullulans* with the olive samples, which led to an assembly of 18 Mb [17], but the K-mer multiplicity indicates that there was also an important proportion of the repetitive content unassembled. However, these three scored an equally high QV, indicating highly accurate consensus calling during contig assembly. Discounting the Oe\_Rao results, Oe451 scored the lowest in terms of completeness and had the most k-mers found only in the assembly, but had the highest QV score. Sequencing error types and error rates vary with the sequencing technology [36,41,46]; these may cause contig misassembles and scaffolding, and thus affect the reliability of a genome assembly for use in the development of genotyping markers, breeding programs, or population studies. A good example of the impact of the genome completeness on olive genetic studies can be found in Kaya et al.: 51% and 75% of the SNPs with strong association signals were mapped to the Oe451 and Oe6 reference genomes, respectively. Although these percentages are not correlated with the estimation of the completeness from Table 2, they are a good indication that the quality of the genome can influence the usability in other analyses.

Finally, the assemblies were evaluated using the LTR Assembly Index (LAI). Genome assemblies based on short-read sequencing technologies such as Oe451, Oe6, and Oe9 presented lower values than the assemblies based on long-read sequencing technologies, Oleur0.6.1 and Oe\_Rao, as was expected and previously described in the use of LAI for the assessment of genome quality [47]. It is interesting to note that the updated version of Oe6, Oe9, has a lower LAI (4.34 compared with 5.10), meaning that some of the transposable element genome space was lost or fragmented during the improvement. From all the genome assemblies used, only the Oleur0.6.1 had a LAI value over 10, such that it is considered a standard value for high-quality genomes with a good contiguity.

Because every sequencing project used a different technological approach, genetic background, and assembly methodology, there is necessarily a great deal of difference between them and their overall quality. While the limitations of short read length can impact the handling of repeat regions during scaffolding with a single short-read mapping at the incorrect or multiple regions, long-read technologies such as Oxford Nanopore and Pacbio have much higher error rates. Further, the choice of assembly tool, corrections, and the use of an alignment tool with and without available references can impact the quality of a genome assembly. Considering the information collected in Table 2, Oleur0.6.1 appears to be the most complete assembly with the highest QV score of the four we could test. However, it remains unclear if the high percentage of duplicated genes is an error or true genome duplication. Furthermore, this genome is not yet anchored to chromosomes, limiting its use in some studies.

#### *2.2. Effect of Reference Genome Choice on Population Analysis*

We explored the potential effect genome selection could have on a typical genetic diversity study such as a population analysis by mapping the same GBS reads to each of the available reference genomes, and performed the same population analysis on all the resulting SNP datasets. This analysis began with a quality control assessement of our dataset, which consisted of unprocessed reads from all 36 of our selected Olea europaea cultivars, amounting to over 326 million raw reads. After pre-processing, mapping, and SNP calling (see Section 4.5 and Supplemental Figure S19) using three different genomes (Sylvestris/Oe451, Farga/Oe6, and Picual/Oleu r0.6.1), samples that had failed at one or more quality control checkpoints on all three genomes were removed from the rest of the analysis (see Sections 4 and 4.5) (Supplemental Table S1). Removed samples are listed in Table 1. As mentioned previously, the GBS samples may have failed during library preparation, or sequencing and would normally need to be repeated; however, 24 samples are sufficient for the focus of this work. The remaining 24 samples were then additionally mapped to the two remaining genomes (Farga/Oe9 and Arbequina/Oe\_Rao), before variant calling, SNP filtering, and population analysis.

#### 2.2.1. Analysis GBS Read Mapping and Variant Calling

Considering only the remaining 24 high-quality samples, the percentage of mapped reads did not vary much between assemblies; indeed, there was only a 2% difference separating the highest and lowest mapping genomes (Table 3). The average number of sites for each sample also showed very little variation, except for Oleur0.6.1, which had around 100,000 more sites, likely due to the increased number of duplicated genes in this assembly. It is important to highlight that a substantial difference was found in the number of variants called. Both Oe451 and Oe6 had double the number of variants called compared to Oe\_Rao, and almost six times that of Oleur0.6.1. However, by far the most SNPs identified were from the improved 'Farga' genome, Oe9 with 23.7 million variants (Table 3). In terms of variants per loci, Oleur0.6.1 and Oe9 had, on average, almost triple the number of Oe451, Oe\_Rao, and Oe6. The total number of SNPs called after filtering (see Section 4.5) was similar for all genomes. This massive drop off from Oe9 was most significantly explained by the 10,000 Kb thinning and 10,000 minimum quality score filtering steps. These two alone removed ~97% of total SNPs. Oleur0.6.1 was notable for having been left with half as many SNPs as the 'Sylvestris' and 'Farga' genomes. Perhaps the Oleur0.6.1 assembly contained several repeat regions which were incorrectly collapsed or, alternatively, reads may have been incorrectly assigned as a repeat(s), both of which could increase or decrease the number and frequency of SNPs, as multiple alleles from the same locus might be mistakenly identified as having come from different loci or vice versa. As can be seen in Table 2, there is indeed variation in the number of duplicated BUSCO genes and the k-mer duplications found across all of these genomes, indicating a difference in the assembly of repeat regions of each of the tested genomes. In such cases, it might be expected to also see variation in the levels of heterozygosity. However, heterozygosity per site was not significantly different, ranging from 0.31 to 0.36 (Table 3), making it difficult to interpret the source of this phenomenon. To test if collapsed regions were indeed the cause of this SNP count variation, we extracted the allele frequency of SNPs from three random samples ('Grappolo', 'Manzanilla de Sevilla', 'Piñonera') called from three of the genomes (Oe451, Oe6, Oleur0.6.1) and compared the frequency of alleles between 0.25 and 0.75 (Supplemental Figure S2). Allele frequencies for an individual different from 0, 0.5, and 1, such as 0.25 and 0.75, may be an indicative of the collapse of four copies into one during the genome assembly. Once adjusted by percentage, we observed similar profiles, with most alleles tending towards 0.25 and a small peak around 0.5. However, from our sub-sampling we did not see any significant difference and so it remains unclear why using Oleur0.6.1 resulted in half the number of SNPs and Oe451.

**Table 3.** GBS and WGS read mapping and variant calling.

**Table 3.** *Cont.*

#### 2.2.2. GBS Population Structure

The filtered SNPs obtained from each genome mapping were analysed using FAST-STRUCTURE (DARTR), ADMIXTURE (LEA), discriminate analysis of principle components (DAPC), and principal component analysis (PCA) to identify the distinct relationships among the cultivars (Supplemental Figures S3–S17).

In two previous studies, there were two distinct clusters identified among cultivars and a third cluster containing wild relatives [13,18]. In this work, the results from FAST-STRUCTURE identify two groups (K = 2, with cross entropy errors ranging 0.7 to 0.8) as being the most probable clustering for all GBS and WGS datasets (Supplemental Figure S8). As wild olive samples were not included in the analysis, two groups are in agreement with the previously published data [13]. When comparing the results of the model-based clustering from ADMIXTURE and genetic-distance-based PCA clustering, it was noted to be in agreement also with FASTSTRUCTURE. The first three principal components of the PCA only explained ~30% of the variation. However, it was sufficient to see a clear separation of the different clusters, a result observed constantly throughout all our SNP datasets (Supplemental Figures S3 and S4). In all cases, the first principal component (PC1) (~14%) separated the 24 cultivars into a group of 13 that could be described as primarily composed of eastern Mediterranean cultivars ('Barri', 'Abou Kanani', 'Abou Satl Mohazam', 'Majhol-1013 , 'Temprano', 'Verdial de Velez-Malaga-1 , 'Uslu', 'Kalamon', 'Morrut', 'Mari', 'Picual', 'Manzanilla de Sevilla', 'Abbadi Abou Gabra-842 ) and a group of 11 mostly northern Mediterranean cultivars ('Mastoidis', 'Klon-14-1812 , 'Grappolo', 'Mavreya', 'Piñonera', 'Leccino', 'Myrtolia', 'Menya', 'Manzanillera de Huercal Overa', 'Koroneiki', and 'Arbequina'), reflecting the complex history of olive domestication [11,48], and aligning with previously published results [13,18]. ADMIXTURE at K = 2 produced a similar grouping as the PCA (Supplemental Figures S5–S7). The levels of estimated admixture in each individual cultivar did not appear to vary depending on the assembly used. It is also interesting to note the near identical results regardless of the number of SNPs. Oleur0.6.1 had almost half the number of filtered SNPs as Oe451 or Oe9, but it appears that this still provides sufficient resolution to evaluate genetic relatedness using these analyses.

In general, the DAPC posterior membership was also estimated to be the same regardless of assembly used during mapping. We could see some variation in group membership when samples were compared by country of origin. There was largely agreement while comparing all genomes; however, we saw some discrepancies, such as in the Oe6, Oe9, and Oe\_Rao genomes, which gave a near identical result in which the cultivar 'Koroneiki' of Greece showed a much higher probability of belonging to the Syrian genetic group. Indeed, this was more in agreement with the comparison of principal components 1 and 2, seen in all datasets belonging to the Syrian genetic group.

We calculated, for each dataset, genetic distances and constructed neighbor-joining distance trees using the R package POPPR (Supplemental Table S4). We can see in all cases there are two clear operational taxonomic units (OTUs), or clusters, defined in each tree (Figure 1) with a smaller third cluster of cultivars, which is more variable depending on the reference genome used. Samples in the third group were those which had higher levels of admixture (Supplemental Figures S5 and S6) and were more difficult to assign. However, the two large groups were in agreement with the eastern Mediterranean and northern Mediterranean clusters identified in the PCA (Supplemental Figure S3). The topology of each grouping varied, but only slightly, and only within an out. There was no placement of individuals in another cluster; even those with greater admixture were either in a seperate group or with cultivars of the same region. The topography of the trees produced from Oe451 and Oe6 were near identical, and even closer matched than Oe6 to its updated Oe9 version. The bootstrapping values were high in all cases with only 2–3 nodes falling below 50 on any tree. Syrian and Italian cultivars were the most consistent in terms of clustering when using any other assemblies. The cultivar 'Mari' was seen in all trees to cluster with Syrian samples, but, interestingly, the DAPC results indicate this same grouping only when using a 'Farga' or 'Arbequina' genome. This might be due to some genetic relationship between the genomes and the 'Mari' sample, such as a possible introgression with some of the SNPs.

**Figure 1.** Neighbour-joining trees constructed using GBS data. Comparison of the neighbour-joining (NJ) trees constructed from SNPs called using GBS data mapped to each of the reference genomes. Using the R package Poppr, the genetic distance between each sample using Provesti's distance and NJ trees is plotted with a bootstrapping support of 100. The genome used to generate each is labelled above the corresponding tree.

Although there was great divergence in the number of SNPs associated with the use of each assembly, we saw little effect on the outcome of the population analysis. As observed in maize, reference genome selection can impact results [49,50]. Gage (2019) reported a reduction in the ability to robustly identify key loci of interest in a genome-wide association study (GWAS). Missing loci and inaccurate structural variations introduced during genome assembly may introduce reference bias, which may be more problematic in such cases. The impact of using any of the five selected assemblies in their current states appears to have been minimal in this population analysis, and indeed each reference genome provides confirmation of the results. However, it is still important to consider the vast differences in variants identified from each assembly and how they could impact other types of analysis. Using more than one reference assembly in this way can remove much of the bias introduced by many of the choices made during a genome sequencing project (see Section 2.1) and is something that could be repeated in future olive breeding projects to improve the reliability of identified structural variations.

#### *2.3. Analysis of WGS Data*

As the population analysis using GBS data showed that the selection of different reference genomes can make a difference in the phylogenetic relationship of one sample to another, but not in the overall grouping of samples, and no significant difference was seen in the levels of estimated admixture, we wanted to see if this was also the case with WGS data. Since WGS would likely generate many more variants and sites, any small differences may be diluted or potentially exaggerated. To process the raw whole-genome sequencing (WGS) read set, we followed the same steps as used for the raw GBS read set (see Section 4.5).

The processed reads were mapped to two assemblies, the purported wild type (Oe451) and the genome that scored the highest overall in our analysis (Oleur0.6.1). First, looking at the number of reads mapped to each, we saw that the WGS read set had a much higher number of mapped reads than the GBS set, as expected (Table 3). In both cases, the percentage of processed reads that was mapped was over 100%. Some double mapping was detected but not enough to cause issue. Variant calling using WGS data produced significantly more variants than GBS; 144 million from WGS-Oe451 and 128 million called from WGS-Oleur0.6.1 before filtering. It was not surprising to have so many more variants called using WGS reads because of the lower amount of missing data [36], but it was interesting to note that when we compared the total number of variants called from Oe451 and Oleur0.6.1 with the GBS data, there was a difference of 3.3 million between the two datasets before filtering. The difference between GBS-Oe451 and GBS-Oleur0.6.1 was almost double after filtering, but there was only a ~10% difference between WGS-Oe451 and WGS-Oleur0.6.1.

We next ran the same population analysis script as previously used for the GBS datasets. The FASTSTRUCTURE, PCA, admixture, and DAPC results were again very similar regardless of the genetic background of the reference genome. As before, there was some difference in the topography of trees but overall, the same OTUs were reproduced. Bootstrapping values were more consistently higher using Oleur0.6.1. The results are near identical to the GBS results, suggesting that it was possible to merge GBS and WGS datasets.

Here, we also compared the genomes that produced the highest and lowest number of SNPs. It is worth noting that, again, despite the higher number of SNPs coming from Oe451, there appeared to be no difference in the outcome of the analysis when using any of the reference genomes or SNP datasets. The reasons for this might be that SNPs provide limited information and so a certain threshold, in terms of reliability and number, is needed for an SNP dataset to be useful. The quality of the SNPs called and filtered from each reference genome was sufficient, even with the lower number of SNPs from Oleur.0.6.1, to extract the same level of detail. We also must consider that while five reference sequences were used, the genetic background of one was a purported wild type and all the others were Spanish cultivars. Our results show no difference when using either sequence, but this may change if a cultivar from another region with a more distinct genetic history is used.

#### *2.4. Analysis of WGS/GBS Merged Data*

The ability of an SNP dataset to differentiate between groups can be quite low because of their biallelic nature [51] and so it could be expected that increasing the number of SNPs from WGS should significantly improve the power of differentiation. However, our results show that GBS and WGS data performed almost identically, despite a difference of roughly 10 times the number of SNPs called using WGS data (Table 3). Given the similarity of the population analysis results, we wanted to know if it was possible to combine variants from GBS and WGS datasets. This would allow for more collaborative opportunities and the supplementation of smaller datasets with already available public data.

The same two genomes (Oe451 and Oleur0.6.1) were chosen for this experiment. The vcf files of GBS-Oe451 and WGS-Oe45 were merged together using bcftools and a list of shared sites (see Section 4), as were GBS-Oleur0.6.1 and WGS-Oleur0.6.1. The merged WGS/GBS-Oe451 file consisted of over 150 million variants and after filtering, 9537 biallelic SNPs, less than 1% of the total variants. The WGS/GBS-Oleur0.6.1 merged file contained only around 11 million variants and less than 1% of total variants were high-quality SNPs. This is in stark contrast to the difference seen with GBS-Oe451 and GBS-Oleur0.6.1. We selected only the SNPs from shared sites and did not allow for any missing data, so this would explain the massive drop off compared with either WGS or GBS after only filtering.

Again, we ran the same population analysis script as with all previous SNP sets. FASTSTRUCTURE found the mostly likely number of clusters, based on lowest crossentropy, to be K = 18–19 for WGS/GBS-Oe451 and K = 20 for WGS/GBS-Oleur0.6.1 (Supplemental Figure S8). This may be an overfit of samples due to there being two of each genotype (Supplemental Figure S11) and the WGS samples being more like their GBS counterpart than anything else, creating the highest probability clusters. However, we could see that the change in cross entropy decreased between possible clusters pastK=2 and K = 3. The sample clustering by PCA was near identical for the WGS/GBS-Oleur0.6.1 and WGS/GBS-Oe451 datasets (Figure 2). Importantly, they were similar to the GBS and WGS only analysis with all the reference genomes. We saw that GBS and WGS samples closely grouped together, forming the same clusters seen previously. The only notable exception was GBS Arbequina. It seems likely that this sample was mislabelled at some point in the process, as it is the only sample that does not pair up with its WGS counterpart. DAPC analysis was analogous to all previous results, with the same groups being identified. Analysis by NJ trees (Figure 3) revealed no major differences in terms of OTUs, but we observed that Oleur0.6.1 had some very low bootstrapping values and Oe451 had the better bootstrapping values, more closely resembling the WGS-only tree.

**Figure 2.** PCA clustering of WGS/GBS datasets. This shows the clustering of WGS and GBS merged samples. Principle components 1 and 2 are shown for each. (**a**) PCA of merged data with reference genome Oe451; (**b**) PCA of merged data with reference genome Oleur0.6.1.

**Figure 3.** Neighbour-joining trees constructed from merged WGS and GBS data mapped to the Oleur0.6.1 and Oe451 reference genomes. Genetic distance between each sample using Provesti's distance and NJ trees plotted with a bootstrapping support of 100.

#### **3. Conclusions**

The results of this work indicate that GBS and WGS sequencing data are highly comparable for population structure analysis. Certainly, more information can be obtained from WGS, and great advances have been made to reduce its cost, but GBS remains much more affordable. For small sampling sizes WGS maybe be preferable, but a large-scale genotyping project can quickly become too expensive for most labs; we show that GBS is a cost-effective alternative capable of providing near identical results and identifying the same genetic relationships at a much lower cost. Further, the ability to successfully combine the two sequencing methods opens opportunities to mine data from a wider range of sources. This would also be a significant cost-saving approach to consider for expanding existing olive tree WGS datasets and past genotyping projects, and allow for collaborations between research groups that have used either WGS or GBS genotyping methods. Although the sequencing cost is decreasing, it is not feasible to re-sequence hundreds of accessions to study the structure of the population or to identify some unknown accessions.

With respect to our comparative analysis between the different reference genomes used, the Oleur0.6.1 sequence may be the most accurately assembled, but lacking chromosomes can limit its application. However, we could not properly test the 'Arbequina' assembly using Mercury's best practices and therefore accurately compare it to the others in all aspects of our quality assessment. Any differences between the current existing genome assemblies are so small that we can say that for SNP-based genetic profiling at least, all are certainly suitable. Thus, the availability of genetic maps and chromosomes could potentially be a more important consideration when choosing which to use in future

projects. Furthermore, we show that while the genetic background of the reference genome plays a role in the possible number of variants identified, there was little overall effect on population structure, genetic clustering, or analysis of genetic relationships. As more reference genomes of existing cultivars are rendered to the chromosome level, it might well be necessary to compare performance before choosing which to use. The current assemblies are of a purported wild type and three Spanish cultivars, so it is still possible that the genetic background may yet be an important factor if using cultivars from other regions as the reference.

Our results highlight the advantages of GBS while at the same time bringing to the table the possible limitations. GBS is a commonly used technique for other crops, but it has not been routinely implemented in olive breeding. Traditional breeding needs these types of tools to accelerate the development of new varieties able to face the important challenges olive production is currently facing.

#### **4. Materials and Methods**

#### *4.1. Plant Material and Genomic DNA Extraction*

In accordance with Belaj et al. (2012), and maximising the genetic diversity in a reduced number of genotypes, 36 olive tree cultivars were selected from the World Olive Germplasm Collection of the Andalusian Institute of Agricultural and Fisheries Research and Training (WOGBC) located in Córdoba Spain (Table 1). Total genomic DNA was extracted from fresh leaves using the Illustra DNA extraction kit Phytopure GE Healthcare (UK) in accordance with the protocol described in the manufacturer's instructions. To ensure high-quality DNA was used for sequencing, the purity was measured with a Qubit 2.0 Fluorometer Life Technologies NY (USA). DNA concentration was then normalised to 20 ng/μL. The DNA had a minimum 260/280 ratio of 1.8.

#### *4.2. GBS Library Construction and Sequencing*

The DNA were sent to DISMED, and the libraries were prepared by BGI. DNA samples were digested with ApeKI (New England Biolabs, Ipswitch MA) for 2 h at 75 ◦C, and T4 ligase (New England Biolabs, Ipswitch MA) was used to ligate to sticky ends, one of 36 unique "barcodes", and the "common adapter". Samples were incubated at 22 ◦C for 1 h and heated at 65 ◦C for 30 min. A set of 36 digested DNA samples, each with a different barcode adapter, were obtained. A total of 7 μL of each component of this set was combined in a unique sample and purified in a final volume of 50 μL with a commercial kit (QIAquick PCR Purification Kit; Qiagen group (Germany), according to the manufacturer's instructions. The result of this library was amplified in 50 μL, containing 5 μL of pool DNA fragments, 1× Taq Master Mix New England Biolabs (UK), and 12.5 pmol each of PCR primers, the sequences of which were: PCR primer: 5 AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT and PCR primer: 5'-CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCT GAACCGCTCTTCCGATCT, containing complementary sequences for amplifying the fragments of DNA with ligated adapters. The PCR conditions were a primer step of 5 min at 72 ◦C; 98 ◦C for 30 s; 25 cycles of 98 ◦C for 30 s, 65 ◦C for 30 s, and 72 ◦C for 30 s; and a final extension step at 72 ◦C for 5 min. The library was purified as above (in a final elution of 30 μL) and 1 μL was used for the quality evaluation of fragment sizes. The library was considered suitable for sequencing if adapter dimers were minimal (~128 pb in length) and the majority of the other DNA fragments were between 170 and 350 bp. Paired-end sequencing of one 48-plex library per channel was performed on an Illumina HiSeq2000 Analyzer by BGI Genomes.

The sequences of the barcode adapter were: 5-ACACTCTTTCCCTACACGACGCTCT TCCGATCTxxxx and 5-CWGyyyyAGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT, where "xxxx" and "yyyy" are the barcode and barcode complement, respectively. The second, or "common" adapter sequence was shared among all samples and consisted of an

#### ApeKI-compatible sticky end: 5´-CWGAGATCGGAAGAGCGGTTCAGCAGGAATGCCGAG and 5´-CTCGGCATTCCTGCTGAACCGCTCTTCCGATCT.

#### *4.3. WGS Library Construction and Sequencing*

Raw sequencing data were obtained directly from the authors of Jiménez-Ruiz et al. (2020). Re-sequencing of all varieties was performed by 2 × 150 paired-end sequencing with Illumina HiSeq 4000, and all sequencing was carried out at the Duke Center for Genomics and Computational Biology (Durham, NC, USA). Raw data are available at NCBI BioProject ID: PRJNA556567.

#### *4.4. Sequence Assembly Assessment*

The five available assemblies used were *Olea europaea* subsp. *sylvestris* (version GCA\_002742605.1/Oe451; NCBI), *Olea europaea* cultivar Picual (version Oleur0.6.1; https: //genomaolivar.dipujaen.es/db/ (accessed on 2 February 2020)), the *Olea europaea* cultivar Farga (version GCA\_900603015.1/Oe6; NCBI) [17] andits updated version (GCA\_902713445.1/Oe9; NCBI) [18], and *Olea europaea* cultivar Arberquina (GWHAOPM00000000/Oe\_Rao; NGDC) [19]. The same three forms of quality assessment were carried out on all the assemblies: contiguity, gene space completeness, and a k-mer-based evaluation.

Assembly and contiguity were calculated using a custom script; FastaSeqStats (https:// github.com/aubombarely/GenoToolBox/blob/master/SeqTools/FastaSeqStats (accessed on 10 June 2013). N50 is a commonly used marker of sequence contiguity, counting the number of bases in the shortest fragments needed to span 50% of the genome. N90 is the same information at 90% of the total genome. Average sequence length provides similar information as the N50. L50 is the smallest number of contigs needed to cover 50% of the genome; L90 provides the same information at 90%.

Benchmarking Universal Single-Copy Orthologs (BUSCO) genes were used to evaluate the completeness of the assembly by searching for a list of known ancestrally conserved genes within the reference genome. All assemblies were assessed using BUSCO v.5 with the eudicot\_db10 gene set containing 2326 eudicot-specific genes.

Merqury v1.3, (https://github.com/marbl/merqury (accessed on 10 January 2021) a k-mer-based method, was used to evaluate both quality and completeness. Using a similar method to KAT, high-quality sequencing reads were decomposed into k-mers datasets, then the k-mer sets were compared to the genome assembly. Merqury generates a k-mer dataset from the Illumina short-read sequencing data used in the genome sequence assembly. By decomposing the original sequencing reads into k-mers, Merqury can count how many times each k-mer appears in the assembly, as well as k-mers from the original Illumina reads not incorporated at all or that appear only in the assembly. These k-mer data were used to generate a completeness score and a phred-scaled consensus quality (QV) score, along with copy number spectra plots to visually inspect the assembly for unassembled reads and artificial duplications [45]. The Illumina short-read datasets were either downloaded from the same location as the assembly or kindly sent by the genome assembly curator. In the case of the 'Arbequina' assembly, only long-read sequencing was performed, which was unsuitable for use with Merqury.

LTR Assembly Index (LAI) was estimated using the LTR\_Retriever tool v2.9.0 with the default parameters [47].

#### *4.5. Read Processing, Mapping, Filtering, and Variant Calling*

Raw reads were first processed and then mapped only to Oe451 to assess quality, using the pipeline described below and in Supplemental Material Figure S10. Samples with low SNP count can affect the overall number of SNPs available for analysis due to the percentage of missing data encountered during the VCF filtering (see Section 4 and Supplemental Figure S10). Before removing poorly performing samples, the full process was repeated with two other assemblies (Oleur0.6.1 and Oe6) to confirm the results (Supplemental Table S1). Total SNPs called from Oe451, Oe6, and Oleur0.6.1 after

filtering were initially in the range 300–900. To increase this number further, low-quality samples were removed. The choice of samples to be included or excluded was based on performance at different stages of processing and the amount of missing data. To achieve this, the number of reads in each sample was counted and there were three clear groups, those with only 50,000 or fewer reads, a second group with 600,000 or fewer, and the majority of samples were within the range of 2–20 million reads. Samples with 50,000 or fewer reads were removed from further analysis, as this is too low to yield any useful information, along with any samples having less than 50% of the average number of sites per sample.

This reduced the sample set from 36 to the current 24, but significantly increased the number of total SNP loci shared across all samples. Sequence read mapping was then performed on all five genome assemblies to evaluate the influence that a particular assembly has on SNP discovery and on a population analysis using the following steps.

After sequencing, the GBS raw reads were demultiplexed by BGI Genomics. Next, the raw reads were processed with FASTQ\_MCF v1.05 [52] to remove Illumina adapters and reads with a phred-scaled quality score of less than 30 and/or shorter than 50 bases. After trimming, the paired reads were aligned to each of the five genome assemblies; Oe451, Oe6, Oe9, Oleur0.6.1, and Oe\_Rao, using BWA (Burrows–Wheeler Alignment Tool) v0.7.17-r1188t with default parameters. Prior to mapping, each reference genome was indexed using BWA. The output file of the mapping was in an unsorted SAM format, and these were converted to bam and sorted to save space and increase SNP calling efficiency using SAMTOOLS v1.7 [53]. Once sorted, the bam files were merged into a single bam file with BAMADDRG (https://github.com/ekg/bamaddrg (accessed on 14 April 2018).

Variants were called with FREEBAYES v1.3.1-16-g85d7bfc [54] using a custom script; MultiThreadFreeBayes, (https://github.com/aubombarely/GenoToolBox/tree/master/ SNPTools/MultiThreadFreeBayes (accessed on 15 May 2018) this script allows FREEBAYES to run faster by using multiple threads on several scaffolds at the same time. Finally, the variant file (VCF) was filtered with VCFTOOLS v0.1.15 [55] using the following parameters: retain only biallelic SNPs, remove indels, a minimum read depth of 5 with a minimum mean depth of 20, a minimum SNP quality of 1000, no missing data, and an MAF of 0.05. Finally, SNPs were thinned to 1 per 10,000 Kb.

To merge the GBS and WGS datasets, we used the BCFTOOLS v1.7 [53] merge function. The WGS VCF file contained many more sites than the GBS, so a bed file containing all sites from the GBS dataset was supplied using the argument regions-file. This was carried out to reduce the final file size, as any site with missing data would eventually be removed during the filtering steps. After the files were merged, they were filtered using the same parameters as above.

#### *4.6. Population Analysis*

Each SNP dataset was analysed in the same way, with the same script to estimate population structure and genetic diversity from each of the datasets. The R script Olea\_pop.R with commands and notations is available at https://github.com/frieljames/Olive\_WGS\_ GBS (accessed on 1 April 2021). Supplemental Figure S19 shows an overview of the R programming pipeline. As a summary of the script, we cross analysed population structure in our datasets with the use of R and three different methods: (1) STRUCTURE [56], a Bayesian-based clustering method assuming Hardy–Weinberg equilibrium and linkage equilibrium between loci within populations, (2) discriminant analysis of principal components (DAPC) [57], which uses sequential K-means and model selection to infer genetic clusters, and (3) principal components analysis (PCA), a method that uses genetic distance to infer clusters. To begin, we used FASTSTRUCTURE [58], part of the DARTR package, a faster and more resource-efficient method of running STRUCTURE. The LEA package [59] performed a STRUCTURE analysis on the SNP data to infer the mostly likely genetic clusters based on allele frequency and clustering probability. ADMIXUTRE was estimated by running LEA's snmf functional analysis for the estimation of ancestral populations

(K) of 1–20. This estimates the admixture coefficient in the selected K range using sparse Non-Negative Matrix Factorisation algorithms. The output of this was used to estimate the admixture of groups at the largest K value. To perform a DAPC and a PCA, we used the R package ADGENETv.1.3-1 [60]. DAPC is different to STRUCTURE in that it does not rely on model assumptions or prior information; instead, it uses a multivariate method with sequential K-means and model selection to infer genetic clusters and assign individuals to clusters. To ensure that maximum variance was being used while attempting to avoid an overfit of the data, DAPC was performed using the suggested optimal number of principle components (PCs) for each dataset (Supplemental Figures S9–S11). These were identified as the number of PCs with the highest a-score, predicted by the optim.a.score() function, which selects an evenly distributed number of PCs in a pre-defined range, computes an a-score for each, and then interpolates the results using splines. For the PCA, the package ADGENET took a genlight object of SNP data to generate a genetic distance matrix for use in the PCA clustering.

Further estimation of possible ancestral populations was achieved with a neighbourjoining (NJ) approach, again based on allele frequencies. This was carried out with the R package POPPR [61]. The POPPR function aboot() allowed for the construction of a dendrogram with 100 bootstrapping support and obtained genetic distance using bitwise.dist, a method that calculates the fraction of different sites between samples equivalent to Provesti's distance. The tree was visualised in R with the package APE [62] using the function plot.phylo(). Afterwards, the tree was exported in the newick format to generate figures using FigTree v1.4.4 [63].

Finally, genetic diversity was calculated using the gl.basic.stats () function with the R package DARTR [64]. Selection and neutrality were estimated by calculating the number of segregating sites, nucleotide diversity, Watterson's theta, and Tajima's D. These values were generated using POPGENOME [65]. Additionally, population diversity was estimated using expected and observed heterozygosity, Fst and Fis values, for each of the assigned groups (POPGENOME) (Table S4).

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/plants10112514/s1, Figure S1: Copy number spectrum plots (spectra-cn plot) of k-mer multiplicity from each set of x50 coverage Illumina reads. Histograms coloured by number of times each k-mer in the read set is found in a given assembly. The result is a set of histograms relating kmer counts in the read set to their associated counts in the assembly, Figure S2: Histograms of allele frequency sampled from 'Grappolo', 'Manzanilla de Sevilla', 'Piñonera' cultivars using vcftools, Figure S3: PCA clustering of GBS SNP sets. Plots show PCs 1 and 2 for SNPs called from a given reference genome, Figure S4: PCA clustering of WGS SNP sets. Plots show PCs 1 and 2 for SNPs called from a given reference genome, Figure S5: Estimated admixture for GBS SNP datasets Oe451, Oe6, Oe9, & Oleur.0.6.1. Ancestry co-efficient (0.0-1.0) estimated using STRUCTURE assigned populations, Figure S6: Estimated admixture for GBS Oe\_Rao SNP dataset and WGS Oe451 & Oleur.0.6.1 SNP datasets. Ancestry co-efficient (0.0-1.0) estimated using STRUCTURE assigned populations, Figure S7: Estimated admixture for merged WGS/GBS Oe451 & Oleur.0.6.1 SNP datasets. Ancestry co-efficient (0.0-1.0) estimated using STRUCTURE assigned populations, Figure S8: Estimated number of populations using STRUCTURE for each dataset. Cross-entropy values along the y-axis, Figure S9: DAPC preliminary optimization of GBS data prior to running DAPC analysis. Reference genome SNP set indicated to the left of each row. Left most column contains the results the function optim.a.score() used to indicate best number of principal components to include. Center column the value of Bayesian Information Criterion (BIC) to assess the best supported model. Right-hand column contains the preassigned group (country of origin) representation using selected PCs and BIC model, Figure S10: DAPC preliminary optimization of WGS data prior to running DAPC analysis. Reference genome SNP set indicated to the left of each row. Left most column contains the results the function optim.a.score() used to indicate best number of principal components to include. Center column the value of Bayesian Information Criterion (BIC) to assess the best supported model. Right-hand column contains the preassigned group (country of origin) representation using selected PCs and BIC model, Figure S11: DAPC preliminary optimization of merged WGS and GBS data prior to running DAPC analysis. Reference genome SNP set indicated to the left of each row. Left most column contains

the results the function optim.a.score() used to indicate best number of principal components to include. Center column the value of Bayesian Information Criterion (BIC) to assess the best supported model. Right-hand column contains the preassigned group (country of origin) representation using selected PCs and BIC model, Figure S12: DAPC group composition and clustering of GBS data set. Reference genome SNP set indicated to the left of each row. Left most column contains the membership probability of each individual sample using discriminant functions, proximal to admixture coefficients used by STRUCTURE. Right-hand column clustering using predefined groups (country of origin) with selected model, Figure S13: DAPC group composition and clustering of WGS data set. Reference genome SNP set indicated to the left of each row. Left most column contains the membership probability of each individual sample using discriminant functions, proximal to admixture coefficients used by STRUCTURE. Right-hand column clustering using predefined groups (country of origin) with selected model, Figure S14: DAPC group composition and clustering of merged WGS and GBS data set. Reference genome SNP set indicated to the left of each row. Left most column contains the membership probability of each individual sample using discriminant functions, proximal to admixture coefficients used by STRUCTURE. Right-hand column clustering using predefined groups (country of origin) with selected model, Figure S15: DAPC posterior membership probability of GBS data set. Reference genome SNP set indicated to the left of each row. Membership probability of each individual at K=2, Figure S16: DAPC posterior membership probability of WGS data set. Reference genome SNP set indicated to the left of each row. Membership probability of each individual at K=2, Figure S17: DAPC posterior membership probability of merged WGS and GBS data set. Reference genome SNP set indicated to the left of each row. Membership probability of each individual at K=2, Figure S18: Genotype-by-sequencing library preparation pipeline, Figure S19: Sequencing read processing and population analysis pipeline. Table S1: GBS read processing and variant calling, Table S2: GBS read processing and variant calling Stats for the raw GBS sequencing reads mapped to all 5 genome assemblies. Low quality individuals have been removed improving the overall number of SNPs available, Table S3: WGS read processing and variant calling. Stats for the raw WGS sequencing reads mapped to two genomes. Only samples matching the remaining GBS samples were selected, Table S4: Results of Popgenome analysis genetic diversity and neutrality.

**Author Contributions:** Conceptualisation, A.M.F.-O., A.B. and F.L.; Methodology, A.M.F.-O., A.B., and J.F.; formal analysis, J.F.; writing—original draft preparation, J.F.; writing—review and editing, J.F., C.D.F., A.B., F.L. and A.M.F.-O.; supervision, A.B. and A.M.F.-O.; project administration, A.M.F.-O.; funding acquisition, A.M.F.-O.All authors have read and agreed to the published version of the manuscript.

**Funding:** H2020 project NUMBER 101000427—GEN4OLIVE: Mobilization of Olive GenRes through pre-breeding activities to face the future challenges and development of an intelligent interface to ensure a friendly information availability for end users. Additionally, the project number UJA 2011/13/13 Detección y aplicabilidad de marcadores SNPs para estudios de diversidad en variedades de olivo y su mapeo genético. PLAN PROPIO DE LA UNIVERSIDAD DE JAEN.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The novel GBS data produced for this study can be found under the NCBI Bioproject accession PRJNA750928. WGS data were already published, associated with the Bioproject accession PRJNA556567.

**Acknowledgments:** The authors want to acknowledge the technical and human support provided by CICT of Universidad de Jaén (UJA, MINECO, Junta de Andalucía, FEDER).

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Morphological and Molecular Characterization of Some Egyptian Six-Rowed Barley (***Hordeum vulgare* **L.)**

**Azza H. Mohamed 1,2,\*,†, Ahmad A. Omar 2,3,\*,†, Ahmed M. Attya 4, Mohamed M. A. Elashtokhy 5, Ehab M. Zayed <sup>6</sup> and Rehab M. Rizk <sup>7</sup>**

**Abstract:** Barley production is essential in Egypt. In the present study, 15 different six-rowed Egyptian barley cultivars were studied. To differentiate between the different cultivars under study in terms of morphological characteristics and ISSR, molecular characterization reactions were carried out. Moreover, four cultivars (Giza 123, Giza 126, Giza 136, and Giza 138) were selected for further studies using scanning electron microscopy (SEM). Computational analysis of the DNA barcoding sequences of the two plastid markers *rbc*L and *mat*K was executed, and the results were deposited in the NCBI database. The morphological traits showed low statistical significance among the different cultivars under study via the data collected from two seasons, suggesting that the mean field performance of these Egyptian cultivars may be equal under these conditions. The results showed that the phylogenetic tree was divided into four groups, one of which contained the most closely related genotypes in the genetic distance, including Giza 124, Giza 130, Giza 138, Giza 136, and Giza 137, which converge in the indicative uses of farmers. The seed coat of the studied cultivars was "rugose". The elevation folding of the rugose pattern ranged from 11 ± 1.73 μm (Giza 126) to 14.67 ± 2.43 μm (Giza 123), suggesting variation in seed quality and its uses in feed and the food industry. According to the similarity matrix of ISSR analysis, the highest similarity value (93%) was recorded between Giza 133 and Giza 132, as well as between Giza 2000 and Giza 126. On the other hand, the lowest similarity value (80%) was recorded between Giza 130 and (Giza 133 and Giza 132), indicating that these cultivars were distantly related. Polymorphism information content (PIC) ranged from 0.26 for the primer ISSR UBC 835 to 0.37 for the primers ISSR UBC 814 and ISSR UBC 840. The current study showed that the *mat*K gene is more mutable than the *rbc*L gene among the tested cultivars.

**Keywords:** plastid markers; DNA barcoding; ISSR markers; Egyptian barley; agro-morphological traits; cluster analysis; genetic variation; biplot

#### **1. Introduction**

Barley (*Hordeum vulgare* L.) is one of the main and oldest cereal crops on Earth. Worldwide, its grain production is ranked fourth after maize, rice, and wheat [1]. Barley is

**Citation:** Mohamed, A.H.; Omar, A.A.; Attya, A.M.; Elashtokhy, M.M.A.; Zayed, E.M.; Rizk, R.M. Morphological and Molecular Characterization of Some Egyptian Six-Rowed Barley (*Hordeum vulgare* L.). *Plants* **2021**, *10*, 2527. https:// doi.org/10.3390/plants10112527

Academic Editors: Milan S. Stankovic, Paula Baptista, Petronia Carillo and Josefina Bota Salort

Received: 22 September 2021 Accepted: 18 November 2021 Published: 20 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

generally considered a poor man's crop because it is easy to cultivate, with few requirements, and has a high capacity for adaptation to harsh environments. Some literature estimates the age of barley at 11,000 years [2]. However, six-rowed barley did not arise until after 6000 BC [3]. Archaeological evidence has dated barley cultivation to 5000–6000 BC in Egypt [4–7]. Barley products, especially bread and beer, comprised a complete diet in ancient Egypt. Based on the health benefits of barley, as well as the need for agricultural development and reducing wheat imports, Egypt is currently examining the return of barley to the bread-making industry as a 30% ingredient. The global production volume of barley reached 142.37 million metric tons in the 2017/2018 crop year [1]. Furthermore, it is expected that barley production will decrease to 140.6 million metric tons in the next crop year [8]. Egypt's barley production has fluctuated substantially in recent years as it increased over the past two decades, ending at 108,000 tons in 2019 [1]. With the new policies for the sustainable development plan and reclaimed land expansion, the area and productivity are expected to increase.

There are approximately 38 Egyptian barley cultivars—two- and six-rowed—but the six-rowed barleys are the most famous and widely used in Egypt. Field evaluations have shown many differences between genotypes [9,10]. Furthermore, [11] found that all studied traits showed significant differences between genotypes, environments, and interactions. Moreover, Sharma, et al. [12] used Euclidean distances based on non-hierarchical cluster analysis to categorize total accessions into diverse clusters, and to determine and select accessions with decent yield and performance for other ancillary traits. The candidate breeding lines can be used in hybridization for barley improvement programs.

Molecular markers are an essential tool used to directly detect the differences between and within genetic materials at the DNA level; they provide a robust estimate of genetic similarity that is not often obtained using morphological data alone [13]. A comparison can also be made to determine the genetic distance between the Egyptian cultivars based on field characteristics and molecular parameters. The inter simple sequence repeats (ISSRs) technique has been successfully applied to many crop species [14–16]. ISSRs demonstrate the specificity of microsatellite markers, and require no specific sequence information for primer synthesis, using the advantage of random markers [17]; thus, they have been widely used for cultivar identification in different crops [18–20]. Moreover, Guasmi et al. [17] found that ISSR primers exhibited variations in the percentage of polymorphism, resolving power (Rp), and band informativeness (Ib); the rate of polymorphism was 66.67%, the Rp ranged from 0.74 to 1.16, and the average Ib ranged from 0.24 to 0.39, suggesting that ISSRs are robust molecular markers that can distinguish between Egyptian cultivars. On the other hand, according to Drine, et al. [21], ISSRs and random amplified polymorphic DNA (RAPD) markers identified 72.2% and 61% of polymorphic bands, respectively. Several parameters were used to compare the relative efficiency of these marker systems, including the effective multiplex ratio (EMR), marker index (MI), and polymorphic information content (PIC); the ISSR system showed higher values for all of the parameters examined. Moreover, Wang, et al. [22] used 10 ISSR primers to investigate the variation between Tibetan and Middle Eastern barley genotypes; the Tibetan genotypes contained 91 allelic variants, of which 79 were polymorphic (86.81%), while the Middle Eastern genotypes contained 82 allelic variants, of which 66 were polymorphic (80.49%). These results suggest that ISSRs are robust molecular markers that can be employed to distinguish Egyptian cultivars.

The morphology of barley grains observed via scanning electron microscopy (SEM) showed starch granules smaller than the standard ones; they appeared in abnormal shapes with a conspicuous peripheral groove and sunken cheeks [23,24]. When SEM was used to study the detailed structure, it proved that starch was degraded by both pitting and surface erosion [25]. The apparent shape of the seed under the electron microscope indicates the quality of its industrial and agricultural importance.

DNA barcoding is a genetic identification technology that uses a genetic region of short DNA sequences, called the DNA barcode [26]; this can be reliably characterized by similar morphological characteristics and chemical compositions [27]; it has two main objectives: identifying organisms, where an unknown sequence matches a known species sequence, and exploring species that are similar in terms of habitat delimitation and description of species [28]. A short DNA sequence obtained from established target regions of the chloroplast genome can be used to classify genera and/or species of plants with respect to orthologous databases, compared to conventional PCR-based markers [29]. DNA barcoding has been proposed as an essential tool for resolving the significant gaps in our current understanding of biodiversity. Furthermore, Barley and Thomson [30] demonstrated that the success of DNA barcoding varies broadly across DNA substitution models, and has a substantial influence on the number of operational taxonomic identified units. Moreover, using recent advances in combinatorial pooling and next-generation sequencing, Lonardi, et al. [31] proposed a new sequencing approach that addresses the challenge of de novo selective genome sequencing in a highly efficient manner. Barcodes can be employed to explain the relationships between Egyptian cultivars, and their relation to sequences within the database.

The main objective of this study was to measure and characterize the differences between the most economically important Egyptian barley cultivars, especially in making bread. The study investigated 15 Egyptian six-rowed cultivars at the field level, the molecular level, via scanning electron microscopic examination, and via DNA barcoding. The results obtained from the present study will potentially enhance breeding programs and lead to the development of new adaptive or high-yield barley cultivars with specific improved traits.

#### **2. Results and Discussion**

#### *2.1. Field Experimental*

#### 2.1.1. Growing across Two Seasons

Figure 1 represents the average values of the field data during the two growing seasons, including grain filing period (day), maturity day (day), and hiding day (day) (Figure 1A); spike height (cm) and plant height (cm) (Figure 1B); the number of spikes per square meter and number of grains per spike (average of 10 spikes per square meter) (Figure 1C); biological yield (ton/ha) and grain yield (ton/ha) (Figure 1D); and weight of 1000 grains (g) (Figure 1E). The average mean values showed low statistical significance among the examined cultivars. There were no significant differences based on the least significant difference (LSD) for any of the studied traits except for biological weight, which showed a significant difference in values between the different genotypes (LSD = 1.78).

Figure 1 shows differences between genotypes in all studied traits, which were divided into three parts according to the convergence of the numerical values of the traits. There were indications of early cultivars being equal through the periods of maturity and seed fullness, and low statistical significance among them (Figure 1A). Plant height indicated vegetative solid growth (Figure 1B), which is sufficient for animal feed. These results are consistent with those of Amer, et al. [32], who found that the average yield of the new cultivar Giza 137 was 16.7 4.95 Ton/ha, while that of Giza 138 was 5.07 Ton/ha. These yields significantly exceeded the national checks Giza 123 and Giza 132 (3.88 Ton/ha). Giza 137 significantly out-yielded Giza 123 and Giza 132, by ~22.4 and ~20.7%, respectively. Furthermore, Giza 138 significantly exceeded the average of national checks Giza 123 (by ~25.6%) and Giza 132 (by ~23.9%).

On the other hand, Noaman, et al. [33] found that biological weight characterized new genotypes. Furthermore, Mariey, et al. [34] considered the Egyptian barley genotypes Giza 123, Giza 131, and Giza 136 to be salt tolerant. It is worthy of note that this study was performed under the conditions and climate of Giza Governorate, Egypt. Furthermore, the behavior of the varieties differs when studied under different environmental conditions such as in the Sinai Peninsula or on the northern coast—even though they have the same genetic background; the same can be said of their other behavior under harsh conditions.

**Figure 1.** The averages of morphological traits in 15 barley genotypes grown in two seasons—2017/2018 and 2018/2019: (**A**) grain filing period (day), maturity day (day), and hiding day (day); (**B**) spike height (cm) and plant height (cm); (**C**) number of spikes per square meter and number of grains per spike; (**D**) biological yield (ton/ha) and grain yield (ton/ha); and (**E**) weight of 1000 grains (g).

#### 2.1.2. Genetic Distance Dendrogram between Genotypes Based on Field Traits

The genetic tree of the genotypes was divided into four groups: Group I contained the most closely related genotypes in terms of genetic distance, including Giza 124, Giza 130, Giza 138, Giza 136, and Giza 1137 (Figure 2); members of this group were characterized by a high maturity day and high grain filing period (days), along with grain yield (Ton/ha) and biological yield (Ton/ha). Group II contained the two cultivars Giza 129 and Giza 133; this group could be described by a low number of grains per spike and high hiding days. Group III consisted of Giza 125 and Giza 2000 on one side of the group, and Giza 132 and Giza 134 on the other side (Figure 2); members of this group were characterized by high plant height and low-to-moderate maturity days. Group IV consisted of Giza 123, Giza 135, Giza 131, and Giza 126; this group could be characterized by height, a moderate number of spikes per square meter, spike height (cm), and the number of grains per spike, along with low-to-moderate weight of 1000 grain (g), hiding days, biological yield (Ton/ha), and maturity days (Figure 2). These results were consistent with the findings of Mariey and Khedr [35]. Moreover, based on their 10 agro-morphological traits, Mareiy, et al. [36] explored biplot and cluster analysis using Euclidean distance matrices and average linkage. According to PCA, all 15 genotypes fell into 4 groups. Cultivars in Group A tend to have higher yields, so they may be considered to be tolerant (Giza 16 and Giza 18). Nevertheless, Giza 124, Giza 132, and Giza 134 are among the cultivars in group D that produce lower grain yields. The characteristics of biological weight and biological yield are used to assess the production of grain in relation to the rest of the plant components, which are used as animal feed in the form of straw. Indeed, increasing the seed yield and decreasing the biological crop is beneficial to grain production, which is the goal of growing barley for nutrition and intensive production.

#### *2.2. Scanning Electron Microscopy (SEM)*

The term seed coat of barley caryopsis includes tissues from three separated organs: the pericarp, the testa, and the semipermeable membrane. Several unique compounds are synthesized in the seed coat, serving the plant's defense and control of its development in different ways. Additionally, many of these compounds are sources of industrial products and components for human consumption or animal feed [37]. The seed coat of the studied cultivars is "rugose" (Figure 3 and Table 1). The elevation folding of the rugose pattern ranges from 11 ± 1.73 μm (Giza 126) to 14.67 ± 2.43 μm (Giza 123). The extension of the rugose pattern (length) ranges from 16.00 ± 2.61 μm (Giza 126) to 18.67 ± 3.13 μm (Giza 136). The frequency pattern in 100 <sup>μ</sup>m2 ranges from 4.67 ± 0.51 (Giza 126) to 12.17 ± 1.69 (Giza 138). Thus, these cultivars could be promising for different purposes in service of contemporary Egyptian interests.

**Table 1.** The seed coat characteristics of four Egyptian six-rowed barley (*Hordeum vulgare* L.) cultivars (Giza 123, Giza 126, Giza 136, and Giza 138).

**Figure 2.** Cluster analysis and heatmap based on agro-morphological traits of 15 Egyptian six-rowed barley cultivars. The heatmap was constructed using JMP®, Version 15 (SAS Institute Inc., Cary, NC, USA, 1989–2019).

**Figure 3.** Scanning electron microscope (SEM) images of four Egyptian six-rowed barley (*Hordeum vulgare* L.) cultivars: (**A**) Giza 123, (**B**) Giza 126, (**C**) Giza 136, and (**D**) Giza 138. Scale bar = 100 μm. White arrows indicate the quality and shape of the rugose.

#### *2.3. Molecular Characterization and Genetic Relationships as Revealed by ISSR Markers*

The ability to effectively utilize genetic variability available to breeders is dependent upon an understanding of population diversity [38,39]. Thus, the primary benefit of cultivar differentiation at the molecular level is to explain with some accuracy the relationships between cultivars, in order to reduce selection costs within breeding programs and provide future breeders with molecular insights. The inter simple sequence repeats (ISSRs) fingerprinting profiles generated by 4 out of the 15 primers used in the present study, targeting 15 Egyptian six-rowed cultivars of barley, are displayed in Figure 4. The polymorphism generated by the 15 ISSR primers is summarized in Table 2. The 15 ISSR primers used in the present study produced a total number of bands (TNB) of 126, and 62 of those were polymorphic with uniqueness (PWU), with a polymorphism percentage (P%) of 50.07%. The TNB ranged from 5 for the ISSR UBC 844A and ISSR UBC 901 primers, to 14 for the ISSR UBC 835 primer. The number of PWU bands also varied, from two in ISSRs UBC 825 and UBC 901, to seven bands in the ISSR 857 and ISSR UBC 835 primers. The average number of PWU bands was 4.13 per primer (Table 2). The polymorphism information content (PIC) values varied between the ISSR primers. PIC ranged from 0.26 for the primer ISSR UBC 835 to 0.37 for the primers ISSR UBC 814 and ISSR UBC 840. Remarkably, some ISSR primers revealed distinct discrimination of 22–80% polymorphism, including ISSR UBC 825 and ISSR UBC 844A (Table 2). The ISSR primer UBC 844A recorded the lowest effective multiplex ratio (EMR) (6.20) and lowest marker index (MI) (0.02), whereas ISSR UBC 814 scored the highest in terms of PIC, resolving power (RP), and MI values (0.37, 12.67, and 0.05, respectively). Furthermore, ISSR UBC 835 scored the lowest values for PIC and MI (0.02 and 0.26, respectively) and the highest value in EMR (12.07).

**Figure 4.** ISSR–PCR product profiles of 15 investigated samples of *Hordeum vulgare L*: (**A**) primer UBC 814, (**B**) primer UBC 826, (**C**) primer UBC 840, (**D**) primer UBC 808, (**E**) primer 807, and (**F**) primer 851. M: molecular size marker (100 bp).

**Table 2.** ISSR marker profiles for 15 Egyptian six-rowed barley cultivars.

Each estimated parameter's minimum and maximum values are highlighted in yellow. MB: monomorphic bands; POU: polymorphic without uniqueness; UB: unique bands; PWU: polymorphic with uniqueness; TNB: total number of bands; P%: polymorphism (%); MBF: mean of band frequency; PIC: polymorphism information content; RP: resolving power; EMR: effective multiplex ratio; MI: marker index.

Genetic diversity in some six-rowed barley cultivars grown in Egypt was assessed using ISSR markers. The 15 ISSR primers produced 97 markers that were utilized to investigate the genetic diversity among the studied cultivars. A polymorphism percentage of 50.07%, with an average of 4.13 markers per primer, was found among the studied cultivars (Table 2). However, this number ranges from two for ISSR UBC 825 and ISSR 857, to seven for ISSR UBC 835. The ISSR primers produced single and unique bands, and four molecular primers had these bands (ISSR UBC 814, ISSR UBC 827, ISSR 807, and ISSR 851). The use of ISSR markers for fingerprinting previously resulted in high polymorphism between species, and reflected intraspecific variations within species [14,15,40]. In addition to the high level of polymorphism observed in the current study by ISSR, this may imply high insertional activity in the genome of the tested barley cultivars [21,41,42].

The genetic diversity parameter data revealed by ISSR markers were utilized to calculate the genetic diversity of the studied cultivars by using multivariate clustering, PCA, and heatmap analyses. In a PCA scatterplot, the ISSR markers reflect the robustness of the markers in categorizing the investigated cultivars. PCA analysis indicated that the four six-rowed Egyptian barley cultivars Giza 126, Giza 2000, Giza 125, and Giza 132 were distinct from the other cultivars (Figure 5). Neighboring affinity was also apparent between the Giza 135, Giza 136, and Giza 130 cultivars (Figure 5). Conversely, the rest of the cultivars—Giza 129, Giza 138, Giza 131, Giza 133, Giza 134, and Giza 137—were scattered at some distance from one another. The cultivars Giza 126 and Giza 2000 were the best foragers, as designated by cluster analysis (Figure 5), which also indicated a significant distance between Giza 123 and Giza 124 (Figure 5), and between Giza 132 and Giza 135, Giza 136, and Giza 130 (Figure 5). The differentiation of the studied cultivars in terms of years of release and pedigree may be due to previous alterations in production conditions. There is a possibility that these morphological characteristics can increase or decrease genetic variation between cultivars. Data from ISSR markers analyzed in this study might be explained by the instability of TNB insertion events, cultivar production, and behavior under environmental conditions [17,35]. There may be a correlation between the high degree of polymorphism observed in ISSR markers and genotype diversity [16,21,43]. Although there were differences between the dendrograms based on field characteristics, and in PCA results based on the molecular parameters, both sorted the cultivars into four groups closer to their uses in Egypt.

Multivariate compound similarity analysis is usually utilized to show more information about the genetic variance of plant breeds, which is detailed in heatmaps [40]. The multivariate compound similarities were presented as a heatmap constructed using R software. As indicated by the columns, 15 Egyptian barley cultivars were clustered into 5 clusters with at least 2 per cultivar (Figure 6). The first cluster included the Giza 134, Giza 133, and Giza 136 cultivars. The cultivars Giza 132, Giza 2000, and Giza 128 were discriminated as two neighboring pairs of cultivars. The third cluster consisted of Giza 126 and Giza 137, while Giza 135, Giza 131, Giza 129, and Giza 130 appeared as two neighboring clusters to make up the fourth cluster. The other cultivars—Giza 124, Giza 123, and Giza 125—were located in one group (Figure 6).

Based on the ISSR marker data for the studied cultivars, a genetic distance tree was constructed using Dice's genetic similarity matrix (Figure 7). In this tree, the two pairs (Giza 126 and Giza 2000) and (Giza 132 and Giza 133) were close to the other cultivars. In Egypt, these cultivars are used in human consumption and animal feed. In addition to the malt industry and the beer industry, the ancient Egyptian barley sector dates back to BC. Meanwhile, Giza 132 with Giza 130 and Giza 133 with Giza 130 were less similar to the rest of the barley cultivars, and have been nominated for a crossbreeding program for Egyptian barley breeders. On the other hand, Giza 129 was separated from the rest of the cultivars. All cultivars were distributed in the three clusters. According to the ISSR molecular marker polymorphism, a similarity matrix among the 15 cultivars was derived based on Dice's coefficient (Table 3). According to the similarity matrix of ISSR analysis, the highest similarity value (93%) was observed between (Giza 133 and Giza 132) and (Giza

2000 and Giza 126). Conversely, the lowest similarity value (80%) was recorded between Giza 130 and (Giza 133 and Giza132), indicating that these cultivars were distantly related, as shown in Table 3 and Figure 7. These distinctive cultivars could be expanded to improve soil properties, reduce fertilizer consumption, increase tolerance to drought and salinity, and facilitate growth in newly reclaimed lands. The results were nearly in agreement with those of previous studies [12,42–44].

**Figure 5.** An illustration of the genetic diversity expressed in 15 Egyptian six-rowed barley cultivars, according to a principal component analysis (PCA) based on polymorphism of ISSR markers, using PAST software.

#### *2.4. Biplots*

Biplots were used to reflect the statistical values and their presentation in order to provide supportive information about all of the investigated parameters. Biplots have been used in previous studies to illustrate and present different types of data [45–47]. Through the different types of data, the information can be dispersed, but the biplot distributes the genotypes based on all of the traits under study, whether morphological data or molecular data. The biplot in Figure 8 shows the differences between the clusters in the morphological data and the clusters of the molecular data, as well as their interaction; it also clearly demonstrates the effects of each field trait on the genotypes, along with the effects of each initiator molecule on the Egyptian barley genotypes.

**Figure 6.** Multivariate heatmap illustrating the genetic diversity of 15 Egyptian six-rowed barley cultivars, based on the 15 ISSR primers for using the module of a heatmap of ClustVis—an online tool for clustering and visualizing of multivariate data [41].

To study the interaction between genotype and environment (GE), biplot analysis was utilized [48]. Using the constructed PCA biplot, it became clear which of the 10 morpho-agronomic traits and 15 ISSR primers contributed most to the discrimination of the examined cultivars (Figure 8). The 15 cultivars were divided into 3 groups based on 10 field traits and 15 molecular ISSR primers. The group that included Giza 130, Giza 136, Giza 138, and Giza 126 was the most influenced by the field and morphological characteristics, as shown in Figure 8. This group was established based on maturity day, biological yield, grain yield, weight of 1000 grains, grain field period, and the number of grains. At the same time, the genotypes Giza 129, Giza 137, and Giza 133 were more influenced by the molecular primers associated with age, including ISSRs 807, UBC 835, UBC 826, 851, and UBC 811, as well as hiding day. On the other hand, the third group was affected by plant height characteristics. The number of spikes per square meter, along with the remainder of the molecular parameters, characterized the cultivars Giza 134, Giza 131, Giza 132, Giza 126 Giza 135, Giza 123, Giza 125, and Giza 2000. Generally, when the cultivar falls on the adjective line, it is more impacted by it. Through the current data, we found that the genetic basis of the ISSR molecular markers is dominant over the morphological traits in the first and second groups, while the effect of field traits is predominant in the third and fourth groups, indicating the merging of the field cluster with the molecular cluster into one form in the biplot. Moreover, the contributions of the genes controlling the traits are shown through the molecular parameters, while the environment is shown by the field traits, and the differences in terms of environment and genetics in this study are united by environmental and genetic data [49,50].

**Figure 7.** Cluster tree of genetic distance between 15 Egyptian six-rowed barley cultivars, based on the analysis of 15 ISSR primers according to Euclidean distance and the UPGMA algorithm in PAST software.

476

**Figure 8.** A biplot cluster tree illustrates the genetic distance between 15 Egyptian six-rowed barley cultivars, based on the analysis of 10 morpho-agronomic traits and 15 ISSR primers according to Euclidean distance and the UPGMA algorithm in PAST software.

#### *2.5. DNA Barcoding Loci of mat*K *and rbc*L *Sequencing*

DNA barcoding is an essential tool for species identification [51]. Genes from the chloroplast genome—such as *mat*K and *rbc*L—were used for DNA barcoding. The genetic diversity and phylogeny of the studied cultivars were determined by amplification and sequencing of both loci. Four barley cultivars were used for DNA barcoding. Two cultivars, marked by an asterisk (\*), had a tough, inedible outer hull around the barley kernel (Giza 123 and Giza 138), while two cultivars marked by two asterisks (\*\*) were characterized by sticks and sprouts that separate from the seed when ripe (Giza 126 and Giza 136) (Table 3). There was 100% amplification success with high specificity of PCR amplification of the *mat*K and *rbc*L regions for all four cultivars, as indicated by sharp DNA bands with no byproducts. The recorded size of the PCR product of the *mat*K region was 900 bp, while for the *rbc*L region it was 600 bp (data not shown). The GenBank accession numbers for *rbc*L in Giza 123, Giza 126, Giza 136, and Giza 138 are MW336986, MW391913, MW336987, and MW391914, respectively. The GenBank accession numbers for *mat*K in Giza 123, Giza 126, Giza 136, and Giza 138 are MW336988, MW336991, MW336990, and MW336989, respectively. To confirm the correct amplification of the *mat*K and *rbc*L sequences, a BLAST function was performed, identifying that all of the sequences were strongly coordinated with *mat*K and *rbc*L of the *Hordeum vulgare* sequences. Sun, et al. [52] assessed the possibility of using five intensively suggested regions (*rbc*L, *mat*K, *trn*H-*psb*A, internal transcribed spacer (*ITS*), and *ITS2*) as DNA barcode candidates to differentiate important species of Brassicaceae in China, in order to establish a new digital identification scheme for economic plants of Brassicaceae. They investigated 58 samples from 27 economic species of Brassicaceae for the success of PCR amplification, intra- and interspecific divergence, DNA barcoding gaps, and identification efficiency. Based on their results, the ITS showed superior discriminative ability, with a rate of 67.2% at the species level when compared with other markers.

Pairwise distances were calculated and evaluated based on the conserved *mat*K and *rbc*L gene sequences, using the WebLogo tool [53]. Additional information on the DNA barcoding regions of *mat*K and *rbc*L in four Egyptian six-rowed cultivars of *Hordeum vulgare* is

provided in the Supplementary Materials (Supplementary Figures S1 and S2, respectively); this includes the alignment length, undetermined characters, missing percentages, and variable sites and their proportions, as well as parsimony-informative sites.

Figure 9A illustrates a phylogenetic tree of *mat*K sequence variation using the UPGMA algorithm to discriminate between the four investigated cultivars. To demonstrate the accuracy and efficacy of the created tree, 10 *mat*K sequences were obtained from NCBI and used as outgroups. The tree has three major clusters: the first includes the cultivars in two groups—Giza 123 and Giza 136—and the second includes Giza 126 and Giza 136 (Figure 9B). The third group comprises NCBI outgroup members (Figure 9B). There was a branch length of 1.5, and the bootstrap value displayed next to the branches designates the bootstrap value supporting the node. The Jukes–Cantor method was used to calculate the evolutionary distances based on the base substitutions per site. The bootstrap values were incredibly high (99%), confirming the validity of the tree branching. In the *rbc*L gene region, the four cultivars were distributed into two groups: the first included Giza 123, Giza 136, and Giza 138, while the other contained only Giza 126 (Figure 9C). Additionally, when 10 versions of the identical gene sequences from NCBI were added to GenBank, they resulted in 14 barley genotypes (Figure 9D). The 14 genotypes were distributed into 2 groups: the first group included the cultivar Giza 126 only, while the second group contained the other 13 cultivars (Figure 9D). There was a branch length of 3.5, and the bootstrap value displayed next to the branches designates the bootstrap value supporting the node. The Jukes–Cantor method was used to calculate the evolutionary distances based on the base substitutions per site. The ambiguous plant pairwise deletion option was used. In addition, the rest of the Egyptian cultivars were compared with GenBank's publications, where it was noted that the closest to Giza 123 were the HQ800432 and MN171390 versions. Using DNA barcoding, species could be classified quickly without relying on morphological characteristics. This technique uses DNA fragments of relatively small size as tags to describe or discover species [54].

On the other hand, the MN171392 version was close to Giza 136. Moreover, Giza 138 fell between two versions MN171388 and MN171387 (Figure 9D). After adding NCBI GenBank accession numbers, sequences had 24 genotypes for each region of the *rbc*L and *mat*K genes. The distribution of the GenBank NCBI accession numbers and the Egyptian cultivars did not differ from that of each gene separately from the regions of the *rbc*L and *mat*K genes; however, the similarity percentage was as follows in the *mat*K gene: The GenBank accession numbers of *rbc*L and *mat*K in Giza 123 and Giza 136 were distributed at 99% similarity, whereas the similarity rate of Giza 138 and Giza 126 reached 66% in the area of the genome. Meanwhile, in the *rbc*L gene region, the similarity rate was 56% for Giza 126, while the rest of the GenBank accession numbers and Egyptian cultivars were distributed at a similarity rate of 99%.

The results of the current study show that *rbc*L is less mutable than *mat*K in terms of sequence variability among the examined cultivars. Previous studies used the *mat*K region in many phylogenetic analyses of flowering plants, due to its conservative mode of evolution [55,56]. Four cultivars were differentiated in the present study according to *mat*K sequence variation, using 10 outgroup sequences from NCBI (Figure 9). The phylogenetic tree created using 10 NCBI-extracted *mat*K sequences of *Hordeum vulgare* confirmed the outstanding finding of separating the four cultivars Giza 123, Giza 126, Giza 136, and Giza 138, along with the *Hordeum vulgare* NCBI *mat*K sequence and its subspecies. Nevertheless, Giza 123 and Giza 136 were separated with the *Hordeum vulgare* NCBI *mat*K sequence of the NCBI accession numbers, suggesting sequence homology. In the second cluster, Giza 138 and Giza 126-super-supreme were in the same group, and shared high homology in *mat*K sequences. The *rbc*L region was used to distinguish between wild parents, as well as being used as precise sequences to distinguish between different degrees of biological diversity [57–59]. In addition to providing potentially helpful information for genomeassisted research, the present study also provides useful information for crop improvement.

**Figure 9.** Phylogenetic trees based on (**A**) the *mat*K DNA barcoding region for four six-rowed barley cultivars (Giza 123, Giza 126, Giza 136, and Giza 138); (**B**) the *mat*K DNA barcoding region for four six-rowed barley cultivars, with 10 additional *mat*K sequences of *Hordium vulgaris* L. used as outgroups; (**C**) the *rbc*L DNA barcoding region for four six-rowed barley cultivars; and (**D**) the *rbc*L DNA barcoding region for four six-rowed barley cultivars, with 10 additional *rbc*L sequences of *Hordium vulgaris* L. used as outgroups, using MEGAX software.

#### **3. Materials and Methods**

#### *3.1. Plant Materials*

This study examined 15 Egyptian barley cultivars (all six-rowed). Those cultivars were selected because they are more critical to the Egyptian barley industry than the two-rowed lines. Viable grains of the studied cultivars were obtained from the Barley Research Department (BRD), Field Crop Research Institute (FCRI), Agricultural Research Center (ARC), Giza, Egypt, during two seasons: 2018/2019 and 2019/2020 (Table 4). These cultivars were chosen based on the recommendations of the barley breeders and the beer industry for their salinity and drought tolerance, high yield, and phytochemical characteristics—such as mineral elements and malt content.

#### *3.2. Morphological Traits and Experimental Design*

Two field experiments were carried out at El-Giza Agricultural Research Station (Giza, Egypt) during the successive winter seasons of 2018/2019 and 2019/2020 to study the morphological traits of the different cultivars. To differentiate between the studied cultivars based on morphological characteristics, the following parameters were recorded: days to 50% heading (HD), days to 50% maturity (MD), grain filling period (GFP) (days), plant height (PH) (cm), spike length (SL) (cm), number of grains per spike (average of 10 spikes per square meter), number of spikes per m2 (No. Sp./m2), weight of 1000 grains (g), biological yield (BY) (t/ha), and grain yield (GY) (Kg/ha). The grain filling period (GFP) was calculated using the following formula:

#### *Grain filing period* = *maturity days* − *f lowering day*

A randomized complete block design (RCBD) with four replications was used. The plot size was 4 rows that were each 3 m long and 20 cm apart. Analysis of variance and least significant difference (LSD) at 5% were used for comparison between the cultivars.

**Table 4.** Name, origin, and year of release of the Egyptian six-rowed barley cultivars as recorded by the Barley Research Department, Field Crops Research Institute, Agricultural Research Center, Egypt, and the GenBank accession numbers for the *mat*K and *rbc*L genes for four six-rowed barley cultivars (Giza 123, Giza 126, Giza 136,andGiza138).

NA: not available; \*: has a tough, inedible outer hull around the barley kernel; \*\*: hull-less—the sticks and sprouts separate from the bean when ripe, and the chromosome number is 2n = 2×=14(mapview:"barleygenomeatncbi.nlm.nih.gov";retrieved6October2014).

#### *3.3. Scanning Electron Microscopy (SEM)*

Four barley cultivars (Giza 123, Giza 126, Giza 136, and Giza 138) were studied using SEM. Those four cultivars were chosen based on the recommendations of the plant breeders. The chosen cultivars have high production demand and can withstand harsh conditions; they also have excellent synthetic qualities, which is the reason for their examination. For example, Giza 123 tolerates harsh conditions and high salinity levels; Giza 126 has excellent drought tolerance, and is grown under the rain on the northern coast of Egypt, while Giza 136 and Giza 138 are characterized by high yield under all conditions. Viable grains of the studied cultivars were obtained during the season of 2019. The clean and dry seed samples of the studied barley cultivars were placed on double-stick tape mounted on a copper electron microscope holder. The specimens were coated with gold, and then investigated and photographed with a JEOL JSM T200 at 25 kV, in the electron microscope unit of Mansoura University, Mansoura, Egypt. Seed coat technical terms were based on the works of Koul, et al. [71],Murley [72],and Stearn [73].

#### *3.4. ISSR Molecular Markers*

#### 3.4.1. Extraction of Genomic DNA

Fresh leaf tissue (0.1 g of combined samples from three different plants) ground in liquid nitrogen with a mortar and pestle was used to extract genomic DNA using the cetyl trimethyl ammonium bromide (CTAB) protocol [74]. DNA concentration and purity for all samples were determined spectroscopically at 260 and 280 nm, respectively. DNA samples were stored at −20 ◦C for subsequent molecular analysis.

#### 3.4.2. ISSR Amplification

ISSR amplification reactions were carried out in equal volumes (15 μL) containing 7.5 <sup>μ</sup>L of 2× Master Mix (OnePCRTM, GeneDireX, Inc., Taipei, Taiwan), 1 <sup>μ</sup>L of DNA template (10 ng/μL), and 1 μL of primer. The names and sequences of the ISSR primers used in the current study are listed in Table 2. The amplification reaction was performed using a T100TM Thermal Cycler (Bio-Rad® Laboratories, Hercules, CA, USA). The polymerase chain reaction (PCR) program was as follows: initial denaturation at 94 ◦C for 4 min, followed by 30 cycles, with the first step at 94 ◦C for 30 s (denaturation), the second step varying between 46 and 52 ◦C—depending on the GC content of each primer—for 45 s (annealing), and the third step (extension) at 72 ◦C for 1 min, followed by a final extension step at 72 ◦C for 7 min. The reaction was stopped by maintaining the tubes at 4 ◦C for at least 30 min. Amplification products were separated via electrophoresis on 1.5% agarose gel in 1× TBE buffer (Tris-borate-EDTA). The gels were stained with 0.5 μg mL−<sup>1</sup> ethidium bromide (EtBr) solution (Thermo Fisher Scientific, Carlsbad, CA, USA). Then, the gel was documented using a Bio-Rad ChemiDocTM MP gel documentation system (Bio-Rad). The primers that gave reproducible results were used for data analysis. Polymorphism indices were calculated using iMEC (Online Marker Efficiency Calculator) (https://irscope.shinyapps.io/iMEC/) [75]. ClustVis, a web tool for visualizing clustering of multivariate data, was used to construct heatmaps (https://biit.cs.ut.ee/clustvis/) [41].

#### *3.5. DNA Barcoding of Plastid Genes rbc*L *and mat*K

DNA barcoding of sequences for the *rbc*L and *mat*K genes was performed using computational analysis. BioEdit software version 7.2.5 (https://bioedit.software.informer. com) was used to analyze and assemble the *rbc*L and *mat*K gene sequences for every cultivar. Using the BLAST function (https://www.ncbi.nlm.nih.gov), the sequences were compared with all accessible sequences in the database. The primers used for barcoding of the *rbc*L and *mat*K genes are listed in Table 5. The PCR program to amplify the two genes was as follows: initial denaturation at 94 ◦C for 4 min, followed by 40 cycles, with a denaturation step at 94 ◦C for 30 s, annealing step at 45 ◦C for 30 s, and elongation step at 72 ◦C for 30 s, followed by a final extension step at 72 ◦C for 7 min, after which it was maintained at 4 ◦C to stop the reaction. The PCR products were subsequently electrophoresed on 1.5%

*w/v* agarose, stained with 0.5 <sup>μ</sup>g mL−<sup>1</sup> EtBr solution (Thermo Fisher Scientific) in 1× TBE buffer, and visualized as described for the ISSR PCR amplification. The PCR products of the *mat*K and *rbc*L genes were recovered from agarose gel and purified using the Monarch DNA Gel Extraction Kit (New England Biolabs, Inc., Ipswich, MA, USA), according to the manufacturer's instructions. The purified *mat*K and *rbc*L amplicons were cloned into pGEM®-T Easy Vector Systems (Promega Corporation, Madison, WI, USA) before sequencing. After being transformed into the competent cells of the *E. coli* strain DH5α (Promega, Madison, WI, USA), the positive recombinants were identified via anti-ampicillin selection and verified by PCR screening. Three of the positive clones were sequenced using the ABI PRISM Big Dye Terminator Cycle Sequencing Ready Reaction Kit (Applied Biosystems, Waltham, MA, USA) in conjunction with ABI PRISM (3100 Genetic Analyzer, Macrogen DNA Sequencing Services, Seoul, Korea), as described by Badr, et al. [76]. Using Gblocks software version 0.91b, the revealed nucleotide sequence was assembled [77,78].

**Table 5.** Primer names, sequences, and product sizes for the *rbc*L and *mat*K genes' DNA barcoding.

Online ClustalW2 software (https://www.ebi.ac.uk/Tools/msa/clustalw2/) was used to align multiple nucleotide sequences, which were double-checked using MEGAX (www.megasoftware.net). Gblocks version 0.91b [77,78] was used to review and assess the gaps in the positions. MEGAX software using the UPGMA algorithm was used to perform the phylogenetic analysis. Confidence of the clustering was attained using SEQBOOT (https://csbf.stanford.edu/phylip/seqboot.html). The sequence logos of the multiple sequence alignments were generated using the WebLogo tool [53]. Additionally, a principal component analysis (PCA) biplot based on the morpho-agronomic data matrix was constructed via multivariate analysis using PAST software versiong 4.02 (https:// www.nhm.uio.no/english/research/infrastructure/past/).

#### *3.6. Statistical Analysis*

Standard analysis of variance (ANOVA) using least significant differences (LSD) was utilized to estimate the significant differences between the 15 cultivars of six-rowed barley [81]. Dendrogram cluster analysis was used to arrange a set of variables into clusters. A cluster analysis was performed using Euclidean distance and similarity levels [82,83]. ISSR markers that generated clear, distinct, and reproducible bands were recorded as (0) for absence or (1) for presence. The ability of ISSR primers to differentiate between investigated genotypes was analyzed by calculating the polymorphic information content (PIC) [84]. Resolving power (Rp) was measured following the formula of Gilbert, et al. [85]. Additionally, marker index (MI) and effective multiplex ratio (EMR) values were calculated. For the calculation of the coefficient of genetic similarity matrix, and for the construction of a distance tree illustrating the relationships between the tested genotypes, the ISSR marker matrices were used in combination with the unweighted pair group method with arithmetic mean (UPGMA) in PAST software version 4.02 [86]. Furthermore, by using PAST software version 4.02 [86], a PCA scatter diagram was constructed based on a Dice coefficient genetic similarity matrix. ClustVis, a web tool for visualizing clustering of multivariate data, was used to constructe heatmaps (https://biit.cs.ut.ee/clustvis/) [41].

#### **4. Conclusions**

Barley plays a vital role in Egypt in terms of agricultural development and added value, as it is used in new and marginal lands. Today, there is expansion in its cultivation due to its high adaptation to water scarcity and other harsh conditions. Thus, barley is added to wheat flour to increase the nutritional value of bread, and is also used in the manufacture of beer and malt. These cultivars are the most common forms of barley in Egypt, so they have been studied for their economic importance in terms of added value and sustainable development. Despite the differences at the molecular level, the examined Egyptian cultivars reflected similarities in terms of field performance under the optimal environment, exhibiting no differences in terms of field characteristics. These cultivars were closely distributed in a genetic tree, similar to the genetic tree based on the molecular description. These differences enable the breeders to choose the best of these cultivars from the most divergent, and to exclude the least different. Moreover, the electron microscope examination reflected differences in the seed surface characteristic, which helps in understanding the chemical content of the Egyptian barley grains and their economic importance. Interestingly, the sequencing results of four cultivars showed that the *rbc*L gene referred to the uniqueness of these four cultivars compared to the sequence database. Nevertheless, the second gene *mat*K revealed that these cultivars are very similar to the GenBank accession numbers. Additionally, the production of new sequences was added to the molecular information about the Egyptian barley cultivars, showing the differences between the Egyptian and European cultivars—especially since Egypt is one of barley's countries of origin. These results will potentially enhance breeding programs and aid in the development of new adaptive or high-yield barley cultivars with specific improved traits.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/ 10.3390/plants10112527/s1: Figure S1: Sequence logos of the multiple sequence alignment of *mat*K; Figure S2: Sequence logos of the multiple sequence alignment of *rbc*L.

**Author Contributions:** Conceptualization, A.H.M., A.A.O. and E.M.Z.; methodology, all authors; software, A.H.M., A.A.O., A.M.A., R.M.R. and E.M.Z.; validation, A.H.M., A.A.O., M.M.A.E. and E.M.Z.; formal analysis, A.H.M., A.A.O., R.M.R. and E.M.Z.; investigation, A.H.M., A.A.O. and E.M.Z.; resources, A.M.A., M.M.A.E. and E.M.Z.; data curation, A.M.A., R.M.R. and M.M.A.E.; writing—original draft preparation, all authors; writing—review and editing, A.H.M., A.A.O. and E.M.Z.; visualization, A.A.O. and E.M.Z.; supervision, A.H.M., A.A.O. and E.M.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** All the obtained data presented in this the article or supplementary material.

**Acknowledgments:** The authors would like to thank Abdoallah Sharaf, Iva Mozgová's Lab, Institute of Plant Molecular Biology, Biology Center ASCR, Budweis, South Bohemia, Czechia for helping with the WebLogo analysis of the *mat*K and *rbc*L genes.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **RWLMod—Potential Model to Study Plant Tolerance in Drought Stress Conditions**

**Florin Sala 1, Mihai Valentin Herbei 2,\* and Ciprian Rujescu <sup>3</sup>**

**Abstract:** Rationale: Water loss by evaporation is a normal physiological process, in order to regulate plant temperature. Under conditions of thermal and water stress, water loss is accelerated compared to normal conditions, and the response of plants is variable. In extreme cases, it can lead to wilting and death of plants. It was found that the phenomenon of water loss behaved as a pattern in different plant species, given by two functions, logistics (first part of water loss) and hyperbola (second part of water loss) in relation to a moment m, at which the rate of water loss (RWL) has reached its maximum value. Method: We studied the water loss process for a series of plant samples on different plant species (*Picea abies* L., H. Karst; *Juniperus communis* L.; *Pinus silvestris* L.; *Thuja occidentalis* L.; *Lamium purpureum* L.; *Veronica hederifolia* L.), measuring the rate of weight loss (RWL) in controlled conditions. The drying of the samples was done in identical conditions (thermo-balance, 100 ◦C, standard temperature for drying the plant samples) with the real-time recording of the drying time simultaneously with the water loss rate (RWL) from the plant samples. The exposure time varied, depending on each species sample, and was approximately 1000 s. Results: The experimental data was recorded at intervals of every 10 s, during the entire drying period. RWL values varied from 0.024 to 0.054 g/min at the beginning of the drying process and reached maximum values after 70–100 s, having values between 0.258 g/min and 0.498 g/min. During the drying period, this indicator presented different graphic evolutions, difficult to be described with a single function. The first segment was described by a logistic function, and the second was described by a hyperbola, resulting in a model (RWLMod) which described the real phenomenon. This model and theoretical calculation were used to quantify the water loss in a time interval and, compared with empirical dates, no significant differences were observed, which indicated an increased degree of accuracy regarding the use of this model. Recommendation and novelty of work: The novelty of the work is given by the obtained model (RWLMod), which makes possible the description of RWL over the entire time interval, and ensures a good fit with the real data. It recommends the method and model in studies of plant behaviour under stress in relation to different influencing factors.

**Keywords:** drought stress; drying processes; mathematical model; plant hydric stress tolerance; rate of weight loss; RWLMod; water evaporation

#### **1. Introduction**

Water has a vital role in plant life, in relation to physiological and metabolic processes, plant nutrition, thermoregulation, plant growth and development, tolerance to stressors, etc. Knowing the dynamics of water loss in plants has multiple applications,

**Citation:** Sala, F.; Herbei, M.V.; Rujescu, C. RWLMod—Potential Model to Study Plant Tolerance in Drought Stress Conditions. *Plants* **2021**, *10*, 2576. https://doi.org/ 10.3390/plants10122576

Academic Editor: Milan S. Stankovic

Received: 24 September 2021 Accepted: 23 November 2021 Published: 25 November 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

beginning with a better understanding of the crop behaviour under conditions of water stress; we also identified technological issues on improving methods for the processing of plant products. The literature contains studies on the physiological mechanisms of water loss [1,2] on increasing tolerance regarding hydric stress [3–5], and on the drying process of some aromatic plants with economic possibilities [6]. In addition, the dynamics of plant humidity was studied under the influence of some external factors, such as air speed [7], or physiological and biochemical factors [8,9]. In the case of crop plants, there was a special interest regarding the increase of plant resistance to water stress and thermal factors; a number of studies have analysed physiological indices [10–12], water efficiency in plants [13,14], photosynthetic capacity and production quality [15,16], regarding the increased demands for food production, in the context of population growth and climate change [17,18]. There is a confirmed existence for the particular well-defined dynamic of some indicators describing humidity, and existing studies are using a mathematical characterization of these processes. Thus, the variation in humidity is behind the creation of some mathematical models, using the study of isothermal curves [19,20], or described using sigmoid curves [21]. In fact, sigmoidal curves that form a part of the mathematical model used in this paper are found in certain distinct classifications of nonlinear models used in agricultural sciences [22]. Phenomena with a downward asymptotic trend have been evaluated in studies regarding the behaviour of biological and biochemical processes [22]. Moreover, the necessity of using functions based on mathematical models which describe specific situations as accurately as possible require the adjustment and modification of classical models. Regarding old exponential models or those with limited growth, there are multiple concerns for setting functions to provide such an attribute. Power-Ricker and modified logistic function [23,24] are such examples, with applications in various fields of biology.

The problem of plant water loss was addressed in this study, in terms of rate of weight loss (RWL), under controlled conditions. By measuring the rate of water loss (RWL) in different plant species, the behaviour of the process was found according to a pattern, regardless of the species studied. A logistic function describes the first part of the process (RWL), and a hyperbolic function describes the second part of the process (RWL), in relation to the RWLmax value, recorded at all samples, but at different times. The typical approach in this study led to the finding of a mathematical model that described RWL throughout the drying process. This model (RWLMod) is the solution of this study in the description of RWL in plants. The model would facilitate the study and better understanding of the behaviour of plants (including crop plants) in conditions of water and heat stress, such approaches all the more necessary in the context of climate change.

#### **2. Results**

For each plant species studied, the series of data recorded in real time (Tables 1–6) on the rate of weight loss (RWL) expressed in grams/minute, respectively weight (g) and drying time, are presented. These values have been obtained for the corresponding time instants of the measurements, carried out at intervals of ten seconds. Similarly, for each sample, the maximum value of RWL has also been distinctly indicated, represented as RWLmax. The total duration of the process was specific to each sample, observing at the end of the tables the moment in time when the experiment stopped automatically, when the amount of water lost had become negligible. However, the upper limit of the time was about 1000 s. These data series were the basis for the statistical determination of the coefficients of the functional models.

**Table 1.** Statistical data series on water loss for the sample *Picea abies* L., H. Karst (Spruce).

\* registered value of RWLmax.

**Table 2.** Statistical data series on water loss for the sample *Pinus silvestris* L. (Pine).

\* registered value of RWLmax.

**Table 3.** Statistical data series on water loss for the sample *Juniperus communis* L. (Juniper).

\* registered value of RWLmax.

**Table 4.** Statistical data series on water loss for the sample *Thuja occidentalis* L. (Thuja).

\* registered value of RWLmax.

**Table 5.** Statistical data series on water loss for the sample *Lamium purpureum* L. (Nettle).

\* registered value of RWLmax.

**Table 6.** Statistical data series on water loss for the sample *Veronica hederifolia* L. (Veronica).

\* registered value of RWLmax.

Rate of weight loss has values between 0.024–0.054 g/min at the initial moment in time (t = 10 s), and the highest value was represented in the nettle. Afterwards, an increased growth rate was observed for this indicator, thus as approximately 70–100 s from the debut; RWL had values between 0.258 g/min for pine and 0.498 g/min for nettle. The limitation stage was observed at the end of this time interval and for a short period (10–20 s). The following stage regarded RWL variation; a rapid decrease on the unit of time, from the

maximum values mentioned above until the values corresponding to the moment of time t of approximately 150 s, after which RWL had a slow trend of decrease until the end of the exposure period.

The following tables, respectively one for each sample (e.g., statistically calculated coefficients, corresponding to the functional model, sample "n" in Tables 7–12, show for each functional model the values of the coefficients, statistically calculated, then the values rsq. (R<sup>2</sup> coefficient of determination), and. sig. (significance probability) for testing the accuracy from the statistical point of view. Sample "n" indicates each plant species studied. For the first branch of the function, all values are superior to 0.990, indicating an almost perfect fitting, and all sig. values are smaller than 0.001, indicating a high degree of accuracy. In addition, for the second branch, the values rsq. are high, rsq. = 0.804 being, in fact, the lowest value. Further, for each individual sample, sig. < 0.001.

**Table 7.** Statistical coefficients calculated, corresponding to the functional model, sample "spruce".

Note. Sig.—significance probability; Rsq.—R2 coefficient of determination; a, b, α, β—coefficients of the function (8); u—upper bound; t1 and t2—time interval; I—values returned by equation (12) for the mentioned time moment; J—values returned by Equation (13) for the mentioned time moment.

**Table 9.** Statistical coefficients calculated, corresponding to the functional model, sample "juniper".

**Table 10.** Statistical coefficients calculated, corresponding to the functional model, sample "thuja".

**Table 11.** Statistical coefficients calculated, corresponding to the functional model, "nettle".

**Table 12.** Statistical coefficients calculated, corresponding to the functional model, sample "veronica".

Moreover, we presented both values of water losses; on the one hand, those resulting from the theoretical integral calculations, and on the other hand, the real values obtained by real-time determinations. For each sample, separately, the results were presented distinctly, for each branch of the function. The data were statistically tested on the differences between the groups determined by the theoretical and empirical method, using the Mann– Whitney test in SPSS. The values obtained are U = 14.5, sig. = 0.575 for the time segment corresponding to the first branch (logistic model), respectively, and U = 15, sig. = 0.631 for the time segment corresponding to the second branch (hyperbole), indicating the acceptance of the null hypothesis. Therefore, it can be considered that the two data sets do not differ.

At the end of the presentation for each sample, the expression of the functional model is represented by a graph (rate of weight loss (g/min), sample, in the Figures 1–6. Here, the value of the theoretical maximum point (t = m) can also be observed, thus indicating the time (theoretical) when the maxim (theoretical) of RWL takes place. Thus, for the

sample "spruce", the nonlinear Equation (1) was solved, resulting in the value: t = m = 81.2 s. Thus, for the sample "spruce", we have the Equation (2).

$$\frac{1}{\frac{1}{0.99} + 99.5066 \cdot 0.9127^t} = -0.0031 + \frac{31.1986}{t} \tag{1}$$

$$\mathbf{f}(\mathbf{t}) = \begin{cases} \frac{1}{\frac{1}{0.99} + 99.5066 \cdot 0.9127 \,\mathrm{t}}, & \mathbf{t} \le 81.2\\ -0.0031 + \frac{31.1986}{\mathrm{t}}, & \mathbf{t} > 81.2 \end{cases} \tag{2}$$

Similar calculations were made solving nonlinear equations corresponding to the other samples and resulting in the functions below. For pine, Equation (3) was found, with m = 80.6 s.

$$\mathbf{f}(\mathbf{t}) = \begin{cases} \frac{1}{\frac{1}{0.55} + 112.688 \cdot 0.9141^{\prime}}, & \mathbf{t} \le 80.6\\ 0.0094 + \frac{26.658}{\mathbf{t}}, & \mathbf{t} > 80.6. \end{cases} \tag{3}$$

For the juniper sample, the resulted model is the Equation (4), and m = 92.5 s.

$$\mathbf{f}(\mathbf{t}) = \begin{cases} \frac{1}{\frac{1}{0.43} + 74.3438 \cdot 0.9239}, & \mathbf{t} \le 92.5\\ -0.0129 + \frac{41.5086}{\mathbf{t}}, & \mathbf{t} > 92.5 \end{cases} \tag{4}$$

For thuja, Equation (5) was found, m = 74.8 s.

$$\mathbf{f}(\mathbf{t}) = \begin{cases} \frac{1}{\frac{1}{0.41} + 109.28 \cdot 0.8973^{\circ}}, & \mathbf{t} \le 74.8\\ -0.005 + 32.2814, & \mathbf{t} > 74.8 \end{cases} \tag{5}$$

Then, for nettle, the model is Equation (6), m = 86.6 s.

$$\mathbf{f(t)} = \begin{cases} \frac{1}{\frac{1}{0.57} + 27.6696 \cdot 0.9262^{\prime}}, & \mathbf{t} \le 86.6\\ -0.0136 + \frac{44.5920}{\mathbf{t}}, & \mathbf{t} > 86.6 \end{cases} \tag{6}$$

Respectively, Equation (7) was found for Veronica, and m = 109.8 s.

$$\mathbf{f}(\mathbf{t}) = \begin{cases} \frac{1}{\frac{1}{477} + 51.086 \cdot 0.9342^{\prime}}, & \mathbf{t} \le 109.8\\ -0.0147 + \frac{49.2996}{\mathbf{t}}, & \mathbf{t} > 109.8 \end{cases} \tag{7}$$

**Figure 1.** Rate of weight loss (g/min) in spruce samples.

**Figure 2.** Rate of weight loss (g/min) in pine samples.

**Figure 3.** Rate of weight loss (g/min) in juniper samples.

**Figure 4.** Rate of weight loss (g/min) in thuja samples.

**Figure 5.** Rate of weight loss (g/min) in nettle samples.

**Figure 6.** Rate of weight loss (g/min) in veronica samples.

#### **3. Discussion**

Multiple concerns have existed to study and describe aspects related to water in plants, or plant products, and the modeling approach has been the basis of many methods and techniques of investigation. The drying process is a method commonly used for conditioning plants. The generally accepted definition in the literature is reducing the moisture content of a certain product. The widespread use appears in preparation techniques for the medicinal herbs in beverages; in the beginning, fresh biomass (herba) will be subject to drying processes, a mandatory stage which reflects directly in the quality of the finished product, the period during which the product can be stored without quality depreciation, and also multiple economic aspects of generating profits [25].

Müller and Heindl [25] have studied the drying parameters for *Salvia officinalis.* One of the conclusions reported was related to the water activity (aw). It correlates with the relative humidity (RH) of the air in the areas adjacent to the material to be analysed, therefore over the limit of RH > 70%, the development of some bacteria, fungus, residue (lees) has been noticed—issues that have a direct effect on the state of product quality.

The transformations occurring at a fixed temperature are also a topic often analyzed. Thus, sorption isotherms are described mathematically for *Artemisia dracunculus* with the help of the Halsey equation [25], but also with the help of other known models: BET, Caurie, GAB, Halsey, Henderson, Lewicki, Modified Mizrahi, Oswin, Peleg [26]. In addition, extensive isotherms for a series of plants were found, such as *Salvia officiinalis* [25], *Artemisia dracunculus*, *Mentha piperita*, *Thymus vulgaris* [26–28], *Mentha crispa* [29]), *Mentha viridis*, *Salvia officinalis*, *Lippia citriodora* [30], *Ficus deltoidea* [31]), *Melisa officinalis* [32]), *Chenopodium ambrosioides* [33]), *Citrus sinensis* [34]), *Ziziphus spina-christi* [21]), *Phylanthus ambelica*, and *Zingiber officinale* [20]).

Furthermore, fruits are often subjected to drying processes. Simal et al. [35] and Kaya et al. [36] describe, using simulation methods, the dynamic water loss from kiwifruit; similar studies regarding various tropical fruits were conducted by Ceylan et al. [37] or Fernando and Amarasinghe [7].

This subject seems to be in full evolution at the moment, especially considering that modern food technology uses a high range of plants, for which the conducted studies are still not sufficiently detailed regarding water loss.

Rate of weight loss (RWL), defined as the amount of water lost in the time interval, can be discussed through physical analogy with the speed–vector size, thereby creating the possibility for it to be described with the help of mathematical physics equations. Jones and Sleeman [38] studied various biological models described by equations of this type.

Moreover, the rate of water loss is not constant; it is different from one species to another, from one organ of the plant to another, the way water content is present in the certain product (free water at intracellular level or present as links with various compounds).

The drying temperature, humidity and air velocity directly influence this parameter, and depending on a specific practical purpose, there are multiple studies that indicate the recommended parameters for drying in order to improve or protect some useful active principles.

Rate of weight loss of water in plants is less studied in a direct manner, but more indirectly by deduction from other calculations. However, here it has presented great interest mainly because of the different mathematical approaches that can be applied to this indicator. For example, determining the total quantity of water lost could be determined by using the properties of the definite integral.

At the same time, studying the RWL, there is the possibility of making direct observations in order to rapidly reduce the content of free water in plants. This fact can be used specifically in order to reduce the risk of developing microorganisms with undesirable effect over the quality of plant products.

In similar research to the one proposed in this study, Kaya and Aydin [39] investigated mint (*Mentha spicata* L.) and nettle (*Urtica dioica*) leaves in order to describe the loss of humidity in the presence of some variable external factors: air temperature, air velocity and relative humidity. Temperatures were lowered closer to the natural environment, with working values being 35, 45, 55 ◦C, and a longer drying period (30–50 min). Study regarding the leaves of lemon grass (*Cymbopogon citratus Stapf*) is also worth mentioning, although this was focused on the decrease of humidity, and various functional models that describe the trend shown by the experimental data were presented [40].