Two cucumber inbred lines, “CG104” (LT-tolerant) and “CG37” (LT-sensitive) were used in this study. “CG37” and “CG104” were crossed to generate an F₁ (“CG104” x “CG37”) and F₁' (“CG37” x “CG104”) populations. The F₁ (“CG104” x “CG37”) was self-pollinated to produce an F₂ population of 189 individuals. A total of 189 F_2:3 families were derived from F₂ populations accordingly. Moreover, 10 core germplasms with high LT-sensitive or LT-resistant phenotype ( Table S1 ) that have been resequenced (Qi et al., 2013) were used for in silico BSA analysis. The F₂ population was used to verify mutation sites in candidate genes. All materials were preserved by the Cucumber Research Group, Institute of Vegetables and Flowers (IVF), Chinese Academy of Agricultural Sciences (CAAS).

LT treatments were carried out twice in the greenhouse and plastic greenhouse of the IVF, CAAS at Changping, Beijing (40°13'N, 116°05'E) in early spring, respectively ( Table S2 ). Three-week-old seedlings of the two parental lines, F₁ and F₂-derived F₃ families were exposed to 15–17°C for 2 weeks, and the LT injury was classified into six grades, based on the degree of yellowing and dryness of the cotyledons and the first true leaf ( Figure 1 ). The index used was as follows: 0: no symptoms on either the cotyledons or the first true leaf; 1: cotyledons were slightly yellow, while the first true leaf was still green; 3: cotyledons were yellow-green with yellowed edges, while the first true leaf was light green; 5: cotyledons were yellowed in large scale, while the first true leaf was yellow-green; 7: cotyledons were completely yellowed, while the first true leaf was yellow-green; 9: cotyledons and the first true leaf dried.

A low-temperature injury index (LTII) was used as an indicator to indicate the LT tolerance of each plant. The formula used for the calculation of LTII refers to Wang et al. (2019): LTII = (0×S₀+1×S₁+3×S₃+5×S₅+7×S₇+9×S₉)/N×9. S₀–S₉ indicates the number of plants corresponding to each grade. N indicates the total number of plants. Three replicates were set for each treatment, and eight plants of each replicate were investigated. The LTII in two experiments was evaluated by two people, respectively.

A modified CTAB (cetyltrimethylammonium ammonium bromide) procedure (Wang et al., 2006) was applied to extract genomic DNA from the second leaf of each two-leaf stage seedlings. The DNA concentration and quality was examined by electrophoresis on a 1% (w/v) agarose gel and NanoDrop ND-1000 spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). A total of 1,288 SSR markers designed based on the genome sequence of cucumber (Ren et al., 2009) were used to screen the polymorphism of the two parents. Then, a linkage map was constructed with selected polymorphic markers. PCR amplification and gel electrophoresis were conducted referring to the method of Song et al. (2016).

