Chromatin-Level Differences Elucidate Potential Determinants of Contrasting Levels of Cold Sensitivity in Maize Lines

Maize is a subtropical, cold-sensitive species. However, some varieties of this species have the potential to withstand long-term low temperatures, even at the seedling stage. The molecular basis of this phenomenon has not been determined. In a chromatin-level study, we compared the cold-stress reaction of seedlings of two maize inbred lines showing contrasting levels of cold sensitivity. The cold-tolerant line was selected based on field data and previous physiological and transcriptomic level studies. The first condition of gene expression—chromatin accessibility—was assessed by formaldehyde-aided isolation of regulatory elements method and DNA sequencing. Potentially expressed genes and cis-regulatory sequences open for interaction with transcription factors have been defined. The results of this study suggest that during cold stress, the tolerant maize line shifted resources from growth to defense. This shift was shown by potential hormone-level events—degradation of growth-promoting gibberellins and synthesis of jasmonic and abscisic acids. This finding is congruent with the xeromorphic morphology of seedlings of the cold-tolerant line and their ability to regrow when stress ceases. It is a common reaction of cold-tolerant maize lines. Moreover, in the cold-tolerant line, several genes from the low-temperature signaling pathways were potentially expressed. Additionally, numerous stress-response AP2/EREBP-binding cis-motifs were accessible in the cold-tolerant line. Differently in the cold-sensitive B73 line, MADS-binding cis-motifs were the most abundant. Development of the photosynthetic apparatus is crucial for the survival of maize seedlings at low temperature. Our results suggest efficient photosynthesis in seedlings of the cold-tolerant line, as was described earlier in physiological-level analyses.


