ddRADseq Applications for Petunia × hybrida Clonal Line Breeding: Genotyping and Variant Identification for Target‐Specific Assays

1. Introduction

Petunia (Petunia × hybrida hort. ex E. Vilm.; synonym: Petunia × atkinsiana (Sweet) D. Don ex W. H. Baxter) is an annual herb belonging to the Solanaceae family, originated in South America [1] and well known worldwide for its ornamental use. A wholesale value of almost $ 160 million was reported in the US in 2020 for petunia ornamental varieties [2,3]. Commercial lines are derived from interspecific hybrids artificially obtained in the eighteenth century [4,5] between P. axillaris (Lamarck) Britton, Sterns and Poggemburg, and one or more species belonging the clade of P. integrifolia (Hooker) Schinz and Thellung [6,7]. Natural hybrids between the two species were never observed, in part due to different pollination syndromes. In particular, P. axillaris is pollinated by hawk moths, while P. integrifolia by bees [8,9]. Nevertheless, hybrids derived from artificial crosses can produce fertile progenies and commercial lines comprehend both seed-propagated F₁ hybrids and cutting-propagated clonal lines [1]. Although varieties are more often diploid, tetraploids have been observed and selected in breeding programs [10,11]. Petunia (2n = 2x = 14) is a facultative allogamous species, in which gametophytic self-incompatibility [12,13] and cytoplasmic male sterility [14,15,16] were extensively studied. It is considered a model plant [17,18], especially for the study of flower development and biochemistry [4], as well as for tissue culture practices [1]. In addition, the application of genetic engineering techniques has led to one of the first and well-known cases of genetically modified petunia, such as the orange-flowered A1-DFR lines [19]. It is also worth noting that the RNA interference mechanism was discovered for the first time in this species [20,21]. Finally, molecular resources developed for petunia comprehend high-density linkage maps with different marker types [22,23,24] and the draft genome assembly of the two of its main species P. axillaris and P. integrifolia [25]. Recently, a chromosome-scale genome assembly and annotation of P. × hybrida was developed by combining short- and long-read sequencing, optical mapping and Hi-C technology [26].

In the last decades, breeding goals for petunia have mainly included aesthetic-sensorial aspects, particularly in relation to flower morphology and phenology [2,27,28]. More recently, breeding has also dealt with new emerging needs, such as ecosystem service provision, adaptability and sustainability [29,30]. The evaluation of these traits through the exclusive use of morphological and phenological measurements is often impossible. Nevertheless, the official protocols for variety registration aimed at assessing distinctness, uniformity and stability (DUS), such as those of the European Community Plant Variety Office (CPVO) and International Union for the Protection of New Varieties of Plants (UPOV) [31,32], are still exclusively based on this type of descriptors. In petunia, 34 and 32 descriptors, including plant growth habit, leaf shape and flower density, are provided respectively by EU and US authorities for DUS tests [33,34]. Since official registration is necessary to achieve Plant Variety Protection (PVP) and to guarantee Plant Breeders' Rights (PBRs) [35], incorporating molecular data into the registered variety profiles, in conjunction with specific phenotypic measurements, could provide significant benefits. For example, a clear enhancement in a sustainability-related trait (e.g., tolerance to low nutrient levels or reduced phytosanitary inputs) assessed through appropriate phenotypic evaluations (e.g., value for cultivation and use (VCU) tests), may not be reflected by visible morphological changes. Hence, traditional DUS tests may be insufficient for a successful evaluation. In this case, DUS tests could be carried out coupling phenotypic measurements and genetic diversity analyses based on variety molecular profiles (value molecular DUS (vmDUS)) [36]. However, the possibility to exploit molecular genotyping for official variety registration is still limited to few specific functional markers, or it is used only to facilitate traditional DUS test procedures [37]. Molecular genotyping can be also useful for directly helping PVP, thus facilitating the identity assessment of plant materials of unknown origin and the identification of commercial frauds [35], as well as for Marker-Assisted Breeding (MAB) and Selection (MAS) [38,39,40].

