Transcriptomic diversity in seedling roots of European flint maize in response to cold

Low temperatures decrease the capacity for biomass production and lead to growth retardation up to irreversible cellular damage in modern maize cultivars. European flint landraces are an untapped genetic resource for genes and alleles conferring cold tolerance which they acquired during their adaptation to the agroecological conditions in Europe. Based on a phenotyping experiment of 276 doubled haploid lines derived from the European flint landrace “Petkuser Ferdinand Rot” diverging for cold tolerance, we selected 21 of these lines for an RNA-seq experiment. The different genotypes showed highly variable transcriptomic responses to cold. We identified 148, 3254 and 563 genes differentially expressed with respect to cold treatment, cold tolerance and growth rate at cold, respectively. Gene ontology (GO) term enrichment demonstrated that the detoxification of reactive oxygen species is associated with cold tolerance, whereas amino acids might play a crucial role as antioxidant precursors and signaling molecules. Doubled haploids representing a European maize flint landrace display genotype-specific transcriptome patterns associated with cold response, cold tolerance and seedling growth rate at cold. Identification of cold regulated genes in European flint germplasm, could be a starting point for introgressing such alleles in modern breeding material for maize improvement.


Background
Maize displays the most widespread geographical distribution of all major crop species [1] with an annual grain harvest of 1135 million tons [2]. In the EU-28 countries, maize is grown second only to wheat by production [3]. Although maize has been adapted to a variety of environmental conditions, traits such as disease, insect resistance and abiotic stress tolerance can be further improved in elite germplasm subjected to a rapidly changing climate [4]. Since the introduction of maize in Europe, geographical separation and natural as well as human selection led to a diversification of landraces. Molecular analyses discovered that traditional flint corn (Zea mays var. indurata) populations of Northern Europe have major contributions from North American flints, which were introduced to Europe during the sixteenth century [5,6]. Adaptations to the North and Central European climate included the development of a shorter growing cycle to avoid cold temperatures during the growth period and as another strategy, higher tolerance to cold temperatures [7]. A high genetic diversity, including favorable alleles to improve elite germplasm are present in European flint populations [4]. In the common dent x flint hybrids, the flint line represents in most instances the cold tolerant parent. However, to date, variation for cold tolerance in elite hybrids is scarce and maize is highly cold sensitive [8].
Due to its tropical origin, the optimum temperature for maize growth ranges from 21 to 27°C [9]. Suboptimal temperatures decrease the capacity for biomass production and lead to growth retardation. Upon exposure to temperatures below 10°C, which often occur at sowing time in Central and Northern Europe, cellular and tissue injuries may cause irreversible damage and may result in plant death [9]. Response to cold stress in maize has been studied broadly and many affected pathways have been identified, revealing the complexity of cold stress response. The immediate reactions to cold involve the decrease of CO 2 assimilation and the down-regulation of photosynthetic electron transport in leaves, inhibiting photosynthesis [9][10][11]. Cell cycle duration and cell proliferation are reduced [12] and cell-wall organization is changed [13]. Furthermore, chilling activates different defense mechanisms. Antioxidant production and activity are altered as a result of increasing levels of reactive oxygen species (ROS) [14]. Changes in gene expression under cold stress include the repression of photosynthesis related genes. An induction is observed in genes related to transcription, phosphorylation, cell-wall organization and expression regulation [13]. Induced regulators are for example many phytohormones [15] and among those in particular salicylic acid (SA) and abscisic acid (ABA) [16,17]. The impact of cold temperatures on root development and function has been less explored. One effect of cold temperatures is the reduction of hydraulic conductance of roots [18], which leads to water deficit of the plants under cold stress. Maize seedlings can acclimatize within 24 h which results in recovery of hydraulic conductance [19]. Further, the lengths of elongation zones and root growth-rates are reduced under cold stress thus affecting root architecture. In particular, the branching angles between primary and lateral roots are reduced upon cold stress [20].
To ensure high yield in temperate climates, a good early seedling vigor during cold temperatures is important. In Germany, maize is typically sown between mid-April and May, where temperatures regularly drop below 10°C [21]. However, early sowing is advised by agricultural consultants in colder regions [22]. This strategy improves the performance throughout the year because the maize plants benefit from a longer vegetation period, improved vegetative growth and earlier ripening and harvest times. By early sowing, plants can avoid summer drought during flowering and ripening [23,24]. Current agronomic strategies to reduce chilling effects in maize involve adaptation of sowing time and soil management such as preparation of a fast warming seedbed or mulching [14]. Breeding for cold tolerance during early development will also be important for no-tilling conservation agriculture, where soil warming is slower [25]. Therefore, inclusion of maize varieties with cold tolerance during early development will be important for environmentally protective agricultural practices.
Maize landraces are a rich source of favorable alleles for broadening the genetic basis of elite germplasm [4].
We hypothesized that European maize landraces display substantial variation for cold tolerance thus carrying beneficial alleles for this trait which might not be present in elite material. In this study, the transcriptomic response of pre-selected doubled haploid (DH) lines derived from the European flint landrace "Petkuser Ferdinand Rot" towards cold treatment was evaluated towards the goal to improve cold tolerance. DH technology is an efficient method to generate homozygous DH inbred lines by chromosome doubling of haploid cells [26].
The aims of this study were to identify genes associated with (i) the general response to cold treatment (ii) cold tolerance and (iii) seedling growth at cold conditions.

