Development of SSR Markers and Construction of a Genetic Linkage Map in Lavender (Lavandula angustifolia Mill.)

1. Introduction

The lavender (Lavandula angustifolia Mill.), lavandin (Lavandula x intermedia Emeric ex Loisel.) and to some extent spike lavender (Lavandula latifolia Medik.) are widely industrially cultivated members of the genus Lavandula worldwide [1,2,3,4,5]. The industrially cultivated lavender and lavandin are mainly used for production of essential oil and to a lesser extent of concrete and hydrolate [2,3,6]. Due to its unique fine fragrance the lavender oil is widely used in perfumery, cosmetics, aromatherapy and medicinal applications [7,8,9,10]. Cultivated lavandin provides higher yield of essential oil per acre than lavender, but its sensory profile is less sophisticated and attractive which assigns the use of lavandin oil mainly in the industrial perfumery. Bulgaria, France and China have been the largest lavender oil suppliers in the world during the last decade [3,11], whereas the production of lavandin oil has been dominated by France [2,4]. The requirement, that the produced essential oil must meet the ISO standards for lavender (ISO 3515:2002) and lavandin (ISO 3054:2017) oil, has resulted in the growing demand for production and use of high-quality planting material from superior vegetatively propagated varieties [2,11]. The increased use and demands of lavender and lavandin oil further attracts farmers and companies from various countries to expand the industrial cultivation and processing of lavender and lavandin, including the start of cultivation in new areas [6,11,12,13,14]. The expansion of lavender and lavandin farming areas is often related to their cultivation in soil and/or climatic conditions that significantly differ from the optimal for the varieties initially developed for the environmental conditions specific for the traditional areas of cultivation [13,15]. Often such cultivation at sub-optimal conditions of the available varieties results in large variations of flower and essential oil yields and quality, and significant economic losses. In addition, the introduction of new technologies for mechanization of lavender and lavandin cultivation, processing and extraction often requires the use of varieties suited for the application of the particular technology [6,16,17,18]. The described above suggests the need for development of genomic resources and related molecular marker tools for better characterization of the available lavender and lavandin genetic resources, identification of QTLs for important agronomic and essential oil quality traits and application of marker assisted selection for efficient breeding of new elite varieties meeting the new challenges of the industry [2,11,19,20].

The biosynthesis and accumulation of volatile organic compounds in lavender and lavandin flowers were the subject of a number of studies which resulted in the identification of a set of genes involved in particular biosynthetic steps, for example [21,22,23,24,25,26], or the regulation of these processes [27,28]. These studies were further supported by transcriptome sequencing [29,30,31] and genome-wide characterization of gene families related to the terpenoid biosynthesis pathways [32,33,34], as well as by the reported nuclear genome sequence assemblies of three L. angustifolia varieties including var. Maillette [35], var. Jingxun 2 [36] and the reference genome sequence of var. Munstead [37] accessible through the GenBank database. Accordingly, lavender was proposed as a model of aromatic and medicinal plants to study the biosynthesis of volatiles and the related ecology and evolution [15,31,38,39,40]. So far, characterization of the genetic diversity and resources of Lavandula species has been the subject of few studies employing EST-SSR [41], SRAP [42], and RAPD [43] markers.

Over the last three decades construction of genetic linkage maps and QTL analysis has been widely applied as an effective tool for the identification of loci for important agronomic traits and direct application of the obtained results in marker-assisted selection programs in various crops including medicinal and aromatic plants [44,45]. Although the construction of genetic maps increasingly involves the application of SNP based markers, the genetic linkage maps based on SSR markers offer the attractive opportunity to be used as core genetic maps for different genotypes of the species of interest, due to the high transferability of SSR markers. In spite the large number of molecular genetics and genomic studies and related data in lavender and lavandin, surprisingly no construction of a genetic linkage map has been reported for them. Moreover, the only SSR markers reported so far is a set of EST-SSR markers, which are not sufficient for genome mapping studies and QTL analysis.

In the present study, we report the development of an extended set of SSR markers for lavender and the construction of the first genetic linkage map for this species following genotyping of an F2 segregating population obtained by self-pollination of the industrial variety Hemus.

2. Results and Discussion

2.1. Development of an F2 Segregating Population by Self-Pollination of L. angustifolia var. Hemus