Double digest Restriction-site Associated DNA sequencing (ddRADseq) is a Reduced Representation Sequencing (RRS) technology exploiting the gDNA digestion with two restriction enzymes, in order to reduce the fraction of the genome to be sequenced [41]. ddRADseq is reported as an approach able to couple reliable genotyping results to relatively low cost and time consumption. Therefore, it can be particularly useful for agricultural species lacking economic convenience in performing deeper and hence more expensive analyses [42,43,44], like the majority of ornamental plants. The technique was successfully exploited in several plant species for different purposes, such as structure evaluation of wild populations [45,46,47], germplasm conservation [48,49,50], Quantitative Trait Locus (QTL) identification [51,52], Genome-Wide Association Study (GWAS) [53], phylogenesis [54], Marker-Assisted Selection (MAS) and Marker-Assisted Breeding (MAB) [38,39,40]. The possibility to further reduce the thousands of variants generated by ddRAD down to an informative core (core set) may additionally simplify future analyses. In particular, SNP core sets may be used for the development of arrays or target-specific PCR (Genotyping By Target Sequencing – GBTS) [55,56,57,58,59]. Besides reducing costs and time of the procedures, this approach can lead to an easier comparison of data from different populations and/or years and allow targeting of specific regions of interest for further functional studies [60,61]. However, the identification of a SNP core set for target-specific assay development can be challenging, especially when dealing with breeding populations characterized by highly related lines with reduced genetic and phenotypic variability.

In this study, a genotyping analysis through a ddRADseq approach was performed in a petunia breeding population of clonal lines. The population included sister line groups, i.e. clones sharing at least one of the two parental lines. The goal was to assess the potential of this genotyping technique in uniquely discriminating each sample and determining the genetic relationships among the accessions considered. The suitability of ddRADseq for variety development through MAB, for variety registration and for protection of already registered varieties was evaluated too. In addition, we determined how much and how the SNP number could be reduced to reach an informative subset of markers (i.e. SNP core set), without altering the clustering of samples. Variants exclusively present in Coding DNA Sequences (CDSs) were selected to assess the possibility of exploiting expressed regions for the genotyping analysis.

2. Materials and Methods

2.1. Plant Materials

Leaf samples from 36 petunia clonal lines were granted by Gruppo Padana S.S. (Paese, TV, Italy), a nursery company specialized in the production, selection and collection of young plants from seeds and cuttings. The analyzed accessions included 36 breeding lines (BLs), including 4 benchmarks (BMs), i.e. already developed varieties exploited to evaluate the distinctness of the former. Genomic DNA (gDNA) was extracted from 100 mg of leaf tissue through DNeasy Plant mini kit (Qiagen, Valencia, CA, USA), following manufacturer instructions. Genomic DNA quality and concentration were assessed using NanoDrop 2000c UV-Vis spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and agarose gel electrophoresis with 1% agarose/ 1× TAE gel and 1× SybrSafe DNA stain (Life Technologies, Carlsbad, CA, USA).

2.2. ddRADseq Library Preparation and Sequencing

