A total of 123 gerbera accessions were collected from different province of China and were grown under identical natural conditions at South China Agricultural University, Guangzhou, China (23.16°N, 113.36°E). The petals of cultivars were collected during the full flowering stage. Three to four biologic replicates of each variety were obtained. For pigment analysis, the collected samples were immediately frozen in liquid nitrogen and stored at –80°C. The detailed information of accessions is listed in Supplementary Table 1.

At full boom, the colors of petals were described using the RHSCC (Royal Horticultural Society Color Chart) with white background and under the same lighting conditions. Furthermore, a CS-210 precision colorimeter (Hangzhou Caipu Technology, China) was used to quantify the color parameters of the same flowers. The parameters of CIEL*a*b* color coordinate (lightness, L*; chromatic components, a* and b*) were also measured. According to the formula: h=arctan(b*/a*),C=(a+*2b)*21/2, the hue angle h and chroma C were measured (Gonnet, 1998). As replicates, five measurements were taken for each accession.

Total carotenoid (Ct) content was measured as described previously (Singh et al., 2015). Briefly, 2 g fresh flower samples were ground in 10 mL of acetone followed by centrifugation at 5,000 rpm for 8 min. The supernatant was repetitively ground and filtered till devoid of the color. The extract was shifted to a 50 mL flask containing 15 g of anhydrous sodium sulfate. Total carotenoid content was measured using the equation: total carotenoid contents (μg/g) = [(A×V(ml)× 10⁴)/A%cm×P(g)]; A = absorbance, V = total extract volume, P = fresh sample weight, A% cm = 2,592 (extinction coefficient). The experiment was performed in three to five replicates.

Total anthocyanin (At) content was measured as explained previously with minor modifications (Ai et al., 2016). In short, ∼300 mg of flower petals per sample was collected and ground in liquid nitrogen into a fine powder. To extract anthocyanins, the fine ground powder was transferred to a 5 mL extraction solution comprising hydrochloric acid and methanol mixture (1:99 v/v). Thereafter, the extract solution was incubated at 4°C for 24 h followed by centrifugation at 13,000 rpm for 20 min. Then, the supernatant was shifted to a new tube, and anthocyanin was measured by calculating the optical density (OD) at A₅₃₀ by using a UV-1200 spectrophotometer (MAPADA, Shanghai, China). The measurement of At was conducted using the equation Q_At = A₅₃₀ × M⁻¹; Q_At = amount of At, M = fresh weight (g) of the plant material used for extraction. Data for each sample are reported as the mean of three replicates.

Total RNA was isolated using the hexadecyl trimethyl ammonium bromide (CTAB) method with minor modification as described earlier (Tel-Zur et al., 1999; Yue et al., 2014). The quality and quantity of RNA were assessed using the Agilent Bioanalyzer (Agilent Technologies, United States). Furthermore, the integrity of total RNA was evaluated via agarose gel electrophoresis. For RNA-sequencing, the cDNA library was generated using the Illumina Hicseq™ RNA sample prep kit (Illumina, United States) following the manufacturer’s protocols. The size and concentration of the library were evaluated using Qubit 2.0 and Agilent 2100. Through Hi-Seq 2500 sequencing machine, pair-end (150 PE) Illumina high-throughput sequencing was conducted. The transcriptome raw data of gerbera have been deposited in the NCBI database (SRR16889011).

From the raw reads, the low-quality reads and adopters were filtered. The clean reads were de novo assembled into contigs with an optimized k-mer length = 25 and group pairs distance = 300 using the Trinity program2. The unigenes functions were predicted via BLAST against the NCBI non-redundant protein (Nr), NCBI nucleotide sequences (Nt), and Swiss-Prot databases (E-value of 10^–5). The resulting datasets were validated to the Protein family (Pfam) database with HMMER (E-value 10^–10). Unigene sequences were aligned against databases of the Gene Ontology (GO) (Ashburner et al., 2000), Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2004), euKaryotic Orthologous Groups (Koonin et al., 2004).

