Controlled conditions and Heliaphen phenotyping platform

To study long-term impacts of early cold stress on late traits and oil yield, plantlets were stressed for three weeks in a growth chamber as described previously (Fig. 1C), and then each was transferred outdoors to a 10 L pot filled with fertile soil (PAM2, Proveen) on the high-throughput phenotyping platform Heliaphen (Gosseau et al., 2019) until the end of its life cycle. The chlorophyll, anthocyanin and flavonoid contents and nitrogen balance of leaves 1, 3, 5 and 10 or more were measured three times per week until senescence of the leaf using a leafclip sensor (Dualex®, Force-A).

Field experiments

Two field experiments were performed to study impacts of the sowing date (i.e., early vs. normal) on sunflower development and oil yield. The first field experiment (21TE01-02) was conducted at the INRAE station in Auzeville, France, in 2021. ES and normal sowing were performed on 1 March and 27 May 2021, respectively. Once half of the plants in each plot had developed a true pair of leaves, unmanned aerial vehicle (UAV) flights were used to take images of the experiment three times per week until half of the plants had three true pairs of leaves.
Analysis was based on that of Madec et al. (2017). Plant height was calculated from a point cloud obtained through structure from motion of UAV RGB (i.e., red, blue and green) images. Photogrammetry was performed using Photoscan software (Agisoft LLC, St. Petersburg, Russia). We then defined two regions: the soil buffer (larger than the original plotmap), used to calculate the soil baseline, and the inner buffer, used to calculate height. Then, outlier 3D points that lay outside the 1st to 99.9th percentile interval for height were identified and removed. Soil height was defined as the lowest point in the point cloud (20thpercentile), while plant height was defined as the distance from the soil to the 99th height percentile in the inner buffer, in a UTM-based reference system. To assess the variability in height, 30 cm × 30 cm cells were analyzed in the same way, and the variability was calculated as the coefficient of variation of height among the cells. The canopy was defined as the percentage of green pixels in the plot. Green pixels were detected using excess green index (ExG = 2 × G - R - B) thresholding. The green cover thus decreased as the density of yellow ligules increased.
The second field experiment was performed at Grisolles (43.815356°N, 1.290086°E), France, in 2022. ES and normal sowing were performed on 1 March and 11 April 2022, respectively. Seeds were sown every 12.5 cm, in 5 m rows separated by 60 cm. The number of plantlets with fully developed cotyledons was recorded three times per week until half of the plots developed them. Similarly, the number of plants with leaves 4 cm long was recorded three times per week until half of the plant on each plot developed them. Visual inspection of growth curves identified the exponential stage of growth for vigor. The last four measurements were used to calculate the growth rate of above-ground cover during the exponential stage, while all measurements were used to calculate the growth rate of height.

Transcriptome analysis

The first two pairs of leaves, hypocotyl and entire root system were harvested after three weeks of growth for plants in both the night chilling and control treatments (Fig. 1C). Organs from two plantlets were sampled between 10:00 and 12:00, pooled and immediately flash frozen in liquid nitrogen. Frozen samples were ground using a grinder (ZM200, Retsch, Haan, Germany) with a 0.5 mm stainless steel sieve. RNA was extracted from 100 mg of frozen powder using an RNA kit (NucleoSpin®, Macherey-Nagel, Düren, Germany) following the manufacturer’s instructions. To eliminate genomic DNA contamination, samples were then treated with DNAse (TURBO, Thermo Fisher Scientific, Waltham, MA, USA).
Concentrations and 260/230 and 260/280 ratios of the RNA were determined using a spectrophotometer (NanoDrop, Thermo Fisher Scientific), and their integrity was verified using a bioanalyzer (2100, Agilent, Santa Clara, CA, USA) and an RNA chip (6000 Nano, Agilent). For each organ per treatment, three replicates were selected for sequencing based on their concentrations and integrity (total: 18 samples). An RNA sequencing kit (Universal Plus mRNA-Seq, Tecan Genomics, Redwood City, CA, USA) was used for library preparation following the manufacturer’s instructions (library type: fr-secondstrand). Paired-end 150 bp sequencing was performed by IGA Technology Services (Udine, Italy) using a sequencer (NovaSeq 6000, Illumina, San Diego, CA, USA), and FASTQ files were submitted to the Sequence Read Archive (www.ncbi.nlm.nih.gov/sra) under accession no. SRP477303.
Raw paired-end reads were processed using the nf-core/rnaseq v 3.0 (doi:10.5281/zenodo.4323183) pipeline and subsequently aligned to the sunflower reference genome HanXRQr2.0 (GCA_002127325.2) usingSalmon v 1.4.0. Read counts were retrieved at the gene level. After a quality-control step based on hierarchical clustering, two samples appeared to be misplaced and were excluded from further analysis.
Read counts were converted to counts per million and then filtered by keeping only genes that had more than three counts per million in at least three libraries. Normalization factors were calculated using thecalcNormFactors function of the edgeR package v 3.40.2 of R Statistical Software (v 4.0.4; R Core Team 2021), and then used to normalize the counts using the voom function of the limmapackage v 3.46.0. The resulting expression matrix was used to perform a scaled principal component analysis using the prcomp function of the stats package v 4.2.2.
Differentially expressed genes (DEGs) were identified using thelimma package v 3.46.0. A linear model was fitted for each gene while setting the treatment as a fixed factor and comparing stressed and control plants. A threshold of 0.05 for Benjamini-Hochberg-correctedp -values was used.
Functional analysis of the transcriptome was performed using two approaches. In the first, the three lists of DEGs were analyzed using ClueGO v 2.5.8 (Bindea et al. , 2009). For hypocotyls and roots, all DEGs were used, while for leaves, only the genes with |Fold Change (FC) | > 1.5 were used. A two-sided test was used to highlight enriched ontology terms. Significance was set at Bonferroni-adjusted p < 0.05, the ‘GO fusion’ option was used, the k-score was set at 0.4 and the analysis considered only gene ontology (GO) levels 3-8. Three custom ontology terms were used: (1) a GO ‘biological process’ term, (2) a KEGG term and (3) a ROS-specific term based on the gene lists of Willems et al. (2016) (Table S5). The three terms were developed as detailed in the supplementary materials (Supplementary Methods). The second approach used the pathifier package v 1.36.0 (with min_std = 0.25 and min_exp = 0.10) for the 5000 genes with the highest variance among all samples. Only the KEGG ontology term was used, and leaf samples in the control treatment were chosen as references. A matrix of pathway deregulation scores (PDSs) was then visualized by producing a heatmap using the gplots package v 3.1.3.
ROS wheel clusters, libraries, differential expression and ClueGO data were submitted to a public portal: https://entrepot.recherche.data.gouv.fr/dataset.xhtml?persistentId=doi:10.57745/4HNS1J