Variety Hemus of L. angustifolia has been established in 1974 and during the years has become one of the most widely industrially cultivated lavender varieties in Bulgaria. It is suitable for mechanized harvesting, tolerates cultivation in different geographic regions in the country [46] and provides high-quality of the distilled lavender oil. For some period of time, var. Hemus was also used as a state standard for evaluation of the performance of newly developed lavender varieties. Since its initial development var. Hemus has been propagated vegetatively, both for maintenance in genetic resources collections and for the production of planting material. Considering all this, var. Hemus was selected for the development of a F2 segregating population and genetic linkage map of lavender. Our early attempt to develop F2 segregating population through spontaneous or mechanical self-pollination of a single var. Hemus plant was not successful. The spontaneous self-pollination set up by covering the entire plant with a cage of fine mesh for isolation of insects resulted in the formation of no seeds, supporting the recently reported lack of spontaneous self-pollination in lavender [47]. The mechanical self-pollination of flowers placed under an isolator produced only a small number of seeds, not sufficient for the development of a large segregating population. Therefore, to generate a larger seed set after self-pollination, a single var. Hemus plants was grown under open field conditions, with no other L. angustifolia detected in the same area. As a result, over 1200 seeds were produced from a single three years old var. Hemus plant during the 2020 flowering season, due to insect-mediated self-pollination. The germination and planting of these seeds resulted in the generation of a segregating population of a thousand plants. The SSR marker genotyping of 110 randomly chosen plants from this population and comparison of the scored alleles with those of the starting var. Hemus plant showed that only one of these plants had alleles not observed in the mother plant. This suggests that the plant most likely originated from cross-pollination of the mother plant with unknown L. angustifolia plant mediated by long-distance pollen transfer by the insects. Accordingly, 95 of these tested plants were selected as a mapping population. The results demonstrated that such open-field self-pollination of a single plant could be successfully employed for the generation of a segregating population. Moreover, the open-field cultivation could be substituted by growing plants in insect restrained room and the use of dedicated greenhouse insects as pollinators. Such an approach could be further applied for the generation of self-pollinated seeds for plants of other species from the Lamiaceae family, already reported as difficult for self-pollination [48,49].

2.2. SSR Marker Development

Three groups of primer pairs were tested for amplification of genomic regions containing SSR sequences and development of a large set of SSR markers for genotyping of the segregating population, including: 471 NGS-SSR primer pairs determined after analysis of Illumina MiSeq generated genomic sequences of var. Hemus, 170 GEN-SSR primer pairs determined after analysis of the reference genome sequence of var. Munstead retrieved from GenBank and a set of EST-SSR primer pairs reported by Adal et al.l [41]. The testing included two rounds, first for successful amplification of the target SSR region from DNA isolated from var. Hemus and production of a distinct profile of the amplified DNA fragments, and second round for generation of segregating alleles after analysis of 8 plants from the segregating population. The results of the tests are demonstrated in Table 1 and show that testing of the NGS-SSR group was the most efficient for development of polymorphic SSR markers. The testing of markers from the EST-SSR group resulted in the identification of a significantly lower portion of polymorphic markers, which support the higher efficiency and informativeness of the genomic SSRs compared to EST-derived SSR markers marker observed also in other studies [50,51]. The testing of SRAP markers showed that this type of markers could be readily applied for simultaneously genotyping of a higher number of loci, due to the higher number of alleles which could be detected and analyzed by a single marker. Thus, SRAP markers could be efficiently applied to increase the marker density of SSR-developed genetic maps.

2.3. Linkage Map Construction