Polymorphic primers were used for genotyping F₂ individuals. QTL analysis was performed with the R/QTL software package (http://www.rqtl.org/). Using individual plant data, a whole genome scan was performed to map the QTLs with composite interval mapping (CIM) procedures (Weng et al., 2015). Genome-wide logarithm of the odds (LOD), threshold values (P < 0.05) for declaring the presence of QTLs were determined using 1,000 permutations. For each detected QTL, a 2-LOD support interval was calculated and defined by left and right markers. The QTL naming format is the abbreviation of the character (i.e., low temperature tolerance—LTT), chromosome (Chr.) number and locus number.

To narrow down the region harboring the major QTL, the association between SNPs and LT tolerance in core germplasm was studied using the in silico BSA strategy. The LT tolerance of 87 sequenced core germplasm (CG) seedlings was previously evaluated in our lab (Wang et al., 2019), two-leaf stage seedlings of CG lines were exposed to 12°C for 11 d and 19.3°C for 14 d, respectively. The LTT of each line was evaluated using the same method as that used to phenotype the F_2:3 families described above. Five highly LT-resistant CG lines (“CG29,” “CG56,” “CG61,” “CG90,” “CG104”) and five highly LT-sensitive CG lines (“CG10,” “CG21,” “CG43,” “CG109,” “CG37”) were selected to generate LT-resistant and LT-sensitive bulks ( Table S1 , Figure S2 ). The 10 CG lines have different geographical origin, and were re-sequenced (Qi et al., 2013). The sequence data was retrieved from National Center for Biotechnology Information (NCBI) Short Read Archive (SRA; Leionen et al., 2011) under accession SRA056480 (https://www.ncbi.nlm.nih.gov/sra/?term=SRA056480). The detailed accession number of each CG line was listed in Table S1 . Then, the SNPs located within the qLTT6.2. region (20,591,185 to 21,186,690 bp) were searched ( Table S3 ) and imported into PilotEdit software. SNPs that were identical within each bulk while different between bulks, were considered to be associated with the LT tolerance.

The parental lines at the two-leaf stage were exposed to 14 h light at 18°C and 10 h dark at 10°C, the second true leaf of each seedling was harvested at 0, 10, 24, 48, 96, 144, 216, 288, 384, and 480 h after LT exposure, respectively, and flash frozen in liquid nitrogen. Total RNA was isolated and the first-strand complementary DNA (cDNA) was synthesized. The qPCR primers were designed with Primer Premier 6.0 (http://www.premierbiosoft.com/primerdesign/index.html). The qRT-PCR was performed using SYBR Premix Ex Taq^™ (Tli RNaseH Plus) (Takara #RR420Q, Takara Bio, Inc. China) in Roche Diagnostics with Light Cycler 480 System, and PCR amplification was conducted following the instructions. Actin1 (Csa3G806800) was used as the reference gene for the normalization of gene expression values (Wan et al., 2010). Each experiment was conducted with three biological replicates and three technical replicates. The relative expression data of candidate gene was analyzed using 2^−ΔΔCt method (Kenneth et al., 2001).

The second true leaf of seedlings at the two-leaf stage was harvested for DNA extraction. Multiple pairs of primers were designed to amplify the full length candidate gene, and neighboring fragments have at least 100 bp overlap. Primers designed were listed in Table S4 . PCR amplification was performed in a 25 µl reaction: 2.5 µl of 10X PCR buffer, 1 µl of 25X Mg²⁺, 1 µl of 10X nucleoside triphosphate, 1 µl of 10 µM forward and reverse primer, 0.25 µl of KOD FX polymerase (Toyobo, Osaka, Japan) and 250–500 ng of DNA. The PCR conditions were as follows: 94°C for 5 min, 30 cycles of (94°C for 40 s, 55°C for 50 s, 72°C for 90 s), 72°C for 10 min, and a final incubation at 4°C for use. The PCR products were then sequenced by Sangon Biotech (Beijing, China). The sequence of the candidate genes of the parental lines were analyzed using SeqMan, and the amino acid sequences were analyzed using the MEGA 6.0 software (Tamura et al., 2013). The promoter sequences of candidate genes were analyzed using the PlantCARE website (http://bioinformatics.psb.ugent.be/webtools/plantcare/html/), and the protein sequences of candidate genes were predicted by SMART software (Ponting et al., 1999).

To further understand the relatedness of the Csa6G445210 and Csa6G445230 sequences, phylogenetic trees were generated. The predicted amino acid sequences encoded by Csa6G445210 and Csa6G445230 in cucumber and orthologues in Cucumis melon, Arabidopsis thaliana, Brassica rapa, Solanum lycopersicum, Oryza sativa cv. Japonica, and Zea mays were obtained using the BLASTP tool (Altschul et al., 1997) in the NCBI database. Multiple sequence alignment was implemented with ClustalW (Thompson et al., 1994) and the phylogenetic tree was structured using the MEGA 6.0 software (Tamura et al., 2013) with a neighbor-joining algorithm.

F₂ individuals were used to identify the correlation between the phenotype (via F_2:3 families) and the SNP mutation sites at the candidate genes. Specific SNP markers at mutation sites were designed ( Table S5 ) to amplify DNA from F₂ individuals, and then the amplified products were digested for genotyping. The individuals that produce the same digestion products as “CG104,” “CG37,” and F₁ were defined as LT resistant (a), sensitive (b), and hybrid (ab), respectively. Then, the genotype and LTII value of each F₂ individual were compared.

All tests for significant differences between “CG104” and “CG37” were done using one-way ANOVA in the R environment (R core Team, 2019).

To phenotype the LT tolerance of the two parental lines, and individuals in the F₁ and F_2:3 populations, 3-week-old seedlings were exposed to 15–17°C for 2 weeks in March 2017 and April 2017, respectively. The LTII was used to indicate the LT tolerance of each plant. Results showed that after 2 weeks of LT treatment, cotyledons and the first true leaf of “CG37” showed severe withering and yellowing, and the LTII score of “CG37” in greenhouse and plastic greenhouse were 59.1 and 57.7 respectively. However, “CG104” had no symptoms and the LTII were 14.9 and 8.7 respectively. The LTII of the F₁ hybrids were 26.4 and 21.7, respectively, with the symptoms more inclined toward “CG104.” Frequency distribution of the LTII among F_2:3 families obeyed a normal distribution ( Figure 1 and Figure S1 ) in two experiments (the coefficient of correlation is 0.87). An additional experiment of phenotyping F₁ and F₁” showed that there is no maternal effect ( Table S6 ), which all together indicate that LT tolerance of cucumber at the seedling stage is a quantitative trait, controlled by multiple nuclear genes.

A total of 1,288 cucumber SSR markers were used to screen the polymorphisms between the parental lines; 509 of them showed polymorphisms, with a polymorphism rate of 39.5%. A total of 190 markers with clear electrophoresis strips were used to generate a linkage map ( Figure 2 ). The markers selected were evenly distributed on seven chromosomes. The total length of the genetic map was 990.8 cM, and the average genetic distance between markers was 5.2 cM ( Table S7 ). Based on the “9930” genome (Huang et al., 2009), the order of all markers on the genetic map was consistent with their physical location, therefore the map was qualified to be used for subsequent QTL mapping.

To detect QTL responsible for LT tolerance in cucumber seedlings, the phenotypic data for LT tolerance (LTII) from the two experiments and the genetic map constructed were used for QTL mapping. Details of each QTL detected, including map location, peak location, peak logarithm of odds (LOD) support value, confidence intervals, and percentages of total phenotypic variances explained (R²) were shown in Table 1 . In total, three QTLs including qLTT5.1 on Chr.5, qLTT6.1 and qLTT6.2 on Chr.6 were repeatedly detected (with a threshold of 3, and the confidence intervals of 2-LOD) in two experiments ( Table 1 and Figure 3 ). Furthermore, LOD scores of qLTT6.2 is 18.2 and 16.8, and accounts for 26.8 and 24.1% of the phenotypic variation. Therefore, it was proposed that qLTT6.2 is a major QTL which is responsible for LT tolerance in cucumber seedlings.

Because qLTT6.2 showed the highest contribution to LT tolerance in two experiments, we further conducted in silico BSA to narrow the region harboring this locus. Our previous study included the LT tolerance evaluation of 87 CG lines ( Figure S2 ). Among these, five lines: “CG29,” “CG56,” “CG61,” “CG90,” “CG104” with strong LT resistance phenotypes, and another five lines: “CG10,” “CG21,” “CG43,” “CG109,” “CG37,” with high LT sensitive phenotypes were used to generate LT-resistant and -sensitive bulks to identify the SNPs variation. The detailed information of these 10 lines including accession name, accession number in NCBI, geographical origin, and their LT tolerance performance are listed in Table S1 . There were 214 non-synonymous SNPs that were identical within each bulk, but different between bulks ( Table S3 ), all of which were distributed within the 42-kb region (20,779,616 to 20,821,620). These data indicated that this 42 kb region was associated with LTII variation among these lines. Annotation of the 42 kb genomic region predicted three genes including one for an auxin response factor, one for a CASP-like protein, and an A. thaliana ethylene insensitive 2 orthologue ( Figure 4 and Table 2 ).

To study the expression pattern of the three candidate genes under LT stress, leaf tissues of the two parental lines exposed to LT stress for 0, 10, 24, 48, 96, 144, 216, 288, 384, and 480 h were harvested, respectively, and then qRT-PCR of each gene was performed ( Figure 5 ). The expression level of Csa6G445210 in “CG104” was significantly higher than that of “CG37” at most time-points, and Csa6G445230 showed a similar expression pattern. However, there was no significant expression difference of Csa6G445220 between the two parental lines. In short, Csa6G445220 did not respond to LT stress, however, both Csa6G445210 and Csa6G445230 were up-regulated by LT stress in “CG104.”

To further analyze the sequence of Csa6G445210 and Csa6G445230, full-length DNA and cDNA of each gene in the two parental lines were sequenced and compared ( Figure 6 , Figures S3 and S4 ). The results showed that two base pair substitutions were detected within the coding sequence region (CDS) in Csa6M445210. The first substitution at position 1733 did not change the amino acid, however a polymorphism at position 1807 resulted in a non-synonymous change (Arg→Cys). Csa6M445230 has seven base pair substitutions within the CDS. Two substitutions were synonymous resulting in no amino acid change, while the rest were non-synonymous (Asp257→Ser, Thr820→Ser, Ser2468→Leu, Leu2678→Pro, Ser3260→Cys). The Asp257 resulted in changes in the transmembrane helix ( Table S8 ).

The consistency between the LT tolerance phenotype, and the genotype of non-synonymous SNP sites within these two genes were examined in F₂ population. Primers were designed based on the non-synonymous SNP site of Csa6G445210 and the first base substitution mutation site of Csa6G445230. By comparing the genotype and LTII value of each F₂ individual, it was found that both non-synonymous SNP sites within Csa6G445210 and Csa6G445230 were associated with the LT tolerance phenotypes ( Table S9 ).

The gene annotation suggests that Csa6M445210 encodes an auxin response factor (ARF), which is a transcription factor that controls the expression of auxin responsive genes. Csa6M445230 encodes a transmembrane protein, ethylene-insensitive 2 (EIN2). To further predict and analyze the function of these two predicted amino acids, the sequences analyzed by BLAST against the NCBI database, and phylogenetic analyses were performed with their orthologues in seven other species ( Figure 7 ). The results showed that the homologs of protein encoded by Csa6M445210 and Csa6G445230 are highly conserved in A. thaliana and other plant species, which suggests that they may share similar functions.