Double digestion Restriction-site Associated DNA sequencing (ddRADseq) libraries were produced following the protocol developed by Abed et al. [62], with some modifications. In detail, extracted gDNA was normalized to 20 ng/µl and a total of 200 ng per sample was used. PstI and MspI (Thermo Fisher Scientific) restriction enzymes were used, with 1µl of each per reaction, for DNA digestion at 37°C for 2 h, following the producer protocol. Adapter ligation was done in 50µl of reaction including 30µl of digested DNA, 1µl of T4 DNA ligase, 5µl of 10× T4 ligation buffer (New England BioLabs, Ipswich, MA, USA), 5µl of barcode adapters (0.1µM) and common adapters (10µM). Ligation was carried out with an incubation at room temperature for 2h, followed by 20 min at 65°C. Ligated DNA samples were grouped into two pools, each including 18 samples, and purified through QIAquick PCR and Gel Cleanup Kit (Qiagen) for the elimination of impurities, non-ligated adapters, small fragments and dimers. Size selection from 100 bp to 600 bp was performed with the CleanNGS magnetic beads (CleanNA, Waddinxveen, The Netherlands). Amplification of the obtained libraries was done with Platinum SuperFi PCR Master MIX (Thermo Fisher Scientific), following the enrichment PCR protocol described by Abed et al. [62]. Amplified libraries were subjected to another step of purification with CleanNGS magnetic beads. Subsequently, they were quantified using the TapeStation System (Agilent Technologies, Santa Clara, CA, USA) with High Sensitivity D1000 DNA ScreenTape assays (Agilent Technologies) and diluted to 200 pM for the Ion Torrent sequencing. Libraries were loaded into Ion Chips 550 (Thermo Fisher Scientific) through the Ion Chef (Thermo Fisher Scientific) automated machine. Single-end sequencing was performed in two runs, one per pool, the Ion GeneStudio S5 Plus System (Thermo Fisher Scientific) setting 200 bp of read length and 500 flows per run.

2.3. ddRADseq Data Processing

Raw reads were quality-checked filtering out those with a Phred Quality score lower than 20, and then demultiplexed according to the adopted barcode recognition motif, through the Torrent Suite platform (Thermo Fisher Scientific) following the native implemented pipeline. After these steps, reads were mapped on the reference P. × hybrida v.01 genome assembly retrieved from the NCBI database (GCA_046563445.1) [26]. Mapping was performed exploiting MEM algorithm of Burrows-Wheeler Aligner (BWA) v. 0.7.12 [63], followed by data conversion into BAM files through SAMtools and BCFtools v.1.13 software [64]. For the variant calling the gstacks and populations commands of Stacks v2.61 software [65] were used. gstacks was launched exploiting the ‘marukihigh’ model, with both a SNP discovery (--var-alpha) and a genotype call (--gt-alpha) alpha thresholds of 0.10, while populations was run in order to keep loci present in at least 30% of individuals in a subpopulation (-r 0.30) and 65% of individuals across subpopulations (-R 0.65). Obtained SNPs were filtered through the poppr v2.9.6 package [66] for RStudio software (Posit, Boston, USA), in order to retain only loci present in all samples and with a minimal Minor Allele Frequency (MAF) of 0.1. The resulting SNPs composed the initial dataset for the genetic analyses.

2.4. Genetic Diversity, Structure and Statistic Analyses

All analyses and graphic processing were carried out through RStudio, unless specified otherwise. Genetic Diversity (GD) matrixes were produced through the Nei’s coefficient [67] considering all pairwise comparisons among sample SNP profiles and were made through poppr package’s nei.dist function. A Kinship-based INference for Genome-wide association studies (KING) matrix based on pairwise comparisons was calculated with the KING-robust method exploiting the PLINK2 [68] --make-king function [69]. 0.4 was chosen as a minimal threshold to hypothesize the same clone membership for a pair. The Unweighted Pair Group Method with Arithmetic Mean (UPGMA) dendrograms for sample clustering were created based on Nei’s coefficient with the aboot function from the poppr package, considering 1000 bootstrap replicates. Clustering visualization was carried out through dendrograms created with plot.phylo function from the ape package [70]. Principal Component Analysis (PCA) based on the GD analysis was performed using the gl.pcoa function from the dartR package [71] Structure v.2.2 [72] software was used for Bayesian clustering to infer the sample genetic structure. The following parameters were set: founding groups ranging from 1 to 10, 10 replicate simulations per each K, burn-in of 20,000, final run of 100,000 Markov Chain Monte Carlo (MCMC) steps. The K values characterized by the highest likelihood were selected using Structure Harvester software [73]. Calculated genetic statistics included: observed (Ho), expected (Hs) and overall (Ht) Heterozygosity; inbreeding coefficient within subpopulations (Fis) and within the total population (Fit); subpopulation fixation index (Fst); gene flow (Nm); average number of observed (na) and expected (ne) alleles per locus; Private Alleles (PA) and Polymorphic Loci (PL) percentage, average locus Polymorphic Information Content (AvgPIC) and MAF. Statistics were measured either with hierfstat package [74] basic.stat function, dartR package loc.metrics function, or with a manually developed R script. Analysis of Molecular Variance (AMOVA) between and within subpopulations was done using the ade4 method with the poppr.amova function and the significance of results was validated through a Monte-Carlo test considering 999 permutations.