A total of 410 loci were analyzed for linkage mapping, including 377 loci obtained through genotyping with three groups of SSR markers and 33 loci identified through analysis with SRAP markers. The first round of map construction involved genotyping of the segregating population with 5 previously reported EST-SSR markers [41] and 255 NGS-SSR markers developed after analysis and testing of SSR regions in genomic sequences of L. angustifolia var. Hemus. This first version of the map was developed shortly before the publication in GenBank of the reference genome sequence of L. angustifolia var. Munstead. Using the latter, BLAST search with the primer and SSR region sequences in the reference genome sequence was used to affiliate the constructed LG maps to the corresponding reference chromosome sequences and for determination the physical position of the mapped SSR loci. The obtained data were used to proceed in a second round of map expansion through identification, selection and testing GEN-SSR primer pairs for amplification of SSR regions located in the chromosome regions outside of the range of the corresponding linkage groups or located in the low marker density regions of the map. All this resulted in mapping of additional 90 loci after analysis with 79 GEN-SSR markers. Finally, the segregating population was genotyped with 11 SRAP markers which resulted in the mapping of 33 SRAP loci. The overall data analysis and quality check of the obtained genetic map resulted in the identification of a total of 25 LGs, which correspond to the number of assembled chromosome sequences in the reference genome, as well as to the chromosome number of var. Hemus used for the development of the segregating population, Figure 1. The main parameters of the map are shown in Table 2 and details for the position of each locus are presented in Table S3, Supplementary Materials. The constructed genetic linkage map of L. angustifolia var. Hemus contains 375 loci covering a total of 2518,9 cM. The average marker distance in the map is 7,05 cM. The number of markers on each LG ranged from 9 (LG 13 and LG 23) to 26 (LG 3), with an average number of 15, and the length per LG ranged from 54,8 cM (LG 25) to 172,2 cM (LG 6), with an average length of 100,8 cM. A total of 35 analyzed loci were not included in the map. Of them, 16 loci had no distinct affiliation to the established LGs, Table S3. The other 19 loci initially affiliated to different LGs were excluded from the final map since they showed segregation distortion and their inclusion resulted in significant changes of the marker order in the maps of the corresponding LG, Table S3. The appearance of loci with segregation distortion which significantly change the established map after their inclusion was reported for different plant species, for example [52,53,54,55,56], and could be expected for continuously vegetatively propagated lavender varieties like var. Hemus. The performed BLAST search revealed that sequences of two NGS-SSR markers, 245 (LG4) and 247 (LG19), were allocated on two scaffold sequences of the reference genome: scaffold 293 (acc. JAPVEC010000093) and scaffold 515 (acc. JAPVEC010000026), which were not included into the chromosome sequence assemblies, Table S3. In reverse, the genotyping of the GEN-SSR marker S246 targeting an SSR region on scaffold 246 (acc. JAPVEC010000206) was affiliated to a locus on LG18, whereas the GEN-SSR marker S229 targeting an SSR region on scaffold 229 (acc. JAPVEC010000199) corresponded to loci on LG10 (S229A) and LG23 (S229B), Table S3. Accordingly, the established genetic map and additional genotyping of selected SSR markers could be applied to clear some ambiguities in the current reference chromosome assembly and affiliate some of the scaffolds to particular chromosomes. The use of high-density genetic maps for improving reference genome assemblies has been reported in other plant species. [57,58,59].

The BLAST search determined positions of SSR sequences in the reference genome and comparison with the respective SSR loci positions in the established genetic map shows that the order of SSR loci in several LGs well resembles those predicted from the corresponding chromosome sequence assemblies, for example LG6, LG9, LG15, LG18, LG19, LG20, LG22, LG24 and LG25, Table S3. For other LGs, the order of SSR loci has a patchy character in comparison with the reference genome predicted order, for example, LG7, LG10 and LG12, Table S3. These LGs contain one or more sub-linkage group(s) in which the loci order corresponds to the predicted from the reference genome, but the order of the rest of the loci differs from the predicted or the sub-LGs have twisted directions and/or positions in comparison to the predicted. The reasons for such differences are not clear and could be a result of chromosome rearrangements in the genome of var. Hemus or/and var. Munstead, local map disorders due to not sufficient marker saturation or even in local genome sequence disorders due to local irregularities of the genome sequence assembly. Further saturation of the map in the border regions of such sub-LGs will provide valuable information to clear these issues. The last together with the appearance of 16 loci with no distinct affiliation to a particular LG and the presence in some LG maps of adjacent markers with a large distance between them suggest that a follow-up genome guided regional saturation of the map will be essential for further improvement of the constructed map. The overall comparison of the established LG maps and corresponding positions of the markers in the reference chromosome sequences makes it possible to estimate the genome coverage of the map. The comparison shows that the genetic map covers between 49,2% and 99,4% of the reference chromosome sequences, with an average coverage of 84,4% +/- 14,2% per chromosome, Table 2. Correspondingly, the current level of marker saturation and genome coverage of the established map is sufficient to carry out efficient QTL analysis, as the early version of the map was already successfully applied for the identification of a QTL controlling the ratio of linalool to linalyl acetate in the flowers of L. angustifolia var. Hemus [60].