From the transcriptome data, SSR motifs were identified using the microsatellite identification online program (Beier et al., 2017). The identification standards for dinucleotides (di-), trinucleotides (tri-), tetranucleotide (tetra-), pentanucleotide (penta-), and hexanucleotide (hexa-) motifs were set to six, five, five, five, and five repeats, respectively. The SSR primers were designed using the Primer 3 software (Koressaar and Remm, 2007). The designed primer pairs were assessed using in silico PCR with unigene sequences as the template. The primers with several amplicons were screened out.

One hundred and thirty-nine EST-SSR markers were used for genotyping 123 gerbera accessions. Genomic DNA was obtained from the fresh and young leaves of each cultivar using the CTAB procedure. DNA quantification was performed via a Nanodrop ND-1000 spectrophotometer (Thermo Scientific, United States). T100™ Thermal Cycler (BIO-RAD, United States) was used for PCR and the conditions were the same as described previously (Zhou et al., 2021). The amplicons were separated using 8% polyacrylamide gel electrophoresis at 300 V in 1.0 × TBE buffer and visualized with silver staining.

SPSS 19.0 program (SPSS Inc., Chicago, IL, United States) was used for statistical analysis and Student’s t-test. Correlation analysis was based on Spearman correlation coefficient using R 4.0.5. PLS-DA was implemented with the “mixOmics” package of R (4.0.5.). The regression coefficient network was performed using Cytoscape version 3.6.2 (Cline et al., 2007).

The number of alleles (Na), number of effective alleles (Ne), Shannon’s index (I), observed heterozygosity; (Ho), expected heterozygosity (He), and the polymorphic information content (PIC) were calculated using GenAlEx 6.5 software (Peakall and Smouse, 2006) and PowerMarker version 3.25 (Liu and Muse, 2005). Based on the admixture model, population structure analysis was performed using STRUCTURE 2.3.4 (Pritchard et al., 2000) and web-based software STRUCTURE HARVESTER version 0.6.92 (Earl, 2012). CLUMPP v1.1.275 was used to assess the 20 replicates for optimal alignment cluster (Jakobsson and Rosenberg, 2007). To cluster the genotypes, the neighbor-joining (NJ) method was employed based on Nei’s genetic distance in PowerMarker software. Principal coordinate analysis (PCoA) was employed using the “vegan” and “ape” package in R 4.0.5.

Association analysis was performed using the mixed linear model (MLM, Q + K model) in TASSEL4.3.6 (Bradbury et al., 2007). The significance of marker-trait associations was determined at P ≤ 0.01 (Jaiswal et al., 2012). To distinguish the true associations, False discovery rate (FDR) was employed (Benjamini and Yekutieli, 2005).

The flower colors of the 123 gerbera accessions ranged from white to purple. The HCA and RHSCC analysis showed that 123 gerbera accessions could be classified into six color groups including white, yellow, orange, pink, red, and purple group (Figure 1 and Supplementary Table 1). Among them, 32 accessions were identified as red group which accounted largest proportion (26.02%) of total gerbera accessions. Similarly, the white and orange groups have 20 gerbera accessions and contributed 16.26% of total gerbera accessions followed by pink (19 accessions, 15.45%), purple (18 accessions, 14.63%), and yellow group (14, 11.38%), respectively (Figure 1).

The L* values of the white and yellow groups were the highest, and the degree of variations among varieties was small, while the L* values of the red and purple groups were the lowest (Figure 2A). The a* value of the white and yellow groups was the smallest and there were negative values in nine accessions, while other four-color groups were distributed in the positive range (Figure 2B). Meanwhile, the red group had the highest a* value with an average of 71.12, followed by the purple, orange, pink, and white group, respectively. The purple group had the smallest b* value, and all of them were less than zero, whereas the b* value of the other color groups was greater than zero (Figure 2C). The b* value of the yellow group was the highest with an average of 58.14, followed by the orange, red, white, and pink groups, respectively. In addition, the C value of the red group was highest followed by orange, purple, yellow, pink, and white, respectively. Moreover, the C value of the white group was the lowest with an average value of 13.14 (Figure 2D). The h value of the purple group was the highest with an average of 358.58, followed by yellow, white, orange, pink, and red, respectively (Figure 2E). The spatial distributions of these six-color groups based on CIEL*a*b* color coordinates were different (Figure 2F) supporting the result of HCA analysis (Figure 1).

