A set of 151 accessions from different geographical origins ( Table S1 ), was used for the GWAS analysis, and was provided by the cucumber research group at the Institute of Vegetables and Flowers, at the Chinese Academy of Agricultural Sciences in China. The 151 diverse plants genotypes were grown, each line includes 5 individuals and self-pollinated under natural field conditions in the Institute of Vegetables and Flowers, Chinese Academy of Agricultural Science at Beijing in 2018 (39°90N, 116°30E; average day/night temperatures 26/14°C; day length ~14 h in May and June). The individuals of each line were collected, and mixed seeds were further used for phenotype identification.

Accessions of the 151 cucumber core germplasms were phenotyped in two different environments: a phytotron in Beijing in summer of 2019 (2019S) and an incubator in Shandong (36°51N, 118°50E) in Shandong province in summer of 2020 (2020S) ( Table S1 ). Fifteen seeds per triplicate were placed at 13°C for 14 days in the dark and at 28 °C for 7 days as control (Kłosińska et al., 2013). For each experiment, three replicates were applied. Each replicate was randomly placed and the core germplasm collection was randomly placed in each replicate (block). Fifteen seeds per line were placed in a 90-mm-diameter Petri dish containing two layers of filter paper soaked 6 h with 30 ml 55°C water, then were exposed to 6 ml of distilled water. Seed germination (radicle ≥ 2 mm long) was recorded daily during the treatment stage (14 days for low temperature and 7 days for control) and radical length were measured at 14 day for low temperature and 7 day for control by Image J (Song et al., 2018). The relative value of four evaluated indices, including relative germination rate (RGR), relative germination energy (RGE), relative germination index (RGI) and relative radical length (RRL) were calculated between 13°C and 28°C (Song et al., 2018). The four indices, i.e., RGR, RGE, RGI and RRL were used for GWAS analysis.

Statistical analysis of the phenotypic data was implemented using SAS v.9.4 (SAS Institute Inc., Cary, NC, USA, 2012). Pearson’s liner correlation coefficients were found among the four indices, and were analyzed using R software (R Core Team, 2019). An analysis of variance (ANOVA), principal component analysis (PCA) analysis was also performed in SAS v.9.4 (SAS Institute Inc., Cary, NC, USA, 2012). Membership function and comprehensive evaluation values (D-value) were analyzed in Excel 2016 (Xie, 1983; He, 2006). Principal component W_j = I_j/∑I_j (j = 1, 2, …, n), where I_j represented the contribution rate of the jth principal component; Membership function value (X_j) = (X_j - X_min)/(X_max-X_min) (j = 1, 2, …, n), where X_j represented the jth principal component value, X_min and X_max represented the minimum and maximum value of the jth principal component in different germplasms, respectively; Comprehensive evaluation values D = ∑ (U_j × W_j) (j = 1, 2, …, n). The D-value as a comprehensive evaluation index was used to represent the LT resistance of various germplasm. The D-value was used as an input for cluster and GWAS analysis. The sum of squares of dispersion (Ward) method was used for cluster analysis and the Euclidean distance method was used to measure in R software (R Core Team, 2019) based on the average D-value of each core germplasms in two experiments.

The re-sequenced genome data and SNPs of the 151 accessions of the core germplasm collection are available in the NCBI Database (PRJNA831637). A FaST-LMM model (Jiang et al., 2021) was used for association tests of LT tolerance with an estimated relatedness matrix as covariate. All of the 151 genotypes were included in the GWAS. P=1.0 × 10⁻⁶ was used as the significant threshold. The linkage disequilibrium (LD) decay coefficients (r² ) between high-density SNPs was calculated using Plink software (Purcell et al., 2007), and was used to evaluate LD decay (Liu et al., 2020b). Association mapping of LTG traits using SNP genotyping, population structure, relative kinship and LD analysis were identified (unpublished). Manhattan plots and the LD heatmaps were drawn by qqman package (Turner, 2014) in the R environment.