Data analysis

All data were analyzed using R with the dplyr (v 1.0.4), purr (v 0.3.4) and ggplot (v 3.3.3) packages to manipulate data and generate graphs. We first calculated the plasticity P of each trait Y as:
\begin{equation} P_{i,j}=\frac{Y_{cold,i,j}}{{\overset{\overline{}}{Y}}_{control,i}}\nonumber \\ \end{equation}
where i represents the experiment, j the individual plant,\(Y_{cold,i,j}\) the value of the trait in the cold treatment and  \({\overset{\overline{}}{Y}}_{control,i}\) the mean of phenotypes in the control treatment in experiment i.
We used the emmeans package (v 1.5.4) to calculate adjusted means of plasticity among all experiments.
For dynamic phenotypes Y at an intermediate developmental stage, the growth rate was calculated using a sigmoid model of Deng et al. (2012):
\begin{equation} Y=\frac{1}{1+e^{(-g\times(t-t_{i}))\ }}\times Y_{\max}\nonumber \\ \end{equation}
where g represents the intrinsic growth rate, t the number of days after sowing, ti the constant for a given variety i and Ymax the maximum value of the trait.

ROS quantification

Eight leaf disks per plant were harvested and placed in 96 well microplates. Disks were stained following the protocol of Kumar et al. (2014). Plates were then scanned using a table scanner (11000XL, Epson). Images were processed using a custom Python script to segment disk interiors and then transformed to grayscale, and peroxide or anion oxide coloration was measured as the mode of pixel intensities.

Lipid composition

Three replicates of leaves, hypocotyls and root systems in both the cold and control treatments were used to analyze lipid composition. Samples were harvested and ground in the same way as for the transcriptome analysis
Lipids were extracted and purified according to Folch et al. (1957). Briefly, the lipids from 20mg of sample were extracted first with chloroform:methanol (1:1, v/v) and 1 volume of H2O then with 2 volume of chloroform. Polar contaminants such as proteins or nucleic acid were removed by adding 2 volume of NaCl (0.9%, w/v) solution. After phase separation, the lower organic phase, which contains lipids, was harvested and the solvent was evaporated.
Transmethylation of fatty acids was performed overnight in 1 mL of methanol:H2SO4 solution (100:2.5, v/v) containing the internal standards C17:0 (57 µg/mL) at 85°C. Aftercooling, 1 mL of hexane:2.5% NaCl (1:1, v/v) was added, and the upper hexane phase containing FAMEs was recovered. The upper hexane phase was harvested, and diluted 20 time before injection in GC-MS.
GC-MS (gas chromatograph – mass spectrometry) was performed using an Agilent 7890B gas chromatograph and coupled MS detector MSD 5977-EI (Agilent). An DB-23 capillary column (60-m, 250-mm, and 0.25-mm film thickness; Agilent) was used with helium carrier gas at 2 mL/min; injection was done in splitless mode; injector and mass spectrometry detector temperatures were set to 250°C; the oven temperature was held at 80°C for 1 min, then programmed with a 65°C/min ramp to 130°C and a 5°C/min ramp to 250°C (4-min hold). The MS instrument us a MS detector MSD 5977-EI (Agilent). Injection (1µl) was done in splitless mode; injector and mass spectrometry detector temperatures were set to 250 °C. The ion source in an electron ionization set at 70eV, the mass range is 40-700 m/z. Quantification of fatty acid methyl esters was based on peak areas, which were derived from total ion current, and using the respective internal standards.
Lipidomic analyses were performed on the Bordeaux Metabolome Facility-MetaboHUB (ANR-11-INBS-0010).