The Ct content in the orange group were significantly higher compared to other groups with an average of 10.27 μg g^–1 FW, followed by yellow (4.72 μg g^–1 FW) and red groups (3.89 μg g^–1 FW), respectively (Figure 3). Whereas, Ct contents in white, purple, and pink groups were low. Furthermore, At content in the red group were highest (49.85 U g^–1 FW), followed by purple (24.54 U g^–1 FW), orange (10.83 U g^–1 FW), and pink groups (10.05 U g^–1 FW). Likewise, At content in the white and yellow groups were the least (Figure 3).

The correlation analysis showed that seven variables had different degrees of correlation, ranging from -0.9624 to 0.9156 (Figure 4A and Supplementary Table 2). L* value negatively correlated with a* value, C value and At content. Moreover, a* value negatively correlated with b* and h value, while positively correlated with C value and At content; b* value positively correlated with C value and Ct contents, while negatively correlated with At contents; C value negatively correlated with h value, while positively correlated with Ct and At contents; h value negatively correlated with At content. The analysis showed that b* and C value were mainly influenced by Ct content, while L* and a* value were mainly influenced by At content.

To analyze the effects of different chromaticity parameters (L*, a*, b*, C, and h) and pigments (Ct and At) on different color groups, a PLS-DA analysis was performed. The cut-off graph shows that seven variables are between two confidence ellipses, and the score graph shows that each color group is distinguished, which indicates that the constructed model can well explain the differences of different color groups (Figures 4B,C). From each color group, the white group was mainly affected by L* value, a* value, C value, and At content with regression coefficients of 0.61, –0.59, –0.65, and –0.50, respectively (Figure 4D and Supplementary Table 3). The yellow group positively correlated with L* value, b* value, and At content, while negatively correlated with a* value, h value, and Ct content, while the absolute value of the regression coefficient of C value is low. In the orange group, b* value and total Ct content were the main factors, with regression coefficients of 0.55 and 0.47, respectively. The pink group was positively correlated with the L* value and h value, and the regression coefficients of all variables were less than 0.3. The red group was mainly affected by L*, a*, C, and At content, with regression coefficients of –0.59, 0.58, 0.57, and 0.50, respectively. Moreover, the red group positively correlated with Ct content. The purple group was mainly affected by L*, a*, b* and At content, with regression coefficients of –0.45, 0.46, –0.55, and 0.43, respectively. Furthermore, the purple group had a negative correlation with Ct content.

After filtering of the raw sequencing data, 57,133,428 clean reads (∼8.57 Gb) were obtained from the transcriptome. RNA-Seq statistics for the tissue are listed in Supplementary Table 4. More than 94.49% of reads had Q30 or higher quality scores. Using the Trinity de novo assembly program, the ∼8.57 Gb high-quality reads were assembled into contigs, transcripts, and unigenes. The statistical data of the assemblies are shown in Supplementary Tables 5, 6. A total of 182,868 transcripts (mean length 486 bp) representing 111,709 unigenes were assembled (mean length 870 bp). The number of transcripts between 200 and 500 bp was the highest, accounting for 51.05%, while the number of genes between 500 and 1,000 bp was the lowest (32.88%).

Searches against the GO, KEGG, KOG, Nr, Nt, and Swiss-Prot databases yielded a total of 111,709 annotated unigenes. 11,295 (10.11%) of them had hits in all seven databases. The unigene hits in each database are shown in Supplementary Table 7 with the number of annotated unigenes being the highest in the Nr (81,057) and the lowest in KOG (21,279). Further BLAST searches against other databases showed that 48,034, 31,484, 60,911, 56,206, and 56,897 unigenes had at least one match in the Nt, KEGG, Swiss-Prot, Pfam, and GO databases, respectively (Supplementary Table 7). In BLASTx homology searches (at cutoff E-value of 10^–5), the five top species with hits in the Nr database were: Vitis vinifera (12,848, 16.03%), Sesamum indicum (6,821, 8.51%), Coffea canephora (6,106; 7.62%), Nicotiana sylvestris (4,078, 5.09%) and Nicotiana tomentosiformis (3,796, 4.74%) (Figure 5A).