Introduction
During seedling establishment, hardy seeds transform into autotrophic plants. This process requires the production of new leaves and the development of the photosynthetic apparatus, as well as the ability to withstand unfavorable conditions. The earliest stages of plant life are critical from the agricultural perspective, as seedling performance is a foundation of successful growth and eventual yield. Therefore, it is necessary to determine both the physiological and molecular bases of the earliest seedling growth. This research is crucial in the case of warm climate species cultivated in the temperate climate. One such species is maize (Zea mays L.), which attracts exceptional interest due to its efficient C4 photosynthesis and notably high productivity under optimal conditions (Osborne and Sack 2012). Maize seedlings at the earliest stages of development are highly cold-sensitive Frascaroli and Revilla 2018). There is considerable variation in cold sensitivity among maize inbred lines (Verheul et al. 1996;Strigens et al. 2013;Sobkowiak et al. 2014Sobkowiak et al. , 2016Grzybowski et al. 2019). However, the molecular-level basis of this trait has not been fully elucidated. Maize is a highly variable species at the genome level, and inbred lines differ in gene content (Swanson-Wagner et al. 2010 and references therein).
Molecular level studies on maize showed that cold-stress leads to profound transcriptome changes. Notably, induction of gene expression dominated over repression in cold-stressed seedlings of tolerant line S68911 relative to control (Sobkowiak et al. 2016). One of the common physiological level changes in cold-stressed maize is inhibition of photosynthesis . Interestingly, seedlings of some maize lines showed induction of expression of chloroplast genes in cold conditions Zhang et al. 2009;Sobkowiak et al. 2016) what could maintain relatively high level of photosynthesis (Sobkowiak et al. 2016;Grzybowski et al. 2019). Protein products of many genes highly expressed in cold in tolerant maize lines have predicted localization in cell wall and membranes, including plasmalemma (Sobkowiak et al. 2014(Sobkowiak et al. , 2016. Plasma membrane is a site of cold reception  and is linked functionally with other highly expressed gene categories in tolerant lines, namely transport and signal transduction Zhang et al. 2009;Sobkowiak et al. 2014Sobkowiak et al. , 2016. Cold stress had strong effect on synthesis and signaling of many hormones, as microarray study has shown (Sobkowiak et al. 2014). However, in most cases, both positive and negative hormonal regulation have been observed in maize seedlings. The exception were genes responsible for ethylene biosynthesis and auxin response factors (ARFs), induced in both two maize lines of different cold sensitivity (Sobkowiak et al. 2014). Notably, the majority of data on the molecular level reaction of maize to cold was obtained in studies on fully developed autotrophic third leaf of seedlings (V3 stage, when the collar of the third leaf is visible, Ritchie et al. (1986)). However, during growth in the temperate climate younger maize seedlings can also experience low temperature, so these stages are equally important in research on maize cold sensitivity (Grzybowski et al. 2019).
Gene expression is regulated at several levels. The first is physical chromatin accessibility. Subsequent levels comprise transcription, mRNA translation, and posttranslational modifications of proteins. Among those levels, the regulation of chromatin accessibility is the least thoroughly characterized. The compact form of chromatin limits DNA accessibility for regulatory proteins and, consequently, the genome available for transcription. The functional genome could amount to only 1% in maize (Rodgers-Melnick et al. 2016). This property is, in part, an effect of a large proportion of noncoding DNA (Schnable et al. 2009), but the functional genome can be strongly influenced by multiple cues, such as developmental transitions (Candaele et al. 2014;Perduns et al. 2015) and stress. It was shown that heat and cold stress have caused histone acetylation and, consequently, chromatin decondensation in maize l e a v e s ( H u e t a l . 2 0 1 2 ; Wang et al. 2 0 1 5 ) . Additionally, water stress has induced chromatin-level changes in various plant species (Han and Wagner 2014). Importantly, chromatin-level regulation has been shown for the crucial element of the drought response in maize, Z. mays dehydration responsive element binding 2A (ZmDREB2A, Zhao et al. 2014). The same finding was observed for Oryza sativa dehydration responsive element binding 1B, which is important in the response of rice to cold (OsDREB1B, Roy et al. 2014). Nucleosome-depleted regions (NDRs) just upstream genes are correlated with their expression level in many eukaryotic organisms, including maize (Hogan et al. 2006;Vera et al. 2014;Rodgers-Melnick et al. 2016;Oka et al. 2017). Determination of NDRs not only provides information on potentially expressed genes but also, more importantly, uncovers the key DNA sequences for gene expression regulation, enabling the construction of gene regulatory networks.
The aim of this study was to determine cold-induced changes in the nucleosome occupancy level of DNA in maize seedlings, and define potentially expressed genes important in cold tolerance. The performance of two maize inbred lines with contrasting cold sensitivity was compared. The coldtolerant S68911 line, which is adapted to the temperate climate (Sobkowiak et al. 2016;Grzybowski et al. 2019), and the reference B73 line, which is cold-sensitive (Joets et al. 2018;Grzybowski et al. 2019), were used. The S68911 line shows high early vigor in cold spring conditions in the field (Grzybowski et al. 2019). When grown from kernels in growth chamber in prolonged cold conditions, S68911 seedlings develop normally, although they exhibit xeromorphic morphology, as has been observed in Grzybowski et al. (2019) and in this study. Xeromorphy has also been shown for other cold-tolerant maize lines grown in cold conditions (Verheul et al. 1996;Sowiński et al. 2003;Strigens et al. 2013;Grzybowski et al. 2019). This growth form has been postulated to prevent water loss and secondary water stress (Strigens et al. 2013). Under cold conditions, B73 plants are retarded, develop only two faint-green leaves, cease to grow, and die (Riva-Roveda et al. 2016;Grzybowski et al. 2019).
The study focused on early maize development; therefore, the first leaves of the seedlings were sampled at three growth stages. Development of the photosynthetic apparatus at the seedling stage determines the survival of the plant in the cold (Grzybowski et al. 2019). The regions of open chromatin (NDRs) were determined using the formaldehyde-aided isolation of regulatory elements (Omidbakhshfard et al. 2014) and sequencing (FAIRE-seq) method.

Growth Conditions and Plant Material
Two maize inbred lines with contrasting cold sensitivity were used. We previously examined the cold-tolerant S68911 line at physiological and transcriptomic levels (Sobkowiak et al. 2016;Grzybowski et al. 2019). The B73 line was chosen as a reference because it is cold-sensitive (Joets et al. 2018;Grzybowski et al. 2019) and its genome has been de novo sequenced (Schnable et al. 2009). Kernels of the S68911 line were provided by Prof. Józef Adamczyk (Plant Breeding Smolice Co. Ltd., Poland). Kernels of the B73 line were obtained from Maize Genetics Cooperation Stock Center (USA) and were propagated in a greenhouse at the Faculty of Biology, University of Warsaw (Poland). In this work, we investigated the regulation-level basis of the outstanding cold tolerance of the S68911 maize line.
Plants were grown from kernels in a growth chamber under a 14 h/10 h photoperiod with light irradiance of 500 μmol quanta m −2 s −1 during the light phase and 70% air relative humidity. The temperature was set according to treatment (control, cold, day/night temperature 24°C/20°C and 14°C/10°C, respectively) from the start of the experiment. For plant growth, special installation was used, consisting of impermeable polyvinyl chloride tubes packed with soil and placed in a container filled with water (for full description, see Grzybowski et al. 2019). The water was cooled by a thermostat, which enabled the soil temperature to be controlled. The thermostat was set so that the soil temperature was the same as the air temperature. There were 10-12 plants per tube, and the plants were watered as needed with tap water. The position of the tubes in the container was randomized several times during the experiment to cancel the effects of the potential temperature gradient in the growth chamber. Three biological replicates of the experiment were conducted. Samples were taken at three growth stages: coleoptile (VE); seedling with the tip of the second leaf visible ("VE/V1" stage); seedling with the first leaf fully developed (V1, when the collar of the first leaf is visible). Stages were described according to Ritchie et al. (1986). The sampling relied on developmental stage and not on time because in cold conditions, the growth of plants is retarded. Each sample consisted of pooled material from at least three plants (coleoptile for VE, first leaf for VE/ V1 and V1 stages). The samples were flash-frozen in liquid nitrogen and stored at − 80°C until further processing.

Isolation of Nucleosome-Free DNA
Nucleosome-free DNA was extracted according to the FAIRE procedure (Omidbakhshfard et al. 2014). The DNA samples were sonicated into 200-800 bp fragments. During sonication, the samples were kept on ice. The procedure consisted of ten 15-s pulses with 2-min breaks between them to prevent heating of the samples. Kontes KT-50 Micro Ultrasonic Cell Disrupter (Walker F.C. Ebel, Germany) was used with an output set at 55. Phase Lock Gel tubes (Eppendorf, Germany) were used to facilitate the retrieval of clean DNA during phenol-chloroform extraction. Precipitated DNA was dried in a vacuum centrifuge (20 min, 65°C) to ensure complete phenol removal. Samples were resuspended in 10 mM Tris-HCl, pH 7.5, treated with 1 μl DNase-free RNase A (30 min, 37°C) and then with 1 μl Proteinase K (1 h, 56°C). After overnight incubation (65°C), the samples were cleaned with DNA Clean & Concentrator-5 (Zymo Research, USA). DNA from each sample was checked for purity with a NanoDrop ND1000 spectrophotometer (Thermo Fisher Scientific, USA). The size range and concentration of the samples were quantified with the Bioanalyzer 2100 system (Agilent, USA). DNA samples (20 ng in 30 μl of nucleasefree water) were sent to Fasteris (Switzerland) for library preparation and sequencing. The libraries were prepared according to ChiP-Seq TruSeq protocol with fragmentation step and sequenced on the Illumina HiSeq 2500 instrument with 150-bp single-end read length.
FAIRE is well suited for finding candidate regulatory regions in poorly annotated complex genomes, e.g., of maize. This method is not biased by chemical DNA fragmentation, which could be problematic in DNase-seq. Additionally, FAIRE, unlike ChIP, does not require high-quality specific antibodies, which are difficult to produce (Tsompana and Buck 2014).
Statistically significant open chromatin peaks were detected with the Homer v4.9 package (makeTagDirectory and findPeaks commands; Heinz et al. 2010). Options optimized for finding nucleosome-free regions ("style factor" and "nfr") and 0.05 false discovery rate (FDR) threshold were used. To minimize false positives, we further considered only peaks present in at least two biological replicates of a given variant (temperature × maize line × growth stage), as determined by the mergePeaks command of Homer. The peaks were annotated relative to gene features (annotatePeaks.pl command of Homer). The − 3 kb to + 100 bp range around the gene transcription start site (TSS) has been defined here as the "promoter." For peaks located in this region, corresponding gene IDs were retrieved and used in all further functional analyses. The peaks were analyzed for overrepresentation of transcription factor binding sites (TFBSs) using CLOVER (Frith 2004) and a nonredundant set of plant motifs (JASPAR 2014, Mathelier et al. 2014AthaMap v7, Hehl and Bülow 2014). In the TFBS overrepresentation analysis, three background sequence sets from maize were used: chromosome 10, 26,300 bp upstream gene region, and randomly generated peaks in the 3 kb upstream gene region. To minimize the number of false positives, only TFBS motifs with score > 0 and p value < 0.05 against all three background sets were retained for further analysis, as recommended by CLOVER developers (Frith 2004).
Gene sets with cold-specific NDR in the promoter were retrieved for each line × growth stage variant. For these gene sets, the protein interaction data were retrieved from the STRING database version 10 (Szklarczyk et al. 2017). Lowconfidence interactions (score < 400) and proteins without interactants were filtered out. Each network was visualized and analyzed with Cytoscape v3.4 (Cline et al. 2007). For functional analyses, the Gene Ontology (GO) annotation of maize genes was used (Wimalanathan et al. 2018). Genes annotated with the "DNA binding transcription factor activity" (GO:0003700) GO term and their interactants were used for enrichment analysis in the BiNGO module of Cytoscape (Maere et al. 2005;Cline et al. 2007). For transcription factors specific for cold treatment in each line × growth stage variant, we searched for more data. Transcription factors with one of the following GO annotations: "response to cold" (GO:0009409), "response to freezing" (GO:0050826), "response to osmotic stress" (GO:0006970), "response to oxidative stress" (GO:0006979), and "response to salt stress" (GO:0009651), were determined. We also used our microarray results (Sobkowiak et al. 2014(Sobkowiak et al. , 2016 to check whether transcription factors were previously determined to be differentially expressed in cold-tolerant maize lines. These datasets were overlaid with networks of direct transcription factor interactions. The most connected nodes were detected by hub analysis using the Network Analyzer tool of Cytoscape.

Verification of NDR Peaks and Transcriptional Activation
Expression of nine selected genes with NDR peak in cold in S68911 plants at V1 stage was validated by RT-qPCR. mRNA was isolated from the first leaf (V1 stage) of coldtreated and control plants with NucleoSpin RNA Plant Kit (Macherey-Nagel, Germany). mRNA was reverse transcribed with the SuperScript First Strand Synthesis System (Invitrogen, USA). For each gene, reactions for all three biological replicates of each temperature variant were done, each reaction was done in triplicate. Three reference genes were used, ubiquitin (GRMZM2G118637), actin depolymerizing factor (GRMZM2G060702), and synthetic luciferase gene (cat. no. L4561, Promega, USA). Luciferase mRNA (1 ng) was added during the reverse transcription step (O'Shaughnessy et al. 2011). The primers are shown in Online Resource 1. The results were analyzed with the ΔΔCt method (Schmittgen and Livak 2008) using average expression of reference genes. Expression values (ΔCt) in cold and control were compared by t test in R v3.6.2 (R Development Core Team 2019) with p value cutoff of 0.05. Reactions were carried out in a CFX96 Touch Real-Time PCR (BioRad, USA) in 10 μl volume using 5 μl SYBR Green JumpStart Taq ReadyMix (Sigma-Aldrich, USA) with primers at 400 nM each and 2.5 μl of 1:50 diluted cDNA. The temperature profile of reactions was as follows: initial denaturation (95°C, 3 min), 40 cycles of denaturation (95°C, 5 s), and annealing/extension (60°C, 20 s).
Numerous studies have signified the link between nucleosome-free chromatin in the upstream region of the gene and the activation of gene expression (Hogan et al. 2006;Vera et al. 2014;Rodgers-Melnick et al. 2016;Oka et al. 2017). To complement RT-qPCR verification of our results, we searched public databases, ArrayExpress (Parkinson et al. 2010) and Gene Expression Omnibus (Barrett et al. 2013). Although those databases contain data from numerous studies on maize, most of them are not comparable to our data. For the S68911 line, there are data only from our previous study, which used older plants than those used in this study; therefore, these data could not be used. Our comparison was limited to the B73 line and control conditions due to the lack of the appropriate cold variant in the databases. Coleoptile samples (VE stage) from our study and external studies were compared. From the remaining studies, we selected those in which leaves of approximately 2-week-old (three leaves visible) seedlings were sampled and compared them to V1 stage samples.
Raw data were downloaded from the ArrayExpress database. Fastq files were quality-controlled with FastQC v0.11.9 (Andrews 2010). Adapters and low-quality bases were removed with Trimmomatic v0.39 (Bolger et al. 2014) and only reads at least 40 nt long were kept for further analysis. Reads were mapped to the B73 AGPv3 maize genome with tophat v2.1.1 (Trapnell et al. 2009). In cases when sequencing was done in paired-end mode, only alignments in which both reads in a pair can be mapped were kept ("-no-mixed" option). Low-quality alignments (score < 10) were filtered-out with samtools v1.10 ). Fragments per kilo base per million mapped reads (FPKM) values were computed with cufflinks v2.2.1 (Trapnell et al. 2010). For each compared project genes were divided to two sets, according to presence or absence of NDR peak in promoter region in our results.
FPKM values for these two sets were compared by t-test in R v3.6.2 (R Development Core Team 2019). Hypothesis that the true difference in FPKM means is not equal to zero was tested against the null hypothesis of equality of means.

Results
Two maize inbred lines with contrasting levels of cold sensitivity were compared in this work. Specifically, information on DNA free from nucleosomes in cold-grown seedlings was obtained. Such sequences may contain cis-regulatory elements determining the expression of downstream genes.
On average, the DNA from each experimental variant was sequenced to comparable depth. In the majority of cases, at least two biological replicates showed a high Spearman correlation (Online Resource 2, Fig. 1). The PCA plot shows that major differences between the variants can be attributed to the maize line followed by temperature conditions (Online Resource 2, Fig. 2). More peaks were detected in S68911 than in B73 samples in both control and cold conditions, the basic statistics regarding peak-calling are presented in Online Resource 3. The number of "promoter" peaks in cold conditions relative to the control increased in S68911 (coldtolerant) and decreased in the B73 (cold-sensitive) inbred line along with the growth stage (Online Resource 2, Fig. 3). Example genes with NDR peak in the "promoter" were visualized in genome browser (Online Resource 4, Fig. 4).

Verification of NDR Peaks and Transcriptional Activation
Expression of nine genes with NDR peak in promoter in S68911 line at V1 stage was checked by the RT-qPCR method. Among chosen genes were those with NDR peak unique to cold and those with NDR peak in both temperature treatments. Expression of six genes was verified successfully. Three genes have NDR peak unique to cold but expression difference between treatments was not statistically significant (Online Resource 1).
To further check the correspondence of NDR presence and transcription of downstream genes, we used public RNA sequencing data. In all comparisons, there was a significantly greater number of reads for genes with NDR peaks in promoters than for genes without NDR peaks. This finding means that genes with NDR were preferentially expressed ( Table 1).

Analysis of Regulatory Elements
Stress can alter the nucleosome occupancy of particular regulatory sequences and induce the expression of a specific set of genes. FAIRE peaks present in the − 3 kb to + 100 bp region around the TSS of genes (designated the "promoter" region) were analyzed for cis-motif (TFBS) overrepresentation with CLOVER software (Frith 2004). The majority of enriched TFBSs were nonunique for variants (temperature × maize line × growth stage). There were no TFBSs enriched specifically in the cold at all developmental stages in one line. However, close to this condition were those for five teosinte branched 1, cycloidea, proliferating cell nuclear antigen factor 1 (TCP) proteins (TCP15, TCP23, OJ1581_H09.2, OsI_08196, and ARALYDRAFT_493022) in S68911 (Table 2(A)). Similarly, in the B73 line, TFBSs for Arabidopsis NAC domain-containing protein 81 (ANAC81), indeterminate growth1 (id1) and TATA-binding protein (TBP) were enriched at all stages (Table 2(B)). During the VE and VE/ V1 stages in the S68911 line, there were many line-specific overrepresented TFBSs targeted by APETALA2/ethyleneresponsive element binding protein (AP2/EREBP) transcription factors (TFs). Such TFBSs were absent at V1 when other motifs were more common (Table 3(A)). Differently in B73, the most numerous TFBSs were targeted by MADS box TFs, which was true for the VE and VE/V1 stages. In this line, only one AP2/EREBP TFBS was enriched at VE/V1 and one at the V1 stage (Table 3(B)). For each overrepresented TFBS, it was determined if its cognate transcription factor had an NDR peak in the same sample. Therefore, transcription factor proteins would potentially be present, and their target genes would be open for regulation. This finding was observed for the S68911 line for auxin response factor 4 (ARF4) and long hypocotyl 5 (HY5) at the VE/V1 stage and for TCP11 at the Table 1 Correspondence of the presence of the nucleosome-depleted region (NDR) in the gene promoter region and transcriptional activation. The promoter was defined as the − 3 kb to + 100 bp range around the gene transcription start site. For comparison, the data from the ArrayExpress database were used (Parkinson et al. 2010

Transcription-Factor Level Analyses
For each variant (line × growth stage), a set of genes with cold-specific NDRs was defined. For these sets, protein interaction networks (STRING v10; Szklarczyk et al. 2017) have been created. These networks were used in enrichment analysis in the BiNGO module of Cytoscape (Maere et al. 2005;Cline et al. 2007). This analysis aimed to identify functional categories overrepresented among proteins directly interacting with transcription factors (sample set) against all proteins in a given network (population set). Transcription factors were selected based on the "DNA binding transcription factor activity" (GO:0003700) GO annotation (Wimalanathan et al. 2018). Categories from the Biological Process GO class were analyzed. For the S68911 line, both negative and positive regulation of gene expression was observed, which was   (Yilmaz et al. 2009) connected with chromatin changes ("chromatin silencing"). Categories related to development were enriched at the VE ("cell differentiation") and VE/V1 ("developmental process") stages. Categories related to stress were shifted forward; at VE/V1, it was "immune system process", whereas at V1, it was "response to stress" and "DNA repair" (Online Resource 6). For the B73 line, there was an enrichment for the V1 sample only. In this line, the most characteristic overrepresented categories were related to gene expression and DNA repair (Online Resource 6). The network of transcription factors for the S68911 line was larger at later stages than at VE. For this line, there were three transcription factors common to all cold variants: knotted 1 (KN1), wuschel 1 (WUS1), and GATA zinc finger family protein 17 (GATA17). In the S68911 line, GATA17 was a highly connected node at all stages, and it was the hub at the Table 3 Transcription factor binding sites (TFBS) that were significantly enriched among nucleosome-free regions located in the promoter region of genes. The results for cold-stressed maize seedlings of line, (A) S68911 (cold-tolerant), (B) B73 (cold-sensitive). Only TFBS specific for one or two consecutive growth stages in a given line in cold are shown. Further description is provided in  VE stage. At this stage, there were two transcription factors described in the literature (Bolduc et al. 2012), but not in STRING, as interacting with KN1. This finding was observed for the DELLA protein Dwarf plant 9 (D9) and the Homeobox transcription factor 102 (HB102) (Fig. 1). The network for the VE/V1 stage for the S68911 line contained few highly connected hubs, and the main hub was histone-lysine N-methyltransferase (SUVR5). This gene interacted with HY5, whose TFBS was enriched at this developmental stage. HY5 was itself connected with four heat shock transcription factors (HSFs). SUVR5 also interacted with several epigenetic and MYB proteins, and it was connected with the second hub, GATA17, interacting with MADS proteins. WUS1 interacted with ARF4, whose TFBS was enriched at this stage (Fig. 2a).
The main hub of the V1 network for the S68911 line was GATA34. This gene interacted with two hubs, GATA17 and the two-component response regulator ARR12. Directly and via three MYB factors, the central hub interacted with the basic region/leucine zipper (BZIP122). The latter was connected itself with several HSFs, some shared with the VE/ V1 stage. GATA34 also interacted with several developmental (MADS and homeotic) transcription factors. As at other growth stages, at V1, there was a subnetwork containing mainly homeotic proteins centered around KN1 and WUS1; however, at each stage, the composition of the subnetwork varied. At the V1 stage, there was also a second small network of three AP2-EREBP proteins joined by interaction with WRKY108 (Fig. 2b).
For the B73 maize line, there was a network of transcription factors for the V1 stage only. The central protein was transcription factor HY5 connected with nine other proteins which, in most cases, did not have further connections.
Among the proteins interacting with HY5 were three HSFs (Fig. 3). TFBS for HY5 was enriched at this developmental stage.
In summary, the investigated maize lines differed markedly at the level of potentially activated transcription factors. Nonetheless, the shared element was HY5, and its TFBS was enriched in both lines, although earlier in S68911 (at VE/V1 stage) than in B73 (at V1 stage). Network files are available in Online Resource 7.

Discussion
Maize exhibits a high level of diversity at both the phenotypic and genomic levels. Variability in cold sensitivity has enabled the successful cultivation of maize in the temperate climate. Although it was possible due to breeding relying on trial and error, as the molecular level causes of lower cold sensitivity of some maize lines are unknown.
In this work, the reaction of the first maize leaf to cold at the, especially poorly known, chromatin level was assessed. Specifically, potentially expressed genes and cis-regulatory motifs accessible for transcription factors were analyzed. DNA free of nucleosomes was isolated and sequenced. Three early growth stages of two inbred lines, cold-tolerant S68911 and cold-sensitive B73, were investigated.

Verification of NDR Peaks and Transcriptional Activation
In RT-qPCR reactions, upregulation of all nine genes in cold relative to control was observed. However, not always the expression change was statistically significant (t test, p < 0.05). Comparing the qPCR and FAIRE-seq results, four cases could be distinguished: (1) clearly verified genes with stronger expression in cold and NDR peak unique for cold; (2) genes with stronger expression in cold but NDR peak in both cold and control, what could take place, when in cold also other factors positively regulate expression of gene; (3) clearly verified genes with nonsignificant expression difference in cold and control and NDR peak in both conditions; (4) unverified genes without significant expression difference between treatments but with NDR peak unique for cold. Taken together, FAIRE-seq and RT-qPCR results are in good correspondence. It is supported by the good correspondence of FAIRE-seq results and RNA-seq expression data from the ArrayExpress database (Table 1). Nevertheless, the differences for some genes mean that FAIRE-seq and other chromatin-level methods need verification by expression assays. Fig. 1 Network of transcription factors potentially upregulated in S68911 maize seedlings at the VE stage in cold conditions. Genes discussed in the article text and some connecting them are labeled. Closely connected members of some transcription-factor families or functional groups are outlined in color. KN1 was shared by all S68911 networks and upregulated in transcriptomic studies. Abbreviations: D9, Dwarf plant 9; GATA17, GATA zinc finger family protein 17; HB102, Homeoboxtranscription factor 102; KN1, knotted 1; WUS1, wuschel 1; MADS18, Zea mays MADS18

Cold-Stressed Plants of the Cold-Tolerant Line Adapt the Photosynthetic Apparatus to Cold and Activate Antioxidative Systems
Potential metabolic-level differences between the S68911 and B73 maize lines were analyzed using the CornCyc database version 7 (Walsh et al. 2016). Potentially activated reactions in each pathway were counted for each experimental variant. This procedure has defined pathways that showed consistent changes in cold-treated vs. control plants in subsequent developmental stages. Among the pathway categories potentially upregulated in the cold-tolerant line grown in the cold, the most numerous were as follows: secondary metabolism, fatty acid biosynthesis, cofactor biosynthesis, reactive oxygen species (ROS) scavenging, and hormone biosynthesis (Online Resources 8 and 9).
Photosynthesis was one of the processes activated in the cold-tolerant S68911 line during cold stress (Online Resource 8). In a previous study, several maize lines were compared under the same growth and stressful temperature conditions used in this work. That study showed that the photosynthetic apparatus of the first leaf was significantly more efficient in S68911 than in B73 and other cold-sensitive maize lines (Grzybowski et al. 2019). Additionally, S68911 seedlings at the V3 stage acclimate to severe cold (< 8°C) during growth in the short-term moderate cold (10-14°C) and successfully Fig. 2 Network of transcription factors potentially upregulated in S68911 maize seedlings at the a VE/V1 and b V1 stages in cold conditions. Genes discussed in the article text and some connecting them are labeled. Closely connected members of some transcription-factor families or functional groups are outlined in color. KN1 was shared by all S68911 networks and upregulated in transcriptomic studies. Abbreviations: AGL16, agamous-like protein 16; ARF4, auxin response factor 4; ARR12, two-component response regulator ARR12; BEL1, BEL1-like homeodomain protein 3; BZIP57, transcription factor PERIANTHIA; BZIP122, bZIP transcription factor superfamily protein 122; DBB4, Bbox zinc finger protein 24; EREB18, AP2 domain-containing protein 18; EREB101, ethylene-responsive transcription factor ERF055; EREB184, AP2-like ethylene-responsive transcription factor ANT; GATA34, GATA zinc finger family protein 34; HB28, pathogenesis-related homeodomain protein; HMG103, high mobility group B protein 2; HSF28, HSF type transcription factor 28; HY5, long hypocotyl 5; IAA27, auxin-responsive protein IAA27; KN1, knotted 1; KNOX1, knotted related homeobox1; MADS18, Zea mays MADS18; MYB36, myb domain protein 109; MYB41, R2R3MYB-domain protein 41; MYB60, myb domain protein 60; MYB112, myb domain protein 112; RAP2-2, ethylene-responsive transcription factor RAP2-2; SUVR5, histone-lysine N-methyltransferase SUVR5; TFIIA, transcription factor IIA alpha/beta subunit; WRKY108, probable WRKY transcription factor 75; WUS1, wuschel 1 recover when stress ceases (Sobkowiak et al. 2016). Acclimation was shown by, among others, a photosynthetic apparatus, relatively efficient in S68911 in comparison to the cold-sensitive lines (Sobkowiak et al. 2016). Similar results were obtained for four cold-tolerant maize lines grown in cold (Verheul et al. 1996;Haldimann 1999). The positive cold response of photosynthesis-related genes in S68911 and other cold-tolerant lines is notable. Photosynthesis is regarded as the most cold-sensitive process in maize (Kingston-Smith and Foyer 2000;Sobkowiak et al. 2014) and in C4 plants in general (Long 1983).
During cold stress, the dark phase of photosynthesis is disturbed earlier, and the light phase is overexcited, generating ROS and leading to photoinhibition (Lidon et al. 2001). One of the causes is a decrease in enzymatic reaction velocity, importantly antioxidative reactions ). This decrease is stronger in C4 plants, such as maize, in which elements of photosynthetic and antioxidative systems are separated into two cell types, and transport of intermediates between them must be operational (Foyer et al. 2002). It has been shown that both short-and long-distance transport is impaired in cold-stressed maize seedlings (Sowiński et al. 2003;Bilska-Kos et al. 2016). As we have found, among the potentially activated pathways in S68911 were those related to ROS. Two of these pathways are located in chloroplasts, "detoxification of reactive carbonyls in chloroplasts" and "zeaxanthin, antheraxanthin, and violaxanthin interconversion." The former is potentially activated in S68911 at all growth stages investigated in this work. Increased activity of the antioxidative system components was shown previously in cold-tolerant maize lines under cold conditions (Hodges et al. 1997;Aroca et al. 2001). Similarly, cold-tolerant Z. diploperennis showed higher activity of antioxidative enzymes than Z. mays (Jahnke et al. 1991).

The Cold-Tolerant Maize Line Activates Stress-Response Processes to Survive Low-Temperature Conditions
Three cold signaling systems are known in plants: inducer of CBF expression-C-repeat binding factor/dehydration responsive element binding (ICE-DREB/CBF), abscisic acid (ABA)-dependent, and mitogen-activated protein kinase (MAPK) cascade (Buti et al. 2018). Many components of these pathways were potentially expressed in S68911, mostly at the VE/V1 and V1 stages (Tables 4 and 5). Buti et al. (2018) postulated that the ICE-CBF/DREB pathway accounts for the difference in cold tolerance between two rice varieties. There was virtually no activation of cold signaling pathways in the B73 line.
Several TFBSs targeted by AP2/EREBP cold-responsive transcription factors were overrepresented in S68911 seedlings grown in cold conditions (Table 3(A)). Members of this family were also present in the interaction network for the V1 stage (Fig. 2b). These findings agree with the high expression of EREB101 and related to apetala 2-2 (RAP2-2) in S68911 in the cold (Sobkowiak et al. 2016). Similar response was found for EREB101 in the cold-tolerant ETH-DH7 maize line (Sobkowiak et al. 2014). The expression level of RAP2-2 was positively related to the freezing tolerance of the two maize lines (Li et al. 2016). In the B73 line, only two AP2/EREBP cis-elements were enriched (Table 3( B)). This fact could indicate delayed and weaker activation of stress-responsive genes in the B73 line.
Other stress-responsive factors, HSFs, were present in the networks for the VE/V1 and V1 stages in the S68911 line (Fig. 2). HSF28 was present in the VE/V1 network in the cold. Expression of this gene was positively related to the cold tolerance level of three maize lines, reaching a maximal value in S68911 (Sobkowiak et al. 2016).

Hormone-Level Changes and Defense Mode in the Cold-Tolerant Line
Marked morphological differences were observed between S68911 and B73 maize seedlings grown under cold conditions (Grzybowski et al. 2019 and this study). This fact prompted us to perform a detailed analysis of potential regulation at the hormone level in cold-stressed maize. The analysis of CornCyc pathways showed that in the S68911 line in the cold conditions the most dynamic changes concerned the metabolism of gibberellins (GAs) with both their inactivation and synthesis. As for stress-related hormones, ABA was synthesized at all stages investigated, and jasmonic acid (JA) was synthesized at the VE/V1 and V1 stages (Online Resource 8).
To complement the analysis of CornCyc pathways, which do not comprise hormone signaling, the Plant Reactome database release 15 (Naithani et al. 2017) was used. In a number of Fig. 3 Network of transcription factors potentially upregulated in B73 maize seedlings at the V1 stage in cold conditions. Genes discussed in the article text and some connecting them are labeled. Closely connected members of some transcription-factor families or functional groups are outlined in color. Abbreviations: DBB8, B-box zinc finger protein 22; HSFs, heat-shock transcription-factors; HY5, long hypocotyl 5 Table 5 Maize homologs of cold-response elements described in rice (Buti et al. 2018) Table 4 S68911 B73 Rice gene name Gene ID Gene name VE_C VE/ V1_C  cases, a single reaction in the pathway could be performed by more than one gene product. To account for this fact, a set of genes responsible for a single reaction is collectively referred to as a "node." For each node in each signaling pathway, the number of genes with NDR peaks in the promoter was counted (Online Resource 10). GA and JA signaling are mutually antagonistic by the action of DELLA proteins. The balance of these hormones determines defense or growth response (Colebrook et al. 2014). The potential prevalence of expression of JA synthesis genes (Online Resource 8) and JA signaling over GA signaling (Online Resource 10) observed in S68911 seedlings suggests their switching to defense mode (Borrego and Kolomiets 2016). A similar reaction has not been observed in B73. Stresses that slow the metabolism, including cold, cause degradation of active GAs. It was postulated that this degradation helps the plant survive (Shan et al. 2013;Colebrook et al. 2014). It was shown that S68911 seedlings were not sensitive to externally applied GA3, what was measured as mesocotyl length. Conversely, GA3 had significant positive effect on B73 seedlings. (Grzybowski et al. 2019). Nevertheless, coldstressed sensitive and tolerant seedlings of rice showed negative and positive GA response, respectively (Buti et al. 2018). The JA-deficient opr7opr8 maize double mutant has longer mesocotyl and coleoptile ). This observation links JA and the xeromorphic phenotype observed in cold stress in S68911 and other cold-tolerant maize lines (Verheul et al. 1996;Sowiński et al. 2003;Strigens et al. 2013;Grzybowski et al. 2019).
Potential activation of JA synthesis and signaling in the S68911 line is in keeping with research on cold-stressed O. sativa (Du et al. 2013). This hormone upregulated the ICE-CBF/DREB cold signaling pathway (Hu et al. 2013). JA increased the cold tolerance of rice seedlings (Lee et al. 1996), and more tolerant lines had an elevated level of this hormone (Buti et al. 2018). JA and cold stress are linked by the chloroplast membrane, which is a site of both JA synthesis and cold reception (Chinnusamy et al. 2007;Borrego and Kolomiets 2016).
JA and ABA synthesis genes were potentially expressed in cold in the cold-tolerant line. It was shown that both hormones induce JA-synthesizing genes (Zhang et al. 2005). Additionally, the ABA-dependent cold signaling pathway is potentially activated, as also have been shown by Zhao et al. (2015). ABA is critical in the maintenance of water relations in maize . The level of this hormone increased in maize lines grown in the cold, regardless of their tolerance level (Ristic et al. 1998;Janowiak et al. 2002). It was shown that treatment with ABA decreased the cold injury of maize seedlings (Aroca et al. 2003;Pál et al. 2011).
In the auxin signaling pathway, the sharpest differences between the lines were visible at the level of auxin reception, potentially activated in all variants of S68911 predominantly in the cold. The potential expression of ARF and transport inhibitor response 1/auxin signaling F-box (TIR1/AFB) suggests the activation of the response to auxin in the coldtolerant line. The pathways are presented in Online Resource 10, and source data are given in Online Resource 11.
It should be noted that the presented metabolic pathway analyses are limited because not all maize genes are annotated and included in pathways (hence absence of ABA, ethylene, and cytokinin signaling in our results).

Homeotic Genes Could Have a Role in the Cold-Stress Response
In the context of hormonal regulation, two potentially expressed genes are notable, KN1 (GRMZM2G017087) and WUS1 (GRMZM2G047448). A subnetwork composed of mainly homeotic genes clustered around KN1 and WUS1 was present at all S68911 stages (Figs. 1 and 2). Potential expression of KN1 is in concert with its nearly fivefold upregulation in cold in ETH-DH7 maize line (Sobkowiak et al. 2014). According to the literature, KN1 and WUS1 are expressed predominantly in the shoot meristem, and KN1 also in generative organs and developing stem (Jackson et al. 1994;Stelpflug et al. 2016). The primary function of KN1 is shoot meristem maintenance by repressing GA synthesis (Kerstetter et al. 1997). The expression levels of KN1 and WUS1 in microarray studies on maize stress response were checked using the Gene Expression Omnibus database (Barrett et al. 2013). GEO2R module and a false discovery rate corrected p value cutoff of 0.05 were used. KN1 expression was activated above twofold during drought stress in leaves (Hayano-Kanashiro et al. 2009) and shoots (Zheng et al. 2010). Similar droughtresponse of maize was shown by Nelissen et al. (2018). The expression of KN1 was also more than fourfold upregulated in maize stems challenged by the European corn borer (Dafoe et al. 2013).
Data on WUS1 expression are scarce. Nevertheless, in rice seedlings, orthologous OsWUS was highly upregulated after 12 h of drought but not after cold or salinity stress (Cheng et al. 2014). These results suggest that KN1 and WUS1 may have a role in the stress response, but its elucidation requires further study. In keeping with growth inhibition and potential GA degradation, KN1 was shown to activate the GA catabolism gene gibberellin 2-beta-dioxygenase 1 (GA2OX1, Bolduc et al. 2012). Moreover, in the network for VE stage, KN1 interacted with a repressor of GA signaling, Dwarf plant 9 (D9, Lawit et al. 2010).
KN1 was shown to repress the auxin signaling pathway (Bolduc et al. 2012), potentially activated in the leaves of cold-stressed S68911 plants. Auxin and KN1 act antagonistically in plant development, what complicates the interpretation of our results. However, a positive role of auxin in the cold response was found in cold-tolerant  and sensitive rice varieties (Jain and Khurana 2009;Du et al. 2013), and cold-tolerant cultivar of the C4 plant Digitaria eriantha (Garbero et al. 2012).

Conclusions
The chromatin-level analysis suggests that the superior performance of the S68911 maize line in cold is dependent on the defense-favoring partitioning of resources and avoidance of water stress. This avoidance is manifested on the phenotype level by the xeromorphic growth form (Grzybowski et al. 2019), which was also found for other cold-tolerant maize lines. It can be hypothesized that the underlying mechanism involves a balance of GAs and JA in favor of the latter. This change would enable seedlings to arrest growth, secure resources, and resume growth after the end of cold stress. The role of sustained developmental processes in cold tolerance is supported by our previous molecular and physiological level analyses (Sobkowiak et al. 2016;Grzybowski et al. 2019). This study concerned the most basic level of gene expression regulation; therefore, it requires further confirmation. Nevertheless, we presented several candidates for further studies on cold stress tolerance in maize and possibly in related species.