Results
Identification of cold tolerant and susceptible maize genotypes from doubled haploid lines derived from the flint landrace Petkuser Ferdinand rot European elite maize germplasm shows only limited variation for the response to cold stress during early seedling development whereas European flint landraces harbor a high genetic diversity for cold tolerance. To untap this underutilized genetic resource, we screened 276 doubled haploid (DH) lines induced from the flint landrace 'Petkuser Ferdinand rot' for the growth response of primary roots to cold.
The primary root of maize was assessed in this study because it is the first organ which emerges after germination. Low soil temperatures at early stages of root development can negatively affect early seedling vigor which is important to ensure high yield in temperate climates. Moreover, the simple structure of the seedling primary root allowed to define the duration of cold stress treatment for the RNA-seq experiments during which no morphological changes of the root were monitored.
In this study we defined cold tolerance as the ratio of growth rate of the primary root at cold versus growth rate of the primary root at control conditions. As a second measure we defined growth rate at cold as absolute growth values under cold conditions in cm/day. These two measures are plotted in Fig. 1a. Cold tolerance and growth rate at cold were correlated significantly (p < 0.001). However, the coefficient of the correlation between the two traits in the DH-population was low with r = 0.22.
Based on the results of this phenotyping experiment, we selected ten cold susceptible (yellow to orange features) and eleven cold tolerant DH-lines (green features) of the 276 genotypes for further analyses based on seed availability (Fig. 1a). The ranking of the genotypes with respect to cold tolerance (Fig. 1b) differed from the ranking of genotypes with respect to the root growth rate under cold conditions (Fig. 1c).  Transcriptome profiling of cold tolerant and susceptible maize DH lines Subsequently, we surveyed how the diversity for cold tolerance is reflected in the transcriptomic landscape of the primary roots in the selected DH lines. The transcriptomes of the eleven cold tolerant and ten cold susceptible lines were investigated after control and cold treatment in four biological replicates per genotype by treatment combination. The RNA-Seq experiments yielded on average~36 million 100 bp paired-end reads per sample (Table S1). The sequencing data has been deposited in the NCBI sequencing read archive (SRA; http://www.ncbi.nlm.nih. gov/sra; BioProject accession number PRJNA556806). Among those, on average 75% of the trimmed highquality reads mapped at unique positions in the gene set of the maize B73 reference genome with 46,272 predicted coding and non-coding gene models (AGPv4 release 36) (Table S1).
We considered a gene active in a genotype if we detected on average ≥ 1 fragment per million reads (FPM) across all eight biological replicate samples of a genotype. The number of active genes ranged from 19,917 in PE0040 to 21,011 in PE0075 ( Figure S1). Overall, 24,448 different genes were active in at least one genotype while 17,204 genes represented the core transcriptome, i.e. genes active in all 21 genotypes.

Kinship relationship among the surveyed panel of maize DH lines
We determined the transcriptomic relationships among the 21 tested maize DH lines under cold and control conditions by a principal component analysis (PCA). In the PCA, the two principal components PC1 and PC2 explained 21% of the total variance (Fig. 2). The samples subjected to control and cold treatment clustered closely together, respectively for each genotype, indicating small overall transcriptomic differences between cold and control treatment. The very cold tolerant genotypes PE0161 and PE0002 clustered closely together and were clearly separated from the other genotypes (Fig. 2). We did not observe separation of the remaining tolerant versus susceptible genotypes. The transcriptomic relationship of the surveyed samples was in all instances mainly determined by the genotype and to a smaller extend by the treatment.

Three types of transcriptomic responses to cold
We identified three types of transcriptomic responses associated with cold stress including treatment, cold tolerance and growth rate at cold (Fig. 3). Treatment, cold tolerance and the interaction between these two were determined with model 1 (see Methods). We identified genes which were differently expressed for treatment and cold tolerance, while for the interaction term no genes were found to be differentially expressed.
To investigate the gene expression changes of the studied genotypes associated with their different growth rates at cold conditions we applied model 3 (see Methods), where we included only samples from the cold treatment.
Finally, model 2 was applied to break down the transcriptomic response associated with the treatment effect (model 1) on the genotype level. To this end, we applied model 2 to each genotype in a separate analysis with the factor treatment alone. Thus, we were able to determine dependence of gene expression on the treatment effect for each genotype and refine the results of model 1.
Differentially regulated genes were computed by pairwise contrasts in the case of two factor levels or significances of continuous factors in the cases of quantitative traits, i.e. cold tolerance and growth rate at cold (FDR < 5% and |log 2 FC| > 1).
First, 148 genes differentially expressed upon cold treatment irrespective of the genotype (factor t j , model 1, see Methods). Second, 3254 genes, which were differently expressed with increasing cold tolerance irrespective of the expression differences between treatments (factor c i , model 1, see Methods) and third 563 genes associated with genotypic growth at cold conditions (factor g i model 3, see Methods). No gene was identified as differentially expressed with respect to the interaction effect ((ct) ij model 1, see Methods). In total, 27 genes were shared between the treatment and cold tolerance effect and 282 genes where differentially expressed with respect to cold tolerance and growth rate at cold. No gene was shared between treatment and growth rate at cold.
Among the 148 genes associated with the treatment effect (Table S2), ten genes were downregulated in ≥15 genotypes and 12 genes were upregulated in ≥15 genotypes (model 2, see Methods, Fig. 4a, Table S3). The genes with the highest number of genotypes with differential expression upon cold treatment were a heat stress transcription factor C-1 (Zm00001d016255) which was upregulated in 20 genotypes (Fig. 4b) and a plant cysteine oxidase 2 (Zm00001d039166), which was downregulated in 19 genotypes (Fig. 4c).
Gene ontology (GO) term enrichment analysis with the 148 treatment-effect genes yielded 59 significantly (p < 0.05) enriched GO terms (Table S4). We identified a connected network of significantly (p < 0.01) enriched GO terms in the set of treatment-associated genes which was related to the response to cold, water deprivation and different organic compounds as well as with lightdependent processes and hormone signaling (Fig. 4d).
The 3254 genes associated with cold tolerance (Table S5) displayed a gradient of gene expression along the susceptible versus tolerant genotypes Fig. 5a-c. Some genes displayed a consistent trend of gene expression change along the susceptible to tolerant genotypes ( Fig. 5a-b). For instance, a gene encoding a bHLH-transcription factor 136 (Zm00001d021019) (Fig. 5a) displayed a decrease in Fig. 3 Venn diagram of the number of differentially expressed genes with FDR < 5% and |log 2 FC| > 1 for the factors treatment (cold vs. control), cold tolerance, and growth rate at cold. Numbers show overlaps of genes between different classes and uniquely differentially expressed genes for each class expression with increasing cold tolerance whereas the increase of expression of a GRAS-transcription factor 82 gene (Zm00001d048682) correlated in general with increasing cold tolerance of genotypes (Fig. 5b). Other genes of that set displayed more pronounced differences between a subset of the susceptible and tolerant Genes differentially expressed upon cold treatment. a Genes preferentially expressed in either cold or control treated plants in ≥15 genotypes (FDR < 5% and |log 2 FC| > 1). b The heat stress transcription factor C-1 (Zm00001d016255) was significantly upregulated upon cold treatment in 20 of 21 genotypes. c The cysteine oxidase 1 (Zm00001d039166) was downregulated at cold treatment in 19 of 21 genotypes. d Network graph of enriched (p < 0.01) GO terms in the set of differentially expressed genes by treatment effect. Node size represents frequency of the GO term in the underlying database (the smaller the more specific). Node fill color represents the log 10 of the enrichment p-value (the darker, the more enriched). Highly similar GO terms are linked by edges (grey lines), where edge width indicates the degree of similarity  Table S6). Values of each sample represent mean across four replicates/total mean across all samples. A gene with unknown function (Zm00001d031037) was significantly differently expressed between cold and control in DH-line PE0161. Read counts were adjusted by dividing the sample mean by the total mean across all samples. e Network graph of enriched (p < 0.01) GO terms in the set of differentially expressed genes by cold tolerance. For node size, fill color and edges refer to scale and description in  genotypes. For instance, the Aux/IAA-transcription factor 14 gene (Zm00001d049141) displayed relatively low expression in the three most tolerant genotypes (Fig. 5c). As a subset of the 3254 genes differentially expressed with increasing cold tolerance, we identified ten genes, which were exclusively expressed (FPM ≥1) in the single most tolerant (PE0161, seven genes) or the two most tolerant (PE0161, PE0328, two genes) or the three most tolerant (PE0161, PE0328, PE0002, one gene) genotypes (Fig. 5d, Table S6). Among the seven genes uniquely expressed in genotype PE0161 the gene of unknown function, Zm00001d031037 displayed a particularly high expression of > 1000 reads at control conditions and was significantly downregulated under cold conditions (p < 0.05; log 2 FC = − 0.82; average reads = 44; Fig. 5d). This was the only of the ten genes exclusively expressed in the most tolerant genotypes which showed differential expression between cold and control treatment in a genotype.
In the set of 3254 differentially expressed genes for the cold tolerance effect, we identified 64 enriched GO terms (Table S7). We detected connected networks of GO terms associated with amino acid transport, catabolic processes and cellular oxidant detoxification (Fig. 5e).
The ranking of the genotypes with respect to growth rate at cold conditions was different from the order with respect to cold tolerance, where the y-axis shows growth rate at cold values, cold tolerance is represented by colors with green as cold tolerance and red as cold susceptible ( Fig. 1b and c). We detected 563 genes with expression patterns associated with root growth under cold conditions (Table S8). Among those, a gene encoding a EXORDIUM protein (Zm00001d018106) displayed in general low expression in faster growing genotypes at cold conditions (Fig. 6a). In contrast, a gene encoding a pyrophosphate-energized vacuolar membrane proton pump 1 (Zm00001d037492) displayed higher expression levels in the genotypes with faster root growth at cold conditions (Fig. 6b) and was in general highly expressed taking into account expression across all genotypes (average read count per sample: 12,424).
In the set of 563 genes associated with growth rate at cold conditions, we detected 29 enriched GO terms (Table S9). Similar to the GO networks of the genes associated with cold tolerance, two connected networks were associated with amino acid transport and catabolic processes (Fig. 6c).