In KOG analysis, 21,279 unigenes were classified into 25 clusters (Figure 5B). Among the 25 KOG categories, the cluster related to general function prediction (2,859 or 13.44%) was the largest, followed by the clusters of posttranslational modification, protein turnover, chaperones (2,671 or 12.55%) and translation, ribosomal structure, and biogenesis (1,743 or 8.19%). The clusters represented by the least number of unigenes were cell motility and extracellular structures (40, 0.19%).

Kyoto Encyclopedia of Genes and Genomes pathway analysis of the 31,484 unigenes revealed that 23,584 (74.91%) could be assigned to 19 pathways belonging to five major categories: metabolism (13,425), genetic information processing (6,261), cellular processes (1,599), environmental information processing (1,227), and organismal systems (1,072) (Figure 5C). The category with the largest number of unigenes was metabolism, in which the most represented four pathways that were carbon metabolism (ko01200) (987, 7.35%), starch and sucrose metabolism (ko00500) (888, 6.61%), biosynthesis of amino acids (ko01230) (832, 6.20%), and purine metabolism (ko00230) (761, 5.67%).

In addition, 56,897 unigenes that matched the GO database were classified into 56 functional sub-groups of the three main GO groups: biological process, cellular component and molecular function (Figure 5D). The majority of the unigenes were assigned to biological processes (140,724, 47.56%), followed by cellular components (86,093, 29.09%), and molecular functions (69,097, 23.35%). Under the category of biological processes, cellular processes (32,193, 22.88%) and metabolic processes (30,201, 21.46%) were predominant. In the cellular component category, cell part (16,753, 19.46%) and cell (16,749, 19.45%) were the most abundant classes. As for the molecular function, catalytic activities (32,996, 47.75%) and binding (26,246, 37.98%) were the top two categories in numbers.

A total number of 14,106 SSR loci with motif types were discovered. Di- repeat motifs (8,452, 59.92%) were predominated among the total SSRs identified, followed by Tri- repeat motifs (5,333, 37.81%), Tetra repeat motifs (138, 0.98%), Hexa- repeat motifs (128, 0.91%), and Penta- repeat motifs (55, 0.39%). For the Di- and Tri- repeat motifs, the frequency was as high as 97.73% in total. The number of motifs types for Di-, Tri-, Tetra-, Penta-, and Hexa- were 8, 30, 34, 18, and 67, respectively (Supplementary Table 8).

The number of SSR loci gradually decreased as the number of repeat times increased (Table 1). Among them, SSR loci with five repetitions were the foremost. The most abundant Di- motif was GA/TC, responsible for 25.21% of Di- motif repeats followed by CA/TG (19.46%). In Tri- repeats, GAA/TTC (10.97%) was the most abundant, followed by AAT/AAT accounting for 10.46%. In Tetra- repeats, TATG/CATA (13.77%) has the largest number and AAAT/ATTT (10.14%) was the second, while the other motif accounts for less than 7%. In Penta- repeats, AATCG/CGATT (14.55%) has the largest proportion, followed by AAAAC/GTTTT (12.73%) and GCTTA/TAAGC (10.91%). In Hexa- repeats, the most abundant motif types are AACCAA/TTGGTT, accounting for 7.03%, while the number of other motifs types was less than 30 (Figures 6A–E and Supplementary Table 9). The five types of motifs with a number of more than 1,000 are shown in Figure 6. The most copious motif (>500) was GA/TC (2,068, 14.66%) followed by GA/TG (1,596, 11.31%), AG/CT (1,470, 10.42%), AC/GT (1,070, 7.59%), AT/AT (1,003, 7.11%), TA/TA (939, 6.66%), GAA/TTC (585, 4.15%) and AAT/ATT (558, 3.96%), respectively. The total number of the remaining 149 motif types was 4,567 only (Figure 6F and Supplementary Table 9).