2.5. SNP Subset Selection for Core Set and Location in CDSs

Further filtering of the initial dataset for the SNP core set identification was done by setting alternative combinations of MAF, AvgPIC and Hs. The filtering continued till the obtained dataset led to a BL clustering consistently different from the one obtained using the entire SNP set. In parallel, Coding DNA Sequences (CDSs) were retrieved from the reference genome annotation and were intersected to the initial SNP dataset using the bedtools intersect command of Bedtools v.2.30.0 software [73], allowing to retain and use variants exclusively present in the coding sequences for clustering analyses. To visually represent them, a map was created using the R package CMplot [74].

3. Results

3.1. Genetic Diversity and Structure Analysis

A ddRADSeq analysis was performed to genotype 36 petunia BLs, comprehending four sister line groups, and 4 benchmarks (BMs). After quality filtering and mapping, more than 80 M reads with an average of 2.2 M reads per sample were obtained. This value corresponded to the 85.7% of the 12.3 Gb sequenced in total. The mean demultiplexed and trimmed read length resulted in 173 bp (Table 1; Supplementary Table S1).

Variant calling resulted in the identification of 42,236 SNPs located in a total of 15,655 loci. SNP filtering allowed the creation of a final dataset of 4,349 variant sites in the 36 analyzed samples.

A GD analysis among samples through the Nei’s coefficient resulted in no pairwise comparison with complete genetic identity, and unique molecular profiles were produced for all the considered BLs. GD among samples ranged from 0.001 (PeHy-BL3-06 vs. PeHy-BL3-08) to 0.51 (PeHy-BL3-09 vs. PeHy-BM03), with an average of 0.23 ± 0.11 (Supplementary Table S2). Considering sister line groups as a whole, the lowest GD value, 0.07, was found comparing PeHy-BL1 and PeHy-BL2 subpopulations, while the highest, 0.29, was observed in the comparison between PeHy-BL2 and PeHy-BL3 (Table 2).

The kinship matrix created with the KING-robust method led to an average value of 0.46 with a SD of 0.52. Three pairwise comparisons exceeded the 0.4 value, used as a threshold for the identifying of monozygotic twins or individuals belonging to the same clone [69,75]: PeHy-BL3-01 vs. PeHy-BL3-07; PeHy-BL3-06 vs. PeHy-BL3-08; PeHy-BL4-03 vs. PeHy-BL4-12 (Supplementary Table S3).

In the UPGMA dendrogram based on the GD analysis, three principal branches were observed: one grouping exclusively PeHy-BL3 lines, a second with PeHy-BL4 lines, and the last comprehending PeHy-BL1, PeHy-BL2 and the four BMs. All bootstrap values at the bases of the branches reached the maximum value (100). Within the last cluster, PeHy-BM3 was highly separated from the rest of the group, while a branch grouped together the other three BMs and two further branches hosted respectively PeHy-BL1 and PeHy-BL2 lines. All the relative nodes were characterized by very high bootstrap values (Figure 1A).