2.3. SSR Marker Transferability

Intra- and inter-species transferability of the employed SSR markers is essential for the wide range application of the established genetic linkage map and the construction of SSR marker anchored reference genetic map of lavender and Lavandula sp. To evaluate the transferability of the SSR markers applied in the established map, the SSR primer pairs were tested for PCR amplification of genomic DNAs isolated from L. angustifolia var. Hidcote Blue, two L. latifolia accessions, and the intersectional hybrid L. × heterophylla var. Big Boy James. A summary of the results is presented in Table 3 and described in details in Table S4, Supplementary Materials. The obtained results show high transferability of the majority of SSR markers used for the development of the genetic map, as above 93% of the tested SSR primer pairs show positive PCR amplification for the L. angustifolia garden variety Hidcote Blue, and above 78% of the SSR markers exhibit positive PCR amplification of DNAs from other Lavandula species. The observed high intra- and inter-species marker transferability suggests that the developed marker set could be efficiently employed for evaluation of the genetic diversity of natural populations and genetic resources collections of Lavandula species, as well as could be further applied for development of an integrated genetic map of L. angustifolia and SSR marker anchored genetic maps in the genus Lavandula. Such SSR marker-anchored genetic maps have been successfully and increasingly applied in modern plant breeding [61,62,63].

3. Materials and Methods

3.1. Plant Material

L. angustifolia var. Hemus plant used for the development of an F2 segregating population was obtained from the Institute of Roses and Aromatic Plants, Kazanlak, Bulgaria. For development of the segregating population, a single plant of var. Hemus was cultivated in the area of the village of Zagore, Thracian lowland, Bulgaria and seeds were collected after the 2020 flowering period. No other lavender plants were detected to grow in the same area during this period. The seeds were stratified for 12 weeks at 4 °C, and germinated in a tray on planting soil. The obtained seedlings were transferred in pots, grown in the greenhouse for 10 weeks and cultivated in the experimental field of the Agrobioinstitute (ABI), Kostinbrod, Bulgaria. The plants of L. angustifolia var. Hidcote Blue, L. latifolia and the intersectional hybrid L. × heterophylla ‘Big Boy James’ were purchased from Bastin Nursery (Kwekerij Bastin), Aalbeek, the Netherlands. The L. latifolia (Ll_abi2) plant was grown and selected from seeds purchased from WeberSeeds Botany Ethnobotany, the Netherlands. The plants were grown in pots in the greenhouse of ABI, Kostinbrod, Bulgaria

3.2. Genomic DNA Isolation

About 20–25 young leaves from every plant were placed in 15-ml plastic containers, immediately frozen in liquid nitrogen and stored at −80 ◦C. The frozen leaf samples were ground to a fine powder by using the Qiagen TissueLyser II Mill at 30Hz for 2min. Genomic DNA was purified according to the cetyltrimethylammonium bromide (CTAB) protocol [64]. Genomic DNA concentration was measured spectrophotometrically by using Nanodrop 2000 (Thermo Fisher Scientific, Waltham, MA, USA) and diluted to a final concentration of 25 ng/μL with Type I ultrapure water.

3.3. Next Generation Sequencing, SSR Identification and Primer Design

Microsatellite sequences from L. angustifolia var. Hemus were identified as a service by Ecogenics GmbH, Switzerland following NGS sequencing of a genomic DNA sample isolated from the var. Hemus plant used for development of the segregating population. The Illumina TruSeq nano DNA library was sequenced by using MiSeq Reagent Nano Kit v2 and Illumina MiSeq sequencing platform. The paired-end reads were passed through Illumina’s chastity filter and were further subject to de-multiplexing and trimming of Illumina adaptor residuals. The quality of the surviving reads was checked with FastQC v0.11.8 [65]. The paired-end reads were next quality filtered and merged with USEARCH v11.0.667 [66] to reform in silico the sequenced molecules. The merged reads were screened with the software Tandem Repeats Finder, v4.09 [67]. Merged reads containing SSR region with at least six repeat units of tri- or a tetra-nucleotide, or at least ten repeat units a dinucleotide of were selected. Primer 3 [68,69] software was used for primer design. Raw NGS sequences can be accessed at the NCBI Sequence Read Archive under accession number PRJNA1207064.