To develop polymorphic SSR markers suitable for 123 accessions, 139 pairs of SSR primers were selected, including 53 previously developed EST-SSR primers and 86 newly developed EST-SSR primers. We used a PAGE analysis to screen out 48 pairs of SSR primers with a good amplification effect and polymorphism to analyze 123 gerbera accessions. Among 48 SSR loci, 3 are Di-, 28 are Tri-, 4 are Tetra-, 3 are Penta-, 8 are Hexa-, and 2 are complex repetition types with lengths ranging from 12 to 54 bp (Supplementary Table 10). The 48 pairs of primers can amplify 155 alleles in the test accessions (Table 2). The PIC values ranged from 0.20 to 0.77 with an average of 0.44, indicating that these plants showed middle genetic diversity. The number of observed alleles (Na) ranged from 2 to 7 with an average allele number of 3.21. The number of effective alleles (Ne) ranged from 1.26 to 5.03 with an average allele number of 2.25. The Shannon information index (I) ranged from 0.38 to 1.76 with an average of 0.86. The observed heterozygosity (Ho) ranged from 0.03 to 0.81 with an average of 0.45. The expected heterozygosity (He) ranged from 0.21 to 0.80 with an average of 0.50. Among 48 SSR primers, GEM65 yielded the highest number of alleles (7 alleles), with a PIC of 0.77.

Population structure analysis (PSA) revealed that when K = 2, ΔK was the largest, and 123 accessions were divided into two clusters (Figures 7A,C) and ΔK was the second largest when K = 3. NJ cluster analysis and principal coordinate analysis (PCoA) also revealed that there were two major clusters among the gerbera population (Figures 7D,E). For each Q-value, genotypes with membership probability >60% were assigned to the same group, while those with <60% probability in any group were assigned as “admixed” (Lu et al., 2009). When K = 2, cluster I include 101 accessions with an average Q value of 0.8516 (72.28% of them with Q > 0.80). Cluster II contains 22 accessions with an average Q value of 0.8659 (72.73% of them with Q > 0.80). The Q value of six accessions was less than 0.6, indicating that they were admixed in which four accessions (G006, G081, G108, and G122) belong to cluster I and 2 accessions (G052 and G099) belong to cluster II. There was no significant difference in petal color group among different clusters for flower color characters (Figures 7C,E and Supplementary Figure 1). However, in cluster I, the flower core color of 58 accessions were dark, while 45 accessions were light (Supplementary Table 1). In cluster II, 86.36% of accessions’ flower core were dark-colored, and only three accessions’ flower core (G099, G106, and G119) were light-colored (Supplementary Table 1). In terms of geographical origin, cluster I include 50 accessions from Yunnan, six accessions from Shaanxi, and 47 accessions from Guangdong province (Figure 7B and Supplementary Table 1). In cluster II, accessions from Guangdong province accounted for 86.36% and only three accessions were collected from Yunnan (Figure 7B and Supplementary Table 1).

Neighbor-joining cluster analysis was consistent with the population structure analysis (Figure 7D). The analysis showed that 123 accessions in the NJ tree were divided into two groups. Group A includes all accessions in cluster I and one admixed accession (G106) in cluster II, while Group B includes 21 accessions in cluster II. The two clusters were also supported by PCoA with 48 EST-SSR markers (Figure 7E). The boundary between the two clusters was obvious. The total proportions of the variation explained by the first and second principal components were 6.48 and 4.74%, respectively.

An association analysis was performed to further obtain the color-related molecular markers of gerbera and to validate the application potential of 48 SSR markers. A mixed linear model (MLM, Q + K model) was used to analyze the correlation of Ct and At contents. The findings are shown in Table 3 and Figure 8. The number of markers significantly associated with L*, a*, b*, C, and h values were 1, 1, 2, 1, and 1, respectively with an average R² of 0.08673. Two markers were significantly correlated with Ct content with an average R² of 0.06206, while one marker was significantly correlated with At contents. Interestingly, some markers can be associated with multiple traits. GEM65c was significantly associated with L*, a*, C, and At, while Gem148c and P3-28b were significantly associated with b* value and Ct, respectively. When using a higher threshold standard (FDR = 0.1), only GEM65c was significantly associated with L* (R² = 0.1143) and At (R² = 0.1034). It can be seen from Figure 8F that the L* value of the non-banded group was significantly higher than that of a banded group, and Figure 8H shows that the At content of the non-banded group was significantly lower than that of a banded group.