The genetic structure analysis showed the most likely K value, namely the most probable number of common ancestors, to be K=3 (Delta K = 10,407.3), followed by K=2 (Delta K = 3,068.5) (Supplementary Figure S1). In both cases, each BL was predominantly associated with a specific ancestral population, with the exception of PeHy-BL3-04 and PeHy-BL3-02, which resulted admixed. PeHy-BM3 and the PeHy-BL4-07 – PeHy-BL4-11 group were also found to be admixed, but only for one of the two K values (Figure 1B). When assuming K=2, PeHy-BL3 was the only group assigned to one of the two inferred ancestral populations, while all other BLs and the four BMs exhibited a shared ancestry. For K=3, in addition to the separate clustering of the PeHy-BL3 group, a further division was observed between the PeHy-BL4 group and the remaining lines (namely PeHy-BL1, PeHy-BL2, and the four benchmark populations) (Figure 1B). Results were in agreement with what was observed with the UPGMA clustering, and therefore with the genealogy data. The PCA analysis based on GD confirmed the findings observed in the Structure analysis and the UPGMA dendrogram. The first two components alone explained 50.3% of the total variability. PeHy-BL3 and PeHy-BL4 clustered apart, while the two other sister line groups and the benchmarks fell into the same group and were found to be poorly distinguishable (Figure 2).

3.2. Genetic Statistics and AMOVA

The sample pool was studied considering several genetic parameters. Average Ho per BL resulted in 26.0 ± 8.5%, with 1.26 na and 1.15 ne. The highest Ho values considering sister line groups were observed for PeHy-BL4 lines (33.6 ± 6.3%), while the lowest values were scored for PeHy-BL1 lines (19.1 ± 3.4 %) and PeHy-BM4 (16.3%). Average Fis was negative and comprehended between -0.2 and -0.3, underlining an overall moderate heterozygosity excess for BLs in relation to the expected for the respective subpopulations (Table 3; Supplementary Table S2).

Fst ranged from 0.22 to 0.50, suggesting a high genetic fixation for sister line groups (Table 2). As a result of these two factors, a slightly positive inbreeding coefficient for BLs in comparison to the expected considering the entire plant collection was observed (Fit: 0.02), in addition to a limited gene flow (Nm: 0.05). The percentage of PA per sample, relative to the total number of PA in the population ranged from 0.4% for PeHy-BL4-12 to 20.1% for PeHyBM3. The other three BM lines also exhibited percentages of private alleles considerably higher than their BL counterparts, ranging from 7.33% to 19.49% (Supplementary Table S2). The PA percentages considering each BL as a whole ranged from 6.4% for PeHy-BL1 to 34.0% for PeHy-BL3. A minor deviation was noted considering PL percentage of the total population loci: from 63.8% (PeHy-BL1) to 81.9% (PeHy-BL4) (Table 3). The AMOVA carried out comparing molecular variance between and within subpopulations (sister line groups and BMs) demonstrated good reliability, showing a p-value of 0.001. The analysis revealed that a substantially larger proportion of the total observed variability was attributable to differences among subpopulations, rather than among samples within the same subpopulation. Respectively, the estimated variance quotes were 58.2% and 41.8% (Table 4).

3.3. GD analysis Comparison between Total and Core Set SNP profiles

The 4,349 SNPs used for the GD analysis (Figure 3A) were filtered multiple times by acting on MAF, PIC and Hs values, to find the lowest number of variants (core set) able to uniquely identify and correctly cluster all the BLs.

The minimum SNP amount by which a unique molecular profile was conserved for all BLs was 1,051, which was obtained retaining only those variants with a MAF greater than 0.3. The UPGMA clustering and the bootstrap values remained unchanged for the majority of the previously highlighted groups (Figure 3B). Only the hierarchical position of the PeHy-BL1 cluster in relation to the PeHy-BL2 group and to the one including PeHy-BM1, PeHy-BM2 and PeHy-BM4, changed from the initial results (i.e. those obtained considering the whole SNP set) (Figure 3A). By increasing the MAF threshold to 0.45, the number of SNPs drops sharply to 103 with an average PIC of 0.3. Based on these settings, PeHy-BL3-06 and PeHy-BL3-08 presented complete genetic identity and minor clustering changes were observed (Supplementary Table S4, Figure 3C). The main difference found in the UPGMA grouping (compared to that obtained using the complete dataset) was a shift in the hierarchical positions of the PeHy-BM3 and PeHy-BL4 clusters. The bootstrap value at the base of the branch separating PeHy-BM3 from the rest of the accessions remained high (82.5), while the other basal nodes showed values lower than 50 (Figure 3C). Beyond the MAF threshold of 0.45, significant grouping distortions occurred. The same occurred when stringent filters related to PIC and Hs were applied to the original dataset (Figure 3D).