3.4. SSR Identification after Search of Reference Genome Sequence and Primer Design

The reference genome sequence of L. angustifolia var. Munstead was downloaded from the NCBI website (https://www.ncbi.nlm.nih.gov/datasets/genome/?taxon=39169, GenBank acc. No GCA_028984105.1) [37,70]. The genome sequence includes 795,075,733 bp chromosome-scale assembly representing 25 chromosomes with a N50 scaffold length of 31,371,815 bp [37]. The Krait Microsatellite Identification and Primer Design tool (v1.5.1) [71] was used for the search of reference genome sequence for the presence of SSR motifs. The following minimal parameters for SSR identification were applied: 12 for mononucleotides, 7 for dinucleotides, 5 for trinucleotides, 4 for tetranucleotides, 4 for pentanucleotides, and 4 for hexanucleotides. The primer pairs flanking the selected SSRs were designed to amplify predicted PCR products ranged from 100 bp to 300 bp. In addition to the primers for SSR regions identified from NGS and genome sequence data, five more primer pairs for PCR amplification of EST-SSR loci LAF1, LAF5, LAF8, LAF9 and LAF15, described by Adal et al. [41] were used within the study. All primers were synthesized by Macrogen Europe BV, the Netherlands. The primer sequences of primer pairs used for genotyping of plants from the segregating population are presented in Table S1, Supplementary Materials. The theoretical positions of the applied SSR markers were determined after BLAST search of the reference genome sequence of L. angustifolia var. Munstead by using the primer sequences and the sequence of the SSR genome region obtained from NGS data for var. Hemus.

3.5. PCR Amplification of SSR Regions

The PCR amplifications of SSR regions from plant genomic DNA, for testing the transferability of SSR markers or for analysis by the Agilent 5200 Fragment Analyzer System, were performed in a volume of 16 µL including: 1 µL of forward primer (10 pmol/ µL), 1 µL of reverse primer (10 pmol/ µL), 8 µL of 2x MyTaqTM Mix (Meridian Bioscience), 4,7 µL ultra-pure water, and 1.3 µL genomic DNA (25 ng/µL). The following PCR conditions were used: 95 ^oC for 3 min followed by 33 cycles of 95 ^oC for 15 s, 57 ^oC for 30 s, 72 ^oC for 30 s and a final elongation at 72 ^oC for 10 min. For testing the SSR marker transferability, the PCR products were resolved after electrophoresis in 1.8% agarose gels and observed under UV light. For fragment analysis and plant genotyping, the PCR products were subject to analysis by the Agilent 5200 Fragment Analyzer System (Agilent Technologies, Inc.) using the dsDNA 905 Reagent Kit (1-500bp).

3.6. PCR Amplification of SSR regions with Tailed Primers

The PCR amplifications of SSR regions from plant genomic DNA, subject of fragment analysis by a fluorescent capillary sequencer, were carried out using the designed reverse and 5′-tailed forward primers. Five different types of tails [72,73] were used including: M13 tail (5′-TAAAACGACGGCCAGT), A tail (5′- GCCTCCCTCGCGCCA), B tail (5′- GCCTTGCCAGCCCGC), C tail (5′-CAGGACCAGGCTACCGTG) and D tail (5′- CGGAGAGCCGAGAGGTG). The tails used for 5′-tailing of the forward primer were: M13 tail when the Tm of the reverse primer was ≤ 58 ^oC, C and D tails when the Tm of the reverse primer was 58 ^oC < Tm ≤62 ^oC, A and B tails when the Tm of the reverse primer was > 62 ^oC. The tails were synthesized and 5′-labeled with FAM, JOE, ROX and TAMRA dyes, as a service by Macrogen Europe BV, the Netherlands. All tails and primer combinations for the used primer pairs are presented in Table S1, Supplementary Materials. Based on the specific tail used, the annealing temperatures (Ta) of the PCR reactions were 54 ^oC for M13 tail, 57 ^oC for C and D tails, and 59 ^oC for A and B tails. The PCR reactions were performed in a volume of 16 µL, containing 0.8 µL of 5′-tailed forward primer (3 pmol/ µL), 1 µL of 5′-labeled tail primer (10 pmol/ µL), 1 µL of reverse primer (10 pmol/ µL), 8 µl 2x MyTaqTM Mix (Meridian Bioscience), 4 µL ultra-pure water, and 1.2 µL genomic DNA (25 ng/µL). The following PCR conditions were used: 95 ^oC for 3 min followed by 33 cycles of 95 ^oC for 15 s, Ta for 30 s (Ta is according to the type of tail used, described above), 72 ^oC for 30 s, and a final elongation at 72 ^oC for 10 min. The PCR amplified DNAs were further subjected to fragment analysis using ABI 3130 Genetic Analyzer (Thermo Fisher Scientific, Waltham, MA, US).

3.7. PCR Amplification of SRAP Fragments

In addition to the SSR marker genotyping, eleven SRAP primer pairs were used for amplification of SRAP fragments from the genomic DNA of L. angustifolia var. Hemus and plants of the analyzed segregating population. The ME and EM type of SRAP primers were designed according to Li and Quiros [74] and are shown in Table S2, Supplementary Materials. The PCR reactions were carried out in a volume of 20 µL containing: 1.5 µL of DNA template (25 ng/µL), 1 µL of 5′-labeled forward EM primer (10 pmol/ µL), 1 µL of reverse ME primer (10 pmol/µL), 10 µL 2x MyTaq HS Mix (Meridian Bioscience), 6.5 µL ultra-pure water. Samples were PCR amplified using the following thermal profile: 5 min at 94 ^oC; 3 cycles of 1 min at 94 ^oC, 1 min at 35 ^oC and 1 min at 72 ^oC; 35 cycles of 1 min at 94 ^oC, 1 min at 50 ^oC and 1 min at 72 ^oC and a final elongation step of 3 min at 72 ^oC. The forward ME primers were 5′-end labeled with FAM. The PCR amplified DNAs were further subjected to fragment analysis using ABI 3130 Genetic Analyzer (Thermo Fisher Scientific, Waltham, MA, US).

3.8. SSR and SRAP Fragment Analysis

ABI 3130 Genetic Analyzer with 36-cm long capillaries, Pop-7 polymer and GeneScan^TM 500 LIZ^TM size standard (all from Thermo Fisher Scientific, Waltham, MA, USA) were used for fragment analysis of the amplified SSR an SRAP fragments, as described earlier [42,50]. The fragment sizing was carried out by using GeneMapper 4.0 (Thermo Fisher Scientific, Waltham, MA, USA). The allele combination for each SSR locus and analyzed plants were determines following comparison with the alleles of the mother L. angustifolia var. Hemus plant. The SRAP fragment analysis was carried out for fragments in the range of 60–600 bp. All SRAP fragments showing a distinct presence or absence in the SRAP patterns of the analyzed plants from the segregating population were scored. Each of these SRAP fragments was considered a dominant allele of a separate locus, and the obtained data were used in linkage analysis.

3.5. Linkage Analysis and Genetic Map Construction

JoinMap 5 [75] was used for linkage analysis. The grouping was done by using the maximum likelihood (ML) mapping algorithm with the Haldane mapping function. The markers were grouped into linkage groups using a LOD score of 4.0 or higher. The values of Nearest Neighbor (NN) Fit and NN Stress from the Fit & Stress tabsheet and –log10P values from the Genotype Probabilities: Locus Means tabsheet were used to inspect quality of the established map. Markers related to poor map quality criteria were excluded and after that re-introduced one by one and mapping analysis was performed again. If re-introduction of a marker confirms its poor fit, the marker was removed from further analysis until a good quality map was established. The obtained order of SSR markers in each linkage group was compared with their theoretical position in the reference genome sequence. When the order of SSR markers in the linkage groups differed from the order of their theoretical position in the corresponding chromosome of the reference genome sequence the possible change of their order was tested through evaluation of the map quality after using the fixed order input in JoinMap 5 and observation of the map quality criteria (high NN Fit, NN Stress or –log10P values). When the fixed order didn’t result in deterioration of the map quality, it was considered as a genetic map of the corresponding linkage group. The genetic map were drawn with Map-Chart v2.32 [76].