The candidate locus was identified by repeatedly detected loci of various indices in two environments. The regions with significant trait-associated SNPs were described as candidate regions by unifying a set of various indices. The annotations of candidate genes in QTL regions were obtained based on “cucumber genome V3 (‘9930’ as reference genome) (http://cucurbitgenomics.org/organism/20)” (Li et al., 2019). The candidate genes were identified by the following steps: first, many genes related to abiotic stresses were predicted using Swiss-Prot (https://www.uniprot.org/)、TAIR (https://www.arabidopsis.org/index.jsp) and the gene ontology (GO) database (http://amigo.geneontology.org/amigo/landing). Secondly, SNPs in the exons or promoters of different genes, that differentiated the 47 LT-sensitive lines from the 21 LT-resistant lines were identified based on sequence alignment using MEGA X (Kumar et al., 2018). Finally, the expression level of these candidate gene was analyzed at varying time points, between germplasm differing in their LTG response, to further identify candidate genes.

Promoter (the 2.0 kb upstream of ATG) were used for searching regulatory elements of the CsPPR gene using PlantCARE online (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/) (Lescot et al., 2002). And the result was visualized by TBtools (Chen et al., 2020). Cucumber orthologs of the Arabidopsis PPR proteins were identified using the Cucurbit (http://cucurbitgenomics.org/feature/gene/CsaV3_1G044080) and Arabidopsis databases (https://www.arabidopsis.org/servlets/TairObject?type=locus&name=AT4G04370). The functional interaction networks of proteins were identified by the online STRING database (https://cn.string-db.org/cgi/network?taskId=baqwil7ulR5D&sessionId=bZBjkvzPG61G).

A time-course expression analysis of the four sensitive lines (‘R13’, ‘R42’, ‘R77’, ‘R137’), together with the four resistant lines (‘R99’, ‘R152’, ‘R167’, ‘R174’) were examined during seed germination. Cucumber seeds exposed to 13°C in an incubator for 0 d,12 h, 1 d, 3 d and 7 days were collected. There were three biological replicates for each treatment. For Arabidopsis thaliana, Columbia-0 ecotype (wild-type) Col-0 and five T₁ transgenic lines were grown 1/2 MS with hygromycin agar medium. Ten-days-old seedlings were collected for RNA extraction.

Total RNA was extracted, and the first-strand complementary DNA (cDNA) was synthesized using RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) and a PrimeScript RT Reagent Kit with gDNA Eraser (TaKaRa, Kyoto, Japan), respectively, according to the manufacturer’s instructions. Quantitative reverse transcription PCR analysis (qRT-PCR) was performed using SYBR Premix Ex TaqTM II (TaKaRa, Kyoto, Japan). The PCR primers ( Table S2 ) were designed with DNAMAN7.0 (Woffelman, 2004). Actin1 (Csa3G806800) was employed as the reference gene to normalize gene expression values (Xie et al., 2018). Each qRT-PCR experiment was performed with three biological replicates. The analysis of gene relative expression data was performed using the 2^−ΔΔCt method (Livak and Schmittgen, 2001). Significant differences between different materials were detected using two-tailed, two-sample Student’s t-test (p<0.05 or p<0.01) in Excel (2016).

Seeds of A.thaliana Col-0 ecotype was obtained from the Arabidopsis Biological Resource Center (https://abrc.osu.edu/). To obtain transgenic lines in which CsPPR was driven by the 35S promoter, the coding sequence of CsPPR was amplified from cucumber and subsequently cloned to the BamHI and PmI I sites of the pCAMBIA1305.4 vector to generate the 35S::CsPPR recombinant plasmid. Agrobacterium tumefaciens strain GV3101 carrying 35S::CsPPR plasmid was used to transform A.thaliana ecotype Columbia (Col-0) plants using the floral dip method (Clough and Bent, 1998). Transgenic seedlings were selected for their resistance to 40 mg/mL hygromycin. The resistant T₁ transgenic seedlings were selected and propagated and T₂ homozygotes selected for further analysis. The expression of CsPPR were compared between the transgenic plants with Col-0 using qRT-PCR.

Seeds of the Col-0 and transgenic plants were harvested at the same time from plants grown under the same conditions. For LT treatment, 50~100 seeds of each genotype were planted on the same 1/2MS medium and exposed to 4°C for 40 days under a 16 h-light/8-h dark cycle. At 4°C for 20 days, the number of seedlings that germinated was recorded, and the germination rate calculated. An obvious emergence of the radicle through the seed coat was used to indicate seed germination. At 4°C for 40 days, the number of surviving individuals was recorded, and the survival rate calculated. For the non-chilling control, treatment and analyses was identical to that described for the chilling treatment, except that the seeds were stored at 22°C for 5 days.

Phenotypic data related to the germination of 151 cucumber accessions of the core germplasm collection were obtained from two environments daily. From these data, four indices, i.e., RGR, RGE, RGI and RRL, were calculated ( Table S3 ). Seedlings were tested in phytotrons in the summer of 2019 and in 2020, i.e., 2019S and 2020S, in Beijing and Shandong province respectively. The phenotypic data showed diverse variation among germplasms ( Figure 1 ). RGR, RGE, RGI and RRL exhibited values of 0-100%, 0-96.59%, 0-75.20% and 0-1.37%, with average values of 29.01%, 21.07%, 16.24%, 0.35%, respectively ( Figures 2A–D ; Table 1 ). There was a clear bimodal distribution of RGR, RGE, RGI and RRL among the accessions, indicating that these traits were quantitative in nature, and were controlled by multiple genes ( Figures 2E–H ). In total, seven cucumber ecotypes were evaluated, including the Indian (n = 3 lines), Japanese (n = 6 lines), American (n = 13 lines), European (n = 17 lines), South China (n = 21 lines), Eurasian (n = 21 lines), and the North China ecotypes (n = 70 lines) ( Table S3 ). Three market-types have higher LT resistant in the two environments, e.g., the North China, the Eurasian, and the South China types ( Figures 2I–L ). The coefficients of variation (CVs) of RGR, RGE, RGI and RRL were 112.05, 130.03, 119.89, and 117.12%, respectively ( Table 1 ). Significant (P < 0.001) positive correlations were observed between each pair of the four traits, i.e., RGR, RGE, RGI, and RRL, when tested in both environments, with the correlation coefficients ranging from 0.56 (RGE_2020S and RRL_2019S) to 0.97 (RGR_2019S and RGE_2020S) ( Figure S1 ). RGR, RGE, RGI in 2019S and 2020S showed high positive correlations, ranging from 0.62-0.97 ( Figure S1 ). These findings indicate that the four traits probably exert a synergistic effect on cucumber seed germination at LT. Meanwhile, data collected from each of the two years were highly correlated and consistent, providing a solid foundation for GWAS analysis.

To evaluate accurately the LT tolerance ability of core germplasms, four principal (prin) components were identified by principal component analysis. Prin1 and prin2 represented 97.59% in 2019S and 98.07% in 2020S of the four indices ( Table S4 ). Weight of prin1 was 95.98%, while prin2 was 4.02% by membership function ( Table S4 ). Finally, a comprehensive evaluation (D-value) of the germination ability of the core germplasm, that encompassed all four germination indices, was calculated. The line with bigger D-value was more resistant for LT. Based on the average D-value per line of two environments, these core germplasms were divided into 5 clusters (I, II, III, IV, V) by the Ward method, when Euclidean distance was 13.5 ( Figure 3A ). Cluster I consisted of highly sensitive accessions, those in Cluster II were sensitive, while Clusters III, IV and V consisted of accessions that were intermediate in tolerance, tolerant, and highly-tolerant respectively ( Table S5 ; Figure 3A ). Each cluster contained different ecotypes ( Figure 3C ). The higher the D-value, the higher the LT resistance of the germplasms ( Figure 3B ). Significant differences were found among various groups ( Figure 3B ). In general, most accessions in the Indian, Japanese and European groups were sensitive to LT, while those in the Chinese group showed higher tolerance ( Figure 3C ). However, the LTG of some germplasms differed greatly among experiments, e.g., the ‘R161’ line grouped in the highly resistant cluster in 2019S, but showed susceptibility to LTs in 2020S ( Figures 2A–D ). Thus, we combined all the phenotypic data from the two experiments to differentiate the susceptible from the resistant lines. Of the 151 tested accessions, we could define 47 as sensitive, i.e., belonging to the highly sensitive group: ‘R13’, ‘R42’, ‘R77’, and 21 as resistant, i.e., belonging to the highly tolerant group: ‘R99’, ‘R152’, ‘R167’.

The various index (D-value, RGR, RGE, RGR and RRL) were used to repeatedly detect loci at the threshold of 5.0 across all 7 chromosomes with well-fitted quantile-quantile (Q-Q) plots in two environments. SNP loci distributed on chromosomes 1, 4, 5 and 6 were detected, and were named gLTG1.1, gLTG1.2, gLTG1.3, gLTG4.1, gLTG5.1, gLTG5.2 and gLTG6.1, respectively ( Figure 4 , Figure S2 , Table 2 ).

Of these seven loci, gLTG5.2 was detected eight times, and two loci (gLTG1.2 and gLTG4.1) were detected four times, using the four indices in each of the two experiments, respectively. These loci were thus considered to have a stable effect on LTG ( Table S6 ). For the rest, four loci (gLTG1.1, gLTG1.3, gLTG5.1, and gLTG6.1) were detected in three times. By comparing the SNP (single nucleotide polymorphism) positions of these 7 loci with QTLs reported in previous studies ( Table 2 ), it was evident that one (gLTG6.1) of the 7 loci in this study colocalized with the known ones ( Figure 5 ), which indicating the locus were reliable via GWAS.

The novel and stable loci on chromosomes 1 (gLTG1.2), 4 (gLTG4.1) and 5 (gLTG5.2) that exhibited significantly association with LT germinability ( Figure 4 S2 ; Table S6 ). So, their candidate regions were further analyzed using LD decay around the peak SNPs, haplotype analysis, functional annotation of the orthologs in Arabidopsis and qRT-PCR. Given that all of them were detected in the GWAS for RGR, therefore RGR is the best choice to further analyze coincidence rate between genotype and loci SNP. Based on RGR in two environments, 21 resistant lines (Group V and IV) and 47 sensitive lines (Group I and II) were used for haplotype analysis ( Table S7 ; Figure S3 ). And ‘R13’, ‘R42’, ‘R77’, ‘R137’ of LT sensitive group (Group I) and ‘R99’, ‘R152’, ‘R167’, ‘R174’ of LT tolerance group (Group V) were used for expression analysis ( Figure 3A ).

For the novel and stable gLTG4.1 locus, the candidate region on Chr.4 was analyzed by pairwise LD correlations leading to a region from 9,703,054 bp to 9,849,931 bp by LD block ( Figure 6A ). Ten annotated genes are located in this region, four of which are associated with different abiotic stresses ( Figure 6B , Table S8 ). CsaV3_4G013370, encodes an enzyme with ATPase and ADPase activity (an apyrase), that are involved in stress-responses in Arabidopsis (Kim et al., 2009); CsaV3_4G013420, encodes a NAC-domain containing protein, which has a role in seed germination (Dekkers et al., 2013) and plant abiotic stress responses (Nakashima et al., 2012); CsaV3_4G013430, encodes a pentatricopeptide repeat-containing (PPR) protein; while CsaV3_4G013480, encodes a RING-type E3 ubiquitin transferase, and takes part in ABA-mediated stress responses (Nakashima et al., 2012). Of these candidate genes, only CsaV3_4G013480 was of interest since five SNPs – three in the exon and two in the intron, were found. Two of the three SNPs are predicted to lead to non-synonymous substitutions ( Figure 6C ), i.e., a Ser →Pro and a Gln→ Glu respectively ( Figure 6C ). We also found that 62% of the 47 sensitive lines carried the Hap1 “CCGCT” sequence, while 86% of the 21 resistant lines had the alternate Hap2 “GGAGA” motif ( Figure 6D ). Interestingly, CsaV3_4G013480 expression was highly induced 1 d after LT treatment in the susceptible lines ‘R13’ and ‘R42’ while in contrast, its expression was lower in the resistant lines ‘R99’ and ‘R152’ ( Figure 6E ). Moreover, after 3 d at 13°C, CsaV3_4G013480 expression in the susceptible lines (‘R13’ ‘R42’ ‘R77’ ‘R137’) was significantly higher than that in the resistant lines (‘R99’ ‘R152’ ‘R167’ ‘R174’) ( Figure 6F ). It is possible that haplotype and the expression of this gene might play a role in cucumber LTG response.

For the novel and stable gLTG5.2 locus, the common peak SNPs were further analyzed and a candidate region from 24,329,867 bp to 24,478,603 bp on Chr.5 was estimated using pairwise LD correlations (r² ≥ 0.6) ( Figure 7A ). In this region, 17 annotated genes were found, and five candidate genes related to abiotic stress identified ( Figure 7B , Table S8 ), they include CsaV3_5G029350 (a serine/threonine-protein kinase), CsaV3_5G029450 (a zinc finger protein), CsaV3_5G029520 (A type I inositol polyphosphate 5-phosphatase), CsaV3_5G029540 (a cysteine proteinase inhibitor), CsaV3_5G029550 (a cytochrome P450). Nevertheless, 84 SNPs significantly associated with LTG in the gene CsaV3_5G029350 ( Figure 7D ). Two SNPs in the third exon caused a base change from A-to-T and from A-to-G, which resulted in amino acid change from Thr→Ser and Ile→Val, respectively ( Figure 7D ). Two SNPs located in the fourth exon, led to a synonymous mutation. All 47 sensitive lines were (Hap1), carried the haplotype “AAGG” while 24% of the resistant lines had the “TGAA” (Hap2) ( Figure 7E ). However, the SNP variation in the resistant lines were not consistent with the phenotype, indicating that cucumber LT germinability could be controlled by multiple genes (such as gLTG1.2, gLTG4.1, and gLTG5.2). The relative expression of CsaV3_5G029350 was higher than that of the susceptible lines ‘R13’, ‘R42’ at 1 d at 13°C ( Figure 7F ). In addition, the expression level of CsaV3_5G029350 in the susceptible lines (‘R99’ ‘R152’ ‘R167’) was significantly higher than in the susceptible resistant lines (‘R42’ ‘R77’ ‘R137’) after 3 d at 13°C ( Figure 7G ).

For the gLTG1.2 locus, the signals were detected that consistently exceeded a significant threshold (−log₁₀ P ≥ 5.0) in 2019_RGR, 2020_RGR, 2020_RGE, 2020_RGI. The haplotype blocks were estimated using Plink tool (Purcell et al., 2007) in the gLTG1.2 locus; this effort further narrowed this region into only one haplotype block harboring significantly associated SNPs (Chr1: 29,185,530—29,524,806 bp) ( Figure 8A ). Fourteen annotated genes were located in this haplotype block ( Figure 8A , Table S8 ). And only one gene (CsaV3_1G044080), encoding pentatricopeptide repeat-containing (PPR) protein (namely CsPPR), was reported to be related to abiotic stress. In Arabidopsis, AtPPR played a central role in abiotic stress during seed imbibition and seed development stages (Liu et al., 2016). AtCsPPR might interact with AtOTP ( Figure S2 ) (Lee et al., 2009). To determine if the candidate gene CsPPR could function in a LTG response in cucumber, the haplotypes of 21 resistant and 47 sensitive cucumber lines were compared ( Figure S3 ). Six variants of the significantly associated SNPs were identified in the 68 lines ( Figure 8B ). For 21 resistant cucumber lines and 47 sensitive cucumber accessions, 21 significantly associated SNPs located at the upstream (< 2 kb) and genic region of CsPRR were classified into 2 haplotypes ( Figure 8B ). Approximate 50% sensitive lines had the haplotype2 (Hap2; “ACTTGT”), while 86% resistant lines carried the Hap1 (“TTCAAC”) ( Figure 8C ). By further analysis of the accessions with different haplotypes, the “ACTTGT” (Hap2) haplotype was only existed in European type, American type, Indian type, South China type and Eurasian type ( Figure 8D ).

Of these significantly associated SNPs, in CDS region, one located at nucleotide 29,278,310 resulted in a C-to-T conversion, but led to no amino acid changes ( Figure 8B ); many cis-elements in CsPPR promoter were related to LT ( Figure S4 ), eg. ABRE (ABA-responsive element), MYB (cold responsive element), ERE (ethylene-responsive element), CGTCA-motif (MeJA-responsive element), TGACG-motif (MeJA-responsive element). Of these, 5 motifs in promoter were altered by 5 SNP. SNP_29,275,643, SNP_29,277,326, SNP_29,277,544 caused change in TATA-box, SNP_29,276,402 bp caused change in TATA-box, and SNP_29276532 caused changes in MYB recognition site, respectively ( Figure 8B ). But the remaining SNPs has no change function of gene. The expression level of this gene in resistant lines (‘R99’ ‘R152’) was significantly higher than that in sensitive lines at 1 d after LT treatment ( Figure 8E ). In a larger panel of 8 randomly selected accessions, a general trend was observed that the expression level of resistant lines (‘R99’ ‘R152’ ‘R167’) were significantly higher than that of susceptible lines (‘R42’ ‘R77’ ‘R137’) at 3 d at 13°C ( Figure 8F ).

To further determine the reliability of the GWAS result, we determined if candidate gene CsPPR, identified at the gLTG1.2 locus, could have a role in seed germination under LT stress. A recombinant plasmid 35S::CsPPR was introduced into A. thaliana by A. tumefaciens-mediated floral dip method ( Figure 9A ). Five homozygous T₂ transgenic lines (OE1-5) were obtained. The transgenic lines (OE-4, OE-5) were randomly selected for further experiments.

The expression level of wild-type (Col-0), the lines ectopically expressing CsPPR, i.e., OE-4 and OE-5, were detected by qRT-PCR. OE-4 and OE-5 had 6-fold and 27-fold higher PPR transcripts compared to those in Col-0 ( Figure 9D ). Between 50-100 seeds of Col-0, OE-4 and OE-5 were germinated at 22°C for 5 days or 4°C for 40 days. Col-0 generally showed germination inhibition, and the germination rate was significantly lower than that of the CsPPR transgenic seeds (OE-4 and OE-5) at 4°C for 20 days ( Figures 9B, G ). In contrast, the germination rate of Col-0, OE-4 and OE-5 were same when exposed to a non-chilling temperature, i.e., 22°C for 5 days ( Figures 9B, F ). Further, the survival rates of the transgenic lines (OE-4 and OE-5) were higher than that of the Col-0 genotypes at 4°C for 40 days ( Figures 9C, E ). Interestingly, CsPPR expression showed a dosage effect for survival rate. For example, the expression of CsPPR on OE-5 had 5-fold higher transcripts than OE-4; Meanwhile, the survival rate of OE-4 was higher than Col-0 by over 63%, and the survival rate of OE-5 was 86%, higher than that of Col-0. These data were consistent between cucumber and Arabidopsis. Viewed collectively, these results support that CsPPR is very likely the causal gene to regulate LTG.