Discussion
In this study, we investigated the genetic diversity for cold tolerance of the European maize flint landrace "Petkuser Ferdinand rot" during early root development. The rationale to survey cold tolerance was that improvement of this trait will allow earlier planting for increased biomass production. At the same time, as a consequence of earlier planting the risk of potential summer drought will be mitigated by an earlier onset of flowering.
Maize landraces are a rich source of favorable alleles for broadening the genetic base of elite germplasm [4]. We hypothesized that European maize landraces display substantial variation for cold tolerance thus carrying beneficial alleles for this trait which might not be present in elite material. This notion is supported by a study which demonstrated that in a panel of 35 European maize flint landraces including "Petkuser Ferdinand rot" used in this study, one landrace covers at least 75% of the genomic variation of all landraces present in this panel [27]. Phenotypic variation in the DH population derived from the European maize flint landrace "Petkuser Ferdinand rot" used in this study was demonstrated by > 4-fold changes in cold tolerance of root growth (Fig. 1). Similarly, extensive diversity for cold tolerance during early development has been observed in other maize landraces for aboveground traits [7,[28][29][30].
Transcriptomic diversity with respect to cold tolerance An RNA-seq experiment resulted in 19,917 to 21,011 active genes in the 21 DH genotypes derived from the European landrace "Petkuser Ferdinand rot" (Fig. S1). In a study surveying a diverse set of seven US inbred lines at a similar stage of root development between 24,923-25,149 genes were active [31]. The lower mapping rate in this study is likely the consequence of mapping the transcriptomes of these Northern flint genotypes on the high quality reference genome sequence of the dent inbred line B73 [32].
In the principal component analysis (Fig. 2), samples from the same genotype subjected to cold and control treatment clustered closely together while cold tolerant and cold susceptible genotypes were not separated clearly, except for the two most cold tolerant lines (PE0002 and PE0161). This indicates that the transcriptomic variation observed in our study mainly reflects general genotypic differences and only to a minor degree differences in cold tolerance or cold treatment effects. Substantial genetic variability between genotypes in contrast to treatment effects, as observed in our experiments, was also described in previous studies on the diversity of European flint landraces accessed through doubled haploids [4,8,27]. Similar results were obtained in a set of diverse European flint and dent inbred lines subjected to heat stress where the transcriptomic responses between heat tolerant and susceptible genotypes were highly divergent [33]. Likewise, in a transcriptome study with the US inbred lines B73 and Mo17 and their reciprocal F 1 -hybrids genotypic differences were much more pronounced than differences between water deficit stress and control conditions (Marcon et al., 2017).
The substantial transcriptomic variation between genotypes in this study provided the opportunity to evaluate the response to cold stress in a genetically diverse germplasm. In contrast, previous transcriptome studies on abiotic stress tolerance such as waterlogging, cold stress, drought stress typically focused on only one or a few maize inbred lines or hybrids under different conditions [15,[34][35][36][37][38][39][40].

Transcriptomic response upon cold treatment
We investigated the effect of mild cold stress on the transcriptomes of 21 DH-lines. Cold responsive genes were defined as genes differentially expressed between control (18°C) and cold (12°C) night temperatures. We identified 22 cold responsive genes which were differentially regulated in ≥15 of the 21 surveyed genotypes (Fig. 4a). Among these, the heat stress transcription factor C-1 (Zm00001d039166) was significantly upregulated upon cold treatment in 20 of 21 genotypes (Fig. 4b). In previous Fig. 6 a and b Examples of gene expression patterns associated with root growth under cold conditions read counts were normalized with plotCounts(). Genotypes are sorted according to their growth under cold conditions (see Fig. 1c). a The EXORDIUM gene (Zm00001d018106) displayed in general lower expression levels in genotypes with faster root growth at cold conditions. b The pyrophosphate-energized vacuolar membrane proton pump 1 gene (Zm00001d037492) displayed higher expression levels in genotypes with faster root growth at cold conditions. c Network graph of enriched (p < 0.01) GO terms in the set of differentially expressed genes expression patterns associated with root growth under cold conditions. For node size, fill color and edges refer to scale and description in Fig. 4d studies this gene was differentially expressed at severe cold treatment and heat stress [41] as well as at a variety of water deficit conditions [35,42]. Hence, this transcription factor responds to several abiotic stress types and might regulate a plethora of downstream target genes. Similarly, the cysteine oxidase 1 (Zm00001d039166), which was downregulated at cold treatment in 19 of 21 genotypes in this study, is also regulated by several other abiotic stress factors such as submergence [37] water deficit [42] and heat stress [41].
GO term enrichment of the 148 genes differentially expressed upon cold treatment, allowed for the identification of molecular processes associated with this type of abiotic stress (Table S4). In total, 20 of 59 enriched GO terms were associated with responses to different abiotic stresses including the responses to light, heat and phytohormones (Table S4). Different abiotic stresses often result in similar responses in plants. They are frequently based on oxidative stress associated with reactive oxygen species such as H 2 O 2 , which serve as inducers of tolerances to abiotic and biotic stresses [43]. GO terms associated with responses to several abiotic stresses were also enriched in an experiment with diverse European maize germplasm response to heat stress [33].
The most enriched GO terms in our dataset associated with response to abiotic stresses were "response to hypoxia", "reponse to karrikin" and "response to cold" (Table S4). Using REVIGO we were able to identify networks of enriched GO terms with high similarity between each other (Fig. 4d). This connected network grouped around the term "response to cold" as the term with highest relevance. The GO terms associated with abiotic responses, which were connected with the "response to cold" were the responses to water deprivation (i.e. drought) and to karrikin. In previous studies, crosstalk between the signaling pathways of drought and cold stress were suggested [44] because many genes are induced by both abiotic stresses. Similarly, an association of cold response with the response to karrikin was observed in Camellia sinensis [45]. Plant-derived smoke, which is the main source of natural karrikin, increases seed germination and the length and fresh weight of maize seedlings, but also regulates reactive oxygen species and their scavenging system [46].

Transcriptomic responses associated with cold tolerance
To investigate the relative primary root growth of the surveyed maize genotypes under control versus cold conditions we determined cold tolerance as the ratio of root growth at cold versus control conditions (Fig. 1b). In total, 3254 genes were associated with this type of cold tolerance. Among those, several genes encoded transcription factors such as bHLH-transcription factor 136 (Zm00001d021019; Fig. 5a), GRAS-transcription factor 82 (Zm00001d048682) (Fig. 5b) and Aux/IAA transcription factor 14 (Zm00001d049141; Fig. 5c). This observation is in line with the notion that a multitude of transcription factors is involved in the regulation of gene expression upon abiotic stress such as cold [47]. Genotype specific expression of Aux/IAA transcription factor 14 (Zm00001d049141) upon cold and heat stress has also been observed between the genotypes B73 and Mo17 [41]. GRAS-transcription factors like Zm00001d048682 (Fig. 5b) have been demonstrated to be involved in the regulation of root and shoot development and in the improvement of plant resistance to abiotic stresses [48]. Genotype-specific upregulation of this gene upon cold stress has been demonstrated for several genotypes [41], although in the present study we did not observe differential regulation upon cold stress, but rather higher expression in cold tolerant genotypes (Fig. 5b).
The subset of 10 cold tolerance associated genes which were uniquely expressed in the one, two or three most tolerant DH-lines (Fig. 5d), might be of special interest for functional characterization, as they might be associated with increased cold tolerance. Remarkably, six out of ten genes have no associated function yet. Our results suggest that these ten genes are highly variable between genotypes. For the gene Zm00001d031037 (Fig. 5d), this is supported by expression analyses of maize leaves of eight diverse inbred lines, where it was expressed in four lines but not in the other four lines [49]. In the single DH-line, where Zm00001d031037 was expressed in our experiment, it was significantly downregulated at cold conditions. Concordantly, Zm00001d031037 was shown to be upregulated at heat in inbred line Mo17 and downregulated at cold in Mo17 and Oh43, but not differentially expressed at either condition in B73 [41].
From the GO terms enriched in the set of genes associated with cold tolerance (Fig. 5e), we identified three networks of connected GO terms with high similarity among each other. The top enriched GO term was "cellular oxidant detoxification" (GO:0098869) (Table S7) which is associated with the reduction of the toxicity of superoxide radicals or hydrogen peroxide. These molecules are associated with oxidative stress occurring when plants are subjected to stresses such as drought, heat or cold.

Transcriptomic patterns associated with growth rate at cold
Another measure of cold performance is the growth rate in cm/day of seedling roots at cold conditions (Fig. 1c). This parameter is of agricultural relevance because it determines the performance of a genotype during early stages of growth under cold spring conditions.
Among the genes related to the trait growth rate at cold, the EXORDIUM protein coding gene (Zm00001d018106) showed consistently low expression in genotypes with high growth rates at cold conditions (Fig. 6a). In previous studies it was demonstrated that this gene is strongly induced by multiple biotic and abiotic stress types including cold and heat [41], submergence [37] and fungal infections [50]. EXORDIUM putatively promotes shoot and root growth in plants [51]. This protein might therefore play a role in the maintenance of growth under adverse environmental conditions in stress tolerant plants to avoid fatal tissue damage.
Among the GO terms enriched in the set of genes associated with growth rate at cold (Fig. 6c) but also with cold tolerance (Fig. 5e) was amino acid transport. Some amino acids can act as precursors for antioxidants or osmolytes or act themselves as osmolytes to prevent the deleterious effect of cold stress. Antioxidants such as glutathione are involved in the response pathways to cold and also in the tolerance to cold as observed in Pinus halapensis seedlings [52]. Glutathione was enriched in cold tolerant seedlings in comparison with cold sensitive seedlings [52]. Furthermore amino acids could act themselves as osmolytes to prevent damages from cold stress [52].

Association between cold tolerance and growth rate at cold
The order of the genotypes with respect to growth rate at cold conditions (Fig. 1c) was different from the order with respect to cold tolerance (Fig. 1b). However, the two parameters of the 276 lines from the 'Petkuser Ferdinand rot' DH population were correlated significantly (p < 0.001) with a low coefficient of correlation of r = 0.22. This implies that genotypes which perform better at cold conditions are also in general more cold tolerant. However, the low coefficient of correlation between the two traits in the DH-population also suggests that the genetic control of these traits is partly independent from each other. This allows combination of both traits during breeding.

Conclusion
The European flint germplasm surveyed in this study showed high variation in their phenotypic and transcriptomic behavior associated with response to cold treatment, cold tolerance and growth rate at cold. Genetic analyses of mutants defective in the cold related genes identified in this study will help to better understand the molecular principles underlying this important trait during early seedling development. This might pave the way for introgressing beneficial alleles of these genes into modern maize by recurrent selection or genome editing thus improving cold tolerance.

Selection of DH lines by phenotypic assessment
We selected 21 doubled haploid (DH) lines from 276 genotypes of a DH population of the flint landrace 'Petkuser Ferdinand rot' [27,53] based on seed availability. This population was generated by the in vivo haploid induction method [54]. The landrace 'Petkuser Ferdinand rot', originating from Northeastern Germany, has been obtained from the federal ex situ Gene Bank for Agricultural and Horticultural Crop Species at the Leibniz Institute of Plant Genetics and Crop Plant Research, Gatersleben, Germany.
Maize roots were imaged after growing five seedlings per DH line under control conditions (26°C/18°C for 16 h light/8 h dark) for 4 days and subsequently for 2 days either under control or under cold (16°C /12°C for 16 h light/8 h dark) treatment to simulate diverging soil temperatures during early establishment of the seedlings. Seedlings were grown in germination paper rolls (Anchor Paper Co, Saint Paul, USA) as previously described [31,55]. Germination paper rolls were randomly arranged in 10 L buckets filled with ca. 3 L distilled water.
We used GiA Roots [56] to determine root network perimeter and calculated root growth-rates by computing a linear regression of root network perimeter from the images before and after treatment. We defined the cold tolerance of each genotype as Cold tolerance = growth rate at cold/growth rate at control.
We selected eleven cold tolerant DH-lines and ten cold susceptible DH-lines based on cold tolerance values, germination rate and availability of seeds for our transcriptome experiment.

Plant growth and RNA isolation
For each biological replicate, 10 seedlings were grown as described above, for 4 days under control condition (26°C/18°C for 16 h light/8 h dark). In this second experiment, the pre-selected DH-lines were only subjected to a treatment (12°C cold and 18°C control) period of 8 h, instead of 2 days. The mild cold treatment of 12°C was chosen to prevent strong damage to root tissue observed below 10°C [9], but to still induce physiological reactions typically occurring at cold stress. After treatment, root tips of ten primary roots comprising the meristematic and elongation zone (tip of the root excluding root hair zone, with maximal length of 1 cm) were excised with a razor blade, pooled, immediately frozen in liquid nitrogen and stored at − 80°C until RNA extraction.
We prepared 168 RNA samples representing 21 genotypes, two treatments (cold and control) and four biological replicates per genotype by treatment combination. Total RNA was isolated using the RNeasy Plant Mini Kit (Qiagen, Hilden, Germany) including on-column DNA removal. RNA quantity and quality were assessed using the 2100 Bioanalyzer and Agilent RNA 6000 Nano Chip (Agilent Technologies, Böblingen, Germany). For all samples, RIN (RNA integrity number, [57]) values ≥8 indicated their high quality and integrity.

RNA-sequencing
cDNA libraries for RNA-seq were constructed according to the TruSeq RNA sample preparation protocol (Illumina, San Diego CA, USA). Library indexing, cluster preparation, and paired-end sequencing were performed according to the manufacturer's instructions (Illumina). Sequencing of 100 bp paired-end reads was performed on an Illumina HiSeq 4000 platform.

RNA-seq data processing
RNA-seq raw reads were quality trimmed using Trimmomatic (version 0.36, [58] The trimmed reads were aligned to the B73 reference genome sequence AGPv.4 (https://www.maizegdb.org/ assembly) using Hisat2 (version 2.1.0, [59]) according to the manual with the parameters for paired-end reads '--min-intronlen 20 --max-intronlen 60000 --knownsplicesites-infile prepared-splice-sites-file'. The file containing the splice sites information was prepared using the Python script hisat2_extract_splice_sites.py included in the Hisat2 binary distribution and an index of the reference sequence assembly was created using hisat2build.
The software samtools (version 1.3.1, [60]) was applied to transform file format .sam to .bam, to order the aligned reads according to their position and divide paired and unpaired read alignments. Reads were counted when they aligned unambiguously with an alignment quality score of > 10 to an exon of a gene as specified in the AGPv4 annotation file (Ensembl release 36, 06/27/17, ftp://ftp.ensemblgenomes.org/pub/plants/ release-36/gff3/zea_mays). Alignment of sequences to the reference genome was performed using HTSeq (version 0.10.0, [61]) with the parameters '-r pos -t exon -i gene_id -s no --secondary-alignments ignore --supplementary-alignments ignore'.
Count tables were produced separately for paired and unpaired read alignments. Count tables were combined by sample adding counts per feature (gene) creating a single count table with 168 samples and 46,272 genes. We removed genes with a sum of less than ten reads across all samples from further analyses.

Principal component analysis
A principal component analysis (PCA) was performed on the expression data using the normalization procedure rlog() implemented in the R package DESeq2 (version 1.18.1, [62]). The plotPCA() function was used for plotting of the data.

Expression analysis
We normalized expression values with library size by calculating fragments per million (FPM) reads using the fpm() function of DESeq2. We declared genes with a mean FPM across samples of a genotype ≥1 as expressed.

Differential gene expression
Expression levels of genes were estimated by the variancemean dependence in the count table based on a generalized linear model using the negative binomial distribution within the R package DESeq2 [62] calculating log 2 fold change (log 2 FC) values between groups of samples specified in the following paragraph. Significance values for log 2 FC values were calculated as Wald test p-value and were adjusted by the Benjamini-Hochberg procedure to obtain false discovery rates (FDR) [63]. To illustrate expression levels of selected genes, raw counts were normalized with DESeq's built-in function plotCounts().
We applied several models to identify differentially expressed genes. First, we used a model with three factors including data of all samples, where Y ij was the expression of respective gene of DHline i at treatment level j. μ was the general mean, c i was defined as the cold tolerance of DH-line i as define above. t j was the treatment level j (cold and control) and (ct) ij was the interaction of cold tolerance and the treatment level. ε ij was the residual error term. The second model was used with samples of each genotype separately to identify genes which are differentially expressed between cold and control conditions in the respective genotype.
where Y j was the expression of respective gene at treatment level j, μ was the general mean, t j the treatment level j (cold and control) and ε j was the residual error term.
In the third model, we input all samples at cold conditions excluding samples at control conditions.
where Y i was the expression of respective gene of DHline i, μ was the general mean, g i the growth rate of DHline i at cold conditions and ε j was the residual error term. With this model, we explained gene expression with relative growth rate of a genotype at cold conditions.

Gene ontology term enrichment
We tested significant (p < 0.05) enrichment of gene ontology (GO) terms in the gene sets assembled using the differential gene expression analysis with FDR < 5% and log 2 FC of > 1 or < − 1, with model (1) with factor treatment and factor cold tolerance as well as with model (3) with factor growth rate at cold [64,65]. GOterms were accessed from the GAMER annotation set for the reference genome sequence AGPv4, a set with high gene coverage for the maize inbred line B73 [66]. We filtered the GAMER annotation for genes expressed in at least one genotype with at least one fragment per million reads on average across all 168 samples (24,448 out of 46,272 annotated genes). The R-package GO.db [67] was used to obtain the biological descriptions for each GO term. The R-package topGO [68] was used to perform weighted Fisher's exact test on the GO terms present in the respective gene set. We pruned the GO hierarchy from terms with less than five annotated genes and used only GO terms associated to biological processes. We did not consider expression values, log 2 FC or FDR values from the differential expression analysis to ensure consistency of the method throughout the study [68]. For each gene set, significantly (p < 0.01) enriched GO terms were then evaluated using REVIGO (http:// revigo.irb.hr/; [69]). Functional redundancies were reduced and the results were visualized based on the semantic similarity, relying on the 'most informative common ancestor' approach. The allowed similarity was set to 0.7 (medium). Zea mays was selected as database, GO term size and all other parameters were left at their default values. Using REVIGO we were able to assess associations and connections between different GO terms enriched in the gene sets and visualize them with network graphs. We visualized the networks using Cytoscape [70], where p-values derived from topGO's enrichment procedure, frequency of the respective term in the database and similarities of the enriched GO terms were depicted as colors and size of connected bubbles (nodes). The placement of the nodes was adjusted to fit the respective figures and does not provide additional information.