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).