3.4. GD Analysis Comparison Between Total and CDS-Located SNP Profiles

The initial SNP set was intersected with the CDSs annotated in the reference genome in order to retain only those variants that were located within coding sequences. Genotyping and clustering the clones using exclusively these variants allowed to assess possible differences with the previous results derived from the same analyses. 1,975 SNPs were retained, corresponding to 45.4% of the total dataset. Their positions resulted uniformly distributed across the 7 chromosomes. Within each chromosome, distal areas presented higher variant site density in comparison to pericentromeric zones (Figure 4).

Pairwise GD comparisons led to values very similar to those observed considering the initial dataset (Avg: 0.22 ± 0.10) (Supplementary Table S5). The Ho, Hs and Fis values also remained largely unchanged. The UPGMA dendrogram led to a completely identical clustering to that obtained with the original dataset. Bootstrap values were slightly lower in the main nodes, but continued to be significant (> 80) (Supplementary Figure S2). The PCA explained 50.4% of the total genetic variation and subpopulation differentiation matched that obtained in the previous results. In particular, PeHy-BL3 and PeHy-BL4 groups were strongly separated, while PeHy-BL1, PeHy-BL2 and BMs samples resulted undistinguishable (Supplementary Figure S3).

4. Discussion

Molecular genotyping exploiting different types of molecular markers proved to be very useful for breeding plan setting and management, in several ways. It can help retrieving genetic information about populations, such as genetic diversity and heterozygosity, planning crosses, reducing inbreeding depression, exploiting heterosis, and favoring traits of interest’s selection [39,76,77,78,79]. In addition, the production of molecular profiles through genomic marker analysis can support PVP and variety registration, simplify DUS tests, help germplasm management, and aid the identification of GMOs, unknown plant materials and Essentially Derived Varieties (EDVs) [35,37,61,80,81]. In petunia, several genotyping studies have been accomplished taking advantage of different molecular marker typologies, including RAPDs [82], AFLPs [83,84], SSRs [8,85], SRAPs [86], and CAPSs [8,87]. GBS was used for genetic structure studies on petunia natural populations, including the parental species used for the production of the commercial interspecific hybrids, P. axillaris and P. integrifolia [88,89]. SNP genotyping was also successfully carried out on breeding collections for linkage map production and QTL analysis of flower-related traits, such as flowering timing, pollinator syndrome, and double flower development. Informative variants were retrieved from transcriptome data [7,90], targeted-sequencing techniques [2,22,91], and other GBS-based protocols [88,89,92,93]. By contrast, the use of molecular assays for the identification of commercial varieties and breeding lines, their registration, and protection remains an underexplored area in petunia, and in general in the ornamental plant sector.

The present study focused on molecular genotyping through ddRADseq of 36 petunia BLs and BMs. The genotyping analysis proved to be effective and useful for PVP-related purposes, assigning a unique molecular profile to all lines of interest and discriminating each genotype from the others. In addition, it allowed to estimate the GD values within the collection and other genetic parameters that are of interest for breeding plan management. BLs resulted subdivided into four sister line groups. As expected, sister lines with a relatively high homozygosity level were also characterized by a high genetic uniformity. Some lines (i.e. PeHy-BL3-01 vs. PeHy-BL3-07; PeHy-BL3-08 vs. PeHy-BL3-06; PeHy-BL4-03 vs. PeHy-BL4-12) turned out to be distinguished by a very limited polymorphism amount (respective GD: 1.7 × 10^-2; 1.2 × 10^-3; 0.29 × 10^-3) suggesting the possibility that they might represent the same clone. In this case, it is crucial to determine standard thresholds of molecular similarity above which two accessions should be regarded as belonging to the same clonal lineage. The exploiting of a kinship analysis based on the probability calculation that a pair of randomly sampled alleles is identical by descent (IBD) (i.e. due to genetic inheritance rather than by chance [94,95]) may achieve the standardization of genetic similarity values. Up to now, kinship matrixes in plant breeding studies have been exploited to identify the most genetically different genotypes in order to plan crosses [96] or to carry out the correction for relativeness in GWASs [97]. A threshold value commonly used for the identification of homozygous twins in animal genetics (0.4) [69,75] might be used to distinguish different clones. Indeed, all the three very genetically similar BL couples presented kinship coefficients higher than the threshold value. Nevertheless, the production of unique genotyping profiles for all lines demonstrated the approach capacity in discriminating also highly related genotypes. The SNP ability in uniquely genotyping very similar lines in petunia was previously shown by GBS studies for QTL analysis exploiting Recombinant Inbred Lines (RILs) and Near Isogenic Lines (NILs) [2,22,91].

The UPGMA clustering, the PCA and the genetic structure analysis led to consistent results and confirmed the genetic uniformity between sister lines, grouping them in the same branch, area of the plot, or ancestral group, respectively. The link with the genealogical data supported the reliability of the genotyping analysis. On the contrary, variability among sister line groups was significant. This was suggested by the AMOVA (variability between subpopulations > 50%), the fixation index calculation (Fst > 0.25 overall and in every sister line group pairwise comparison) and the genetic structure. However, despite the PCA explaining a percentage of total variability > 50%, groups PeHy-BL1, PeHy-BL2 and BMs resulted undistinguishable and were assigned to the same putative ancestor group. Therefore, although showing the highest presence of PAs, BMs resulted not entirely distinguishable at the genetic level from all the other lines. This implies that, in case a breeding aim is to obtain new varieties that are highly genetically distinguishable from those already present in the market, PeHy-BL1 and PeHy-BL2 lines did not match this goal. This statement highlights the importance in defining also genetic similarity thresholds, applicable at least at the species level, which allow to determine whether experimental materials present distinctness or should be considered as EDVs.

The collection under study showed a relatively low average level of Ho, especially if we consider its origin (breeding population of clones with a consistent structure) and its mode of reproduction (facultative allogamous species). This finding would suggest that the occurrence of crosses between closely related lines within the population may be possible, hence that the population would be at least partially self-compatible [98,99,100]. Nevertheless, inbreeding coefficient within subpopulations was generally low, highlighting a breeding selection effort to increase heterozygosis. The fact that the BL inbreeding coefficient within the entire population (Fit) did not considerably deviate from the equilibrium is the result of the high subpopulation differentiation, counteracted by low inbreeding coefficient within subpopulations. Heterozygosis levels (10-15%) even lower than those found in our study were reported in a genotyping analysis performed through RNAseq on petunia cultivars, mostly reported as F₁ hybrids. In the mentioned study, lines were representative of the total variation observable in the market considering the growth habit trait. Authors suggested a quite limited genetic diversity for the total commercial P. × hybrida germplasm pool [7].

The obtainment of a SNP core set to be used in the development of target-specific molecular assays for genotyping plant collections is notably important to facilitate the analysis of breeding populations and the data comparisons between years [55,56,57,58,59]. In our sample group, the reduction of the number of variants by maintaining a unique molecular profile for every BL was only possible by keeping a final set of more than one thousand SNPs. Also, in this case, clustering was not completely respected, thus making the dataset suitable for PVP-related purposes, but not ideal for MAB goals. The high genetic uniformity observed within sister lines may be the reason why a high minimum number of variants was needed to maintain their unique discrimination. Nevertheless, even in the case of breeding planning, there is often no need to distinguish very genetically similar genotypes as sister lines. This is because, in most cases, crosses are programmed between the most differentiated genotypes, in order to achieve a phenotypical variability or to combine multiple traits of interest. Considering this aspect, the variant reduction could be extended up to an SNP number around one hundred. In fact, with this variant set no significant clustering change (e.g. sister line separations) in comparison to the one-thousand SNP dataset was observable, and the hierarchical order of the four sister line groups remained unchanged. In addition, the non-discriminated genotypes only included a couple of sister lines previously hypothesized as belonging to the same clone. This number of variants might then be suitable for developing innovative target-specific genotyping assays. Studies aimed at selecting SNP core sets were done in collections of several species of agronomic interest, including lettuce, corn, barley, cotton, tobacco [61,101,102,103,104]. Core sets capable of discriminating plant accessions are generally characterized by only a few dozen markers. However, in the majority of the studies, only genetic diversity measurements were performed, without carrying out clustering analyses and hence not assessing the core set potential for MAB purposes. Furthermore, the accessions analyzed in the above-mentioned studies were part of germplasm collections, characterized by a genetic diversity far greater than that observed in our breeding population. This explains why a higher number of variants was required in our study. In addition, in the other works PIC is considered a key parameter for defining a core set of markers to discriminate between different accessions (e.g. [61]). However, in our study we were able to identify a suitable SNP dataset by primarily adjusting the MAF thresholds. Our results would suggest that MAF might be more efficient for the discrimination of strongly related lines (i.e. characterized by low genetic diversity), whereas it may be less useful than PIC if the genetic differentiation degree of the entire population is higher.

The selection of SNPs located in coding sequences led to retaining less than half of the initial dataset variants. This result is in agreement with what reported for ddRADseq analyses exploiting the same restriction enzymes for library preparation [105]. Indeed, PstI and MspI are methylation-sensitive and therefore tend to favour the selection of coding sequences with greater probability than methylation-unsensitive enzymes. The fact that CDS-located SNPs were more present in the distal zones of chromosomes can be due to a higher hypermethylation degree, which is typical of pericentromeric regions [106]. Indeed, centromere positions reported for P. × hybrida chromosomes [26] reflect the SNP lacking areas in all chromosomes. Variant density along chromosomes was also in agreement with what already reported for ddRADseq analyses with methylation-sensitive enzymes [50]. Coding sequences are more conserved than the non-coding ones and hence less polymorphic [107]. However, the exclusive use of SNPs located in CDS regions led to results that were entirely comparable to those obtained with the initial dataset (considering the GD, UPGMA clustering, PCA, and Structure analysis) and to the genealogy data (i.e. similar genotyping outputs for sister lines). The use of coding sequences for genotyping analyses may enhance the identification of variants linked to traits of interest and the results obtained in this study might represent a clue for the reliability of this approach. In addition, the high number of variants located in CDSs retrieved (1,975) may suggest ddRAD as a potential technique for the identification of loci to be used for functional gene variant identification.

5. Conclusions

In conclusion, ddRADseq was demonstrated to be effective in uniquely discriminating petunia experimental clonal lines, in assessing their genetic relationship and in providing informative data for breeding purposes. In addition, the selection of a SNP core set was achieved, in order to allow the development of target-specific methods leading to simplified and more comparable genotyping results. Furthermore, the adopted method demonstrated to be exploitable in detecting polymorphisms within coding sequences, which may provide suitable variants for targeted in-depth analyses to assess their putative involvement in the traits of interest’s manifestation. Results encourage the adoption of GBS techniques for PVP to guarantee PBRs, for the registration of novel varieties whether coupled with DUS testing and other phenotype measurements, for MAB approaches and for the development of targeted-specific assays. The production of molecular profiles relative to experimental plant materials can help the achievement of innovative traits for cultivars, especially related to sustainability cultivation and management, which is aimed at reducing the agronomic input needs of plants and therefore enhancing the economic and environmental security. Further studies are needed toward two goal directions: the enhancement of basic molecular knowledge relative to petunia, and the development of useful applications for breeders, as assays able to predict traits of interest’s values. Eventually, GBS analyses need to be extended to understudied species, as it is for the majority of the ornamental plants nowadays commercialised worldwide.