We used two near isogenic lines (paired NILs) that were genetically identical, differing only for the presence or absence of a chromosome 9 introgression from S. habrochaites acc. LA1778 in an otherwise S. lycopersicum cv. T5 genetic background. Development of the NILs (designated NIL175 and NIL163) and their graphical genotypes are described in Easlon et al. (2013). Briefly, NIL175 contains a S. habrochaites introgression on chromosome 9 spanning markers T1670 to T0532, an interval which includes QTL stm9, and this NIL maintains shoot turgor (i.e., exhibits tolerance) when exposed to root chilling. NIL163 does not contain any introgressions from S. habrochaites and is chilling sensitive, exhibiting shoot wilting when its roots are exposed to chilling stress (Easlon et al., 2013).

All plant flats and hydroponic tanks used in our experiments were located in a greenhouse at UC Davis in Davis, California, with ambient daytime conditions of 25–37°C and 55–80% relative humidity, and nighttime conditions of 18–25°C and 20–55% relative humidity. Seeds of NIL175 and NIL163 were seeded into 72-cell flats containing SuperSoil potting media (Rod McLellan Co.), watered daily, and fertilized with a 10:30:20 NPK solution 10 days after planting. After 2 weeks, roots of seedlings of NIL175 and NIL163 were carefully washed free of potting media in deionized water and transferred to a hydroponic growth tank containing 20% strength modified Hoagland's solution (Epstein and Bloom, 2005). Plants were maintained in the tank containing the aerated, circulating root solution at 20°C under ambient illumination. After 1 week, plants were transferred and randomized (see Section Experimental Design and Sample Collection) in a separate refrigerated hydroponic tank at 20°C containing fresh 20% strength modified Hoagland's solution. Plants were acclimated overnight prior to initiation of the root chilling experiment the next day.

A completely randomized design was employed for conducting a time-course experiment, with two genotypes (NIL175 and NIL163), five sampling time-points (0, 1, 1.5, 2, and 4 h from initiation; see description below), and three biological replicates per time-point/genotype combination. Experimental conditions were similar to those used in our prior experiments on shoot turgor maintenance under root chilling (Truco et al., 2000; Goodstal et al., 2005; Arms et al., 2015). At the initiation of the experiment, the temperature setting on the hydroponic tank was decreased from 20 to 6°C, and a supplemental cooling wand was placed in the solution to aid chilling until the target temperature of 6°C was reached at 2 h. At each of the five time-points, three plants of each NIL were collected (for a total of 30 samples). A table describing each sample is available (see Supplementary Table 1). Plants were removed from the hydroponic tank, roots briefly and quickly dried on absorbent paper towels, and whole plants placed in pre-labeled and autoclaved aluminum foil envelopes before being flash frozen in liquid nitrogen. Samples were transported to the lab in liquid nitrogen, transferred to a −80°C freezer, and stored at −80°C until processed for mRNA extraction (see Section Library Preparation and Sequencing). Several additional plants of NIL175 and NIL163 were also included in the experiment to observe and confirm the expected phenotypes at 4 h. Those are: shoot wilting in NIL163 and maintenance of shoot turgor in NIL175 (Arms et al., 2015).

The five time-points were selected to capture biologically important time-points in the plant stress response to root chilling, including the transition from no stress (above 10°C) and the start of chilling stress (below 10°C; Geisenberg and Stewart, 1986; Allen and Ort, 2001). At 0 h, the plants were at ambient greenhouse temperatures with a root solution temperature of 20°C. Plants collected at 0 h served as non-stressed controls for all subsequent time-points. It should be noted that changes in other environmental conditions such as light quantity and quality during the course of an experiment performed over time may affect transcriptional regulation of some genes. At 1 h, the root solution decreased to ~11–12°C, just above the 10°C threshold for root chilling stress. At 1.5 h, the root solution decreased to ~8–9°C, below the 10°C threshold for root chilling stress. At 2 h, the root solution reached 6°C, where it was maintained for the duration of the experiment. Susceptible NIL163 showed first signs of wilt (loss of turgor in leaf tips) at 2 h. At 4 h, NIL163 was fully wilted, while NIL175 maintained shoot turgor. In our prior experiments, phenotyping for shoot turgor maintenance was performed at 4 h (Truco et al., 2000; Goodstal et al., 2005; Arms et al., 2015). After sample collection and phenotyping at 4 h, the experiment was terminated. For a more in depth description of phenotypic responses of the NILs and experimental conditions, we refer the reader to Arms et al. (2015).

Approximately one quarter of the root mass of each sample (30 in total) was ground in liquid nitrogen to a fine powder with a mortar and pestle. Following pulverization, 100 mg of each sample was processed using the Qiagen RNeasy mini kit (www.qiagen.com). Isolated RNA was quantified using Life Technologies Qubit RNA Broad Range Assay kit (www.lifetechnologies.com), and quality control was performed with the Agilent Bioanalyzer RNA 6000 Nano kit (www.agilent.com). Messenger RNA (mRNA) from each sample was isolated from 10 μg of total RNA as input using the Bioo Scientific NEXTflex Poly(A) Beads isolation kit (www.biooscientific.com). Isolated mRNA was quantified using the Life Technologies Qubit RNA High Sensitivity Assay kit (www.lifetechnologies.com).

Sequencing mRNA libraries were prepared for each sample using the Bioo Scientific NEXTflex Rapid RNA-Seq kit with NEXTflex RNA-Seq barcodes (1–30 were used from the 48 barcode set; www.biooscientific.com; see Supplementary Table 1). Library quality control was performed on the Agilent Bioanalyzer with the DNA High Sensitivity DNA kit (www.agilent.com), and quantified via qPCR at the UC Berkeley Vincent J. Coates Genomics Sequencing Laboratory (GSL) (http://qb3.berkeley.edu/gsl/). All barcoded libraries were pooled and sequenced at the GSL on 12 lanes of single-end 50 sequencing on the Illumina Hi-Seq 2000 (www.illumina.com). The barcodes enabled sorting of the resulting reads from each sequencing lane with CASSAVA v1.8.2 (www.illumina.com) at the GSL before further analysis.

All read data were uploaded to and analyzed on the iPlant Discovery Environment using supported bioinformatics tools (https://de.cyverse.org/de/). For each of the 30 samples, all reads were concatenated into a single file, and quality control was performed. Scythe adapter-trimming was used to remove 3′-end adapter contamination, and Sickle quality-based-trimming was used to remove low quality reads. Read mapping statistics are provided (see Supplementary Table 1; Buffalo, 2010; Joshi and Fass, 2011). FastQC 0.10.1 (multi-file) was then used to check the quality of the read collection. Reads from each library were aligned to the S. lycopersicum cv. Heinz 1706 reference genome (SL2.50; solgenomics.net) with GSNAP (Genomic Short-read Nucleotide Alignment Program) using default settings and the quality format option Sanger (Supplementary Table 1; Wu and Nacu, 2010). GSNAP output SAM files were converted to BAM files using SAM-to-BAM, and then BAM-to-Counts was used to quantify the number of reads mapping to each S. lycopersicum ITAG v2.40 annotated transcript in FPKM (fragments per kilobase of transcript per million mapped reads). Number of reads mapped to each transcript per sample are available in Supplementary Table 2. To identify differentially expressed (DE) transcripts (genes), pairwise comparisons of different time-point/NIL combinations were made using DEseq1.0, with a P = 0.05 significance threshold. Pairwise comparisons were done between NILs at each time-point, as well as between the 0 h control and each subsequent time-point within each NIL. DEseq controls for false discovery rate due to multiple testing using the Benjamini-Hochberg procedure. All significantly DE genes were identified using this corrected P-value (Anders and Huber, 2010).

The datasets supporting the conclusions of this article are included within the article and its supplementary files. The raw Illumina sequences are publicly available via the NCBI Sequence Read Archive as of January 1st, 2017, under the following accession number: SRP081139.

Clustering of transcript expression patterns across the five time-points for NIL175 and NIL163 was performed using Mclust v5.0.1 (www.stat.washington.edu/mclust/) in R64 v3.2.0 (www.r-project.org/). Mclust uses a mixed method that combines model-based hierarchical clustering, EM for Gaussian mixture models, and BIC to optimize clustering results, including number of clusters (Fraley et al., 2012). The log average of FPKM values across the three biological replicates of a NIL at each time-point for each gene was used as input. After Mclust analysis, trend lines for each cluster were produced by averaging the values for all genes in that cluster at each time-point. Within each NIL, trend lines were grouped according to transcription pattern to produce “Trend Groups,” each of which included all individual clusters with the same relative transcriptional pattern regardless of the level of transcription. Trend Groups may therefore represent more than one cluster, but each cluster will only be included in a single Trend Group. The cluster and Trend Group assignment for each gene (for both NIL175 and NIL163) are included in a Supplementary File (Supplementary Table 2). Trend Groups were then compared between the two NILs. Trend Groups identified in both NIL175 and NIL163 were denoted with capital letters A–M, and Trend Groups identified in a single NIL were indicated with lower case letters n–v. Trend Group A consisted of a single cluster for each NIL, which represented no detectable level of transcript at any time-point (i.e., FPKM = 0).

Results of the differential expression and Trend Group analyses were visualized via area-proportional Venn diagrams using BioVenn (www.cmbi.ru.nl/cdd/biovenn/). Visualizations and comparisons between NIL175 and NIL163 for DE genes identified in time-point comparisons to the 0 h control were used to identify overlapping and unique DE genes. BioVenn was also used to visualize and compare gene lists from NIL175 and NIL163 of Trend Groups A–M. Trend Groups (n–v) unique to either NIL175 or NIL163 were not included in Venn diagram visualization. Output lists of unique and overlapping genes from Venn diagram analyses were used as input gene lists for GO enrichment analysis using g:Profiler vr1435_e80_eg27 (biit.cs.ut.ee/gprofiler/). S. lycopersicum ITAG v2.40 was specified, and default settings were used, with the exception of a user-defined P-value of 0.05 to return only statistically significant results (used to produce tables provided as supplementary files due to size). The default g:SCS (Set Counts and Sizes) threshold method was used to control for false positives (inflated type I error) from multiple testing. This correction method takes into account the effects of multiple tests and GO terms that are not independent, but are a hierarchically organized collection of related general and specific terms (Reimand et al., 2007). Gene lists that returned significantly enriched GO terms were rerun using the “Best per parent group” option, where “parent” was defined as the highest-level GO term in a hierarchy of related GO terms (used to produce in-text tables). This setting returned results that grouped related GO terms and represented them with the most statistically significant parent GO term to aid in parsing of the data.

Differential expression analysis identified DE genes for all time-point/NIL combinations, with generally increasing numbers of DE genes as the experiment progressed. A list of differentially expressed genes for all time-point/NIL combinations is available in a table (see Supplementary Table 3). In comparisons between the NILs at each of five time-points, the number of DE genes ranged from 20 at 0 h to 63 at 4 h. Chromosome 9 had the greatest number of DE genes, ranging from 13 at 0 h to 22 at 2 h. Interestingly, DE genes located across all 12 chromosomes were identified only at 2 h (Supplementary Table 3).

DE genes were identified within each NIL by comparing subsequent time-points to the 0 h control (Supplementary Table 3). In both NILs, the number of DE genes increased as the experiment progressed, and DE genes were located across all 12 chromosomes. In tolerant NIL175, DE genes increased from 371 at 1 h to 2677 at 4 h. In susceptible NIL163, DE genes increased from 228 at 1 h to 2238 at 4 h. In NIL163, only one DE gene was located within the introgressed region, and only at 4 h (Supplementary Table 3). Venn diagrams revealed overlap between the NILs for DE genes when time-point DE gene lists were compared (Figure 1). A majority of DE genes in NIL175 were shared with NIL163 only at 1.5 and 4 h. In contrast, the majority of DE genes for NIL163 were also DE in NIL175 for all time-points. NIL175 also had more DE genes than NIL163 at all time-points except at 1.5 h (Figure 1).

GO enrichment analysis for the unique and overlapping gene lists from the time-point to time-point comparisons between the NILs identified significantly enriched GO terms for the majority of gene lists. The best by parent grouping and the full results of GO term enrichment are available in Table 1 and Supplementary Table 4. All overlapping gene lists from each time-point comparison had significantly enriched GO terms, which increased in number as the experiment progressed. The overlapping gene lists were predominantly enriched for GO terms related to protein modifications or transcriptional regulation. At 4 h, the comparison also showed enrichment for GO terms related to photosynthesis and plant response to light.

Significantly enriched GO terms were identified for unique gene lists of NIL175 at all time-points. Similarly to the overlapping gene lists, enriched GO terms were related to protein modification and regulation. In addition, GO terms related to cell communication were identified at 1 h, and response to nutrient levels (specifically nitrate and nitrate transport) were identified at 2 h. NIL163 had significantly enriched GO terms in the unique gene lists only at 1 h (nucleic acid binding) and 4 h (nucleosome and response to blue light). Unique parent GO terms identified at 1 h for both NIL175 and NIL163 were significantly enriched in the list of overlapping genes at 1.5 h. Interestingly, tolerant NIL175 had enrichment of GO term GO:0004674 (protein serine/threonine kinase activity) starting at 1.5 h and through the end of the experiment. This GO term was also present in the overlapping gene list at 2 h, and was not represented in any other overlapping DE gene lists, or for any unique gene lists of NIL163 (Table 1, Supplementary Table 4). It is also interesting to note that no GO terms directly related to ABA signaling were significantly enriched, including those for membrane transport or signaling pathway regulation (Table 1, Supplementary Table 4).

Cluster analysis resulted in 48 unique expression pattern clusters for NIL175 and 146 unique clusters for NIL163. Cluster assignments for each gene are available for each NIL (see Supplementary Table 2). Visual examination revealed that clusters were based on both transcriptional pattern (i.e., changes in transcription across time-points), as well as general transcript level. For example, two genes with the same transcript expression pattern but a three-fold difference in relative level of transcription were sorted into different clusters. Cluster analysis of NIL175 and NIL163 did not show any evidence of genes that are expressed only under chilling stress conditions but not at the 0 h (non-stressed) control (Supplementary Table 2), as exemplified by the lack of clusters composed of genes with zero or near-zero expression levels at 0 h, followed by expression at subsequent time points as chilling stress was imposed.

Comparisons of Trend Groups revealed 13 were common to both NILs (designated Trend Groups A–M; Figure 2). In addition, NIL175 had one unique Trend Group (designated n), and NIL163 had eight unique Trend Groups (designated o–v). Figures showing each trend line are provided for A–M (see Figure 2) and n–v (see Supplementary Figure 1). For both NILs, Trend Groups A and C contained the largest number of genes, although the number of overlapping genes within each shared Trend Group (A–M) varied greatly (Figure 2). Trend Group A (no expression) and Trend Group C (very slight decrease to 1.5 h, peak at 2 h, and decrease at 4 h) had almost completely overlapping gene sets between the NILs. Trend Group E (increase to 1.5 h, decrease at 2 h, and highest level of expression at 4 h) contained only 13 genes in NIL175, all of which were represented in the larger gene list for NIL163 (Figure 1, Supplementary Table 2). For most Trend Groups, NIL175 had a larger set of genes than NIL163, with the exception of Trend Group E (mentioned previously), and Trend Group J (peaks at 1 and 2 h).

All Trend Groups (with the exception of Trend Groups B and L) had significantly enriched GO terms in at least one of the three gene lists (unique to NIL175 or to NIL163 or overlapping). The best by parent grouping GO term analyses are given in Table 2. Full results of GO term enrichment are provided in an additional table (see Supplementary Table 5). The number of significantly enriched GO terms mirrored the number of genes in a list, with more populous gene lists containing more significantly enriched GO terms. Trend Group A was enriched for GO terms relating to DNA metabolism and transcription, including the transcriptional complex Mediator (GO:0016592), and DNA helicase (GO:0003678). Trend Group C was enriched for GO terms related to transcription factor regulation, protein modification, and transport. Trend Group E was unique in that it was the only Trend Group for which significantly enriched GO terms were identified for the unique set of genes in NIL163, but not the overlapping list or unique list for NIL175 (Table 2, Supplementary Table 5).

Multiple gene networks have been identified in crop plant responses to abiotic stresses, including changes in stress signaling, protein metabolism, energy metabolism, ROS scavenging enzymes, and photosynthetic enzymes (Allen and Ort, 2001; Des Marais and Juenger, 2010; Kosová et al., 2015). In our study, we looked at genome-wide patterns of transcriptional regulation in a set of paired NILs that differed in their tolerance to rapid-onset water stress induced by root chilling. We observed overlap between root chilling-tolerant NIL175 and chilling susceptible NIL163 for significantly enriched GO terms related to cell communication, transcription factor activity, various metabolic processes, response to light, and response to oxygen-containing compounds. In gene lists unique to NIL175, significantly enriched GO terms were detected for cell communication, protein modification, and response to nutrient levels, while in NIL163 many of the significantly enriched GO terms were related to photosynthesis (Tables 1, 2, Supplementary Table 4). Although there are few published studies on transcriptional regulation in tolerant and susceptible tomato lines exposed to chilling stress, our results were generally consistent with, but also distinct from, other studies. Liu et al. (2012) compared gene expression in S. lycopersicum acc. LA4024 with wild S. habrochaites acc. LA1777 and LA3969 (a cultivated tomato line containing an introgression from S. habrochaites on chromosome 12). Whole plants were exposed to chilling stress at 4°C, and leaf samples were collected each day for 7 days. The authors reported that GO terms related to response to stimulus, response to stress, and metabolism were significantly enriched across all three lines. Metabolism and stress response-related GO terms were significantly enriched in tolerant LA1777 and LA3969, while in susceptible LA4024, defense response and circadian rhythm terms were significantly enriched (Liu et al., 2012). In a study on a chilling-tolerant tomato species, S. lycopersicoides, Chen et al. (2015a) identified GO term enrichment for response to stimulus, signaling, and cell death in leaves after whole plant exposure to chilling stress at 4°C compared to non-stressed controls. GO enrichment analysis provides a translation of a gene list into a functional profile of biological processes, cell components, and molecular functions overrepresented in the gene list. It may be negatively affected by the quality of the genome annotation, as well as the quality of the GO term annotation and association (Khatri and Draghici, 2005; Reimand et al., 2007; Gene Ontology Consortium et al., 2013). Therefore, it is possible that additional biological, cellular, and molecular processes involved in the root response to chilling-stress are not captured in this analysis.

An interesting discovery in our analysis was the lack of clusters in either NIL that showed evidence of transcriptional activation in response to root chilling (Table 2, Supplementary Table 4 and Supplementary Figure 1), and the significant enrichment of GO terms related to Mediator and DNA helicase in Trend Group A (no expression) for both NILs (Table 2, Supplementary Figure 1). Both Mediator and DNA helicase are essential for transcriptional activation of genes that are not already primed for transcription in eukaryotes, including plants (Lohman and Bjornson, 1996; Allen and Ort, 2001; Patel and Donmez, 2006; Wu and Nacu, 2010; Hentges, 2011; Wu, 2012; Poss et al., 2013). Our findings suggest that while de-novo transcription may be occurring, the response to chilling stress in roots may be more dependent on transcriptional regulation of genes already primed for expression, as well as post-transcriptional/post-translational modifications of transcripts or proteins already present in root cells. The importance of post-transcriptional and post-translational regulation of proteins in regulating transcriptional changes caused by osmotic stress has been reported previously and was recently reviewed by Guerra et al. (2015).

When plants are exposed to osmotic stress, a major mechanism of stress signaling within cells is increased enzymatic activity of kinases, referred to as osmotic stress-activated kinases (Kulik et al., 2011; Fujii and Zhu, 2012). The phytohormone abscisic acid (ABA) has been identified as fulfilling a critical role in coordinating plant response to osmotic stress in conjunction with various kinases and other regulatory molecules (Fujita et al., 2013). The serine/threonine kinase type 2 (SnRK2) family of kinases (plant equivalents of the Snf1 kinases in yeast) are predominantly involved in the cellular response to osmotic stress in both an ABA-dependent and ABA-independent manner (Coello et al., 2011; Yunta et al., 2011). SnRK2 kinases regulate ABA-responsive elements (ABRE) such as ABRE-binding proteins and ABRE-binding factors, transcription factors that regulate the ABRE-mediated regulation of downstream targets. Members of the SnRK2 family also have been shown to control stomatal closure in response to osmotic stress (Yunta et al., 2011; Fujita et al., 2013). The significant enrichment of protein serine/threonine kinase activity (GO term GO:0004647) in tolerant NIL175 in response to root chilling stress, but not in susceptible NIL163 (except for the overlapping gene set at 2 h), suggests that the SnRK2 kinase pathway plays a significant role in the NIL175 response to root chilling. This is further supported by the significant enrichment of GO:0004647 in only the unique gene list for NIL175 for Trend Group C in the Trend Group analyses. As our analyses did not show evidence of enrichment for ABA-specific GO terms, it is possible that the signaling of root-chilling stress in the roots occurs in an ABA-independent manner. These results do not preclude the involvement of ABA in the shoot response to root chilling stress, which was not investigated in these experiments.

The enrichment of protein serine/threonine kinase activity (GO:0004647) is also interesting in the context of the significant enrichment in NIL175 of GO terms related to cell communication at 1 h and response to nutrient levels at 2 and 4 h. Two other classes of SnRK kinases (SnRK1 and SnRK3) have been implicated in stress signaling in addition to more primary roles in cellular metabolism, including carbon, nitrogen, sucrose, and lipid metabolism (Laurie and Halford, 2001; Hey et al., 2010; Kulik et al., 2011). It is possible that SnRK kinases and their associated pathways may be serving multiple functions during plant exposure to root chilling, initially serving as part of the signaling cascade sensing and responding to root chilling stress, and eventually as modulators of metabolism throughout the duration of the stress episode.

Our previous research documented that the physiological response to rapid-onset water stress induced by root chilling causes distinct phenotypic responses in susceptible and tolerant tomato lines that differ for the presence or absence of a S. habrochaites introgression on the short arm of chromosome 9 (Truco et al., 2000; Bloom et al., 2004; Goodstal et al., 2005; Easlon et al., 2013; Arms et al., 2015). In this study we discovered that while the roots of both NILs exhibited some overlap in their transcriptional response to root chilling, unique responses were also evident. Plant transcriptional responses to abiotic stresses are complex and overlapping, and a core set of stress response genes are rapidly induced in response to a range of stresses (Kültz, 2005; Walley and Dehesh, 2010). In general, plant responses to abiotic stresses can be divided into overlapping responses that are exhibited in response to a variety of stresses, vs. specific responses that result in re-establishment of homeostasis in a manner that is specific to the imposed stress (Kültz, 2005). Overlap observed in the DE time-point comparisons between the NILs as well as in Trend Groups A–M may represent a general response to stress that is common to multiple types of abiotic stress, and may also be shared across tomato species. Transcriptional patterns not shared by the two NILs may be the result of each NIL responding to root chilling in a genotype- and stress-specific manner.

A pattern of overlapping and unique transcriptional responses has also been observed in other studies of tolerant and susceptible tomato genotypes exposed to chilling or drought stress (Gong et al., 2010; Chen et al., 2015b). Chen et al. (2015b) compared transcriptional regulation in leaves of chilling-tolerant S. habrochaites LA1777 to sensitive S. lycopersicum cv. Glamour exposed to 4°C for 0, 1, and 12 h. The majority of gene regulation patterns were dissimilar, yet some genes showed overlapping patterns in both susceptible and tolerant genotypes. Gong et al. (2010) reported similar transcriptional patterns as Chen et al. (2015b) when comparing S. lycopersicum cv. M82 with tomato lines containing S. pennellii introgressions on chromosome 2 (line IL2-5) and chromosome 9 (line IL9-1, with an introgression syntenic to the introgression in our NIL175) exposed to drought stress for several days. The findings of these two studies are in parallel with our observations of a combination of unique and overlapping transcriptional patterns in our two NILs. Furthermore, these studies compared a tomato cultivar to either a wild species or to cultivated tomato containing a wild species introgression, thus the observed overlapping transcriptional responses imply that cultivated tomato shares portions of general stress response pathways with its wild relatives but without exhibiting the same phenotypic tolerances to some abiotic stressors.

We had hypothesized that chilling-tolerant NIL175 would exhibit a more complex transcriptional response to root chilling in comparison to susceptible NIL163, however a larger number of clusters was identified for NIL163 than for NIL175 (Supplementary Table 2). This implies greater transcriptional complexity in chilling-sensitive NIL163, while the greater enrichment of significant GO terms (general and stress-specific) in unique gene lists for NIL175 in both the time-point comparisons for DE genes as well as for the Trend Groups imply a more focused transcriptional response to root chilling stress in NIL175 compared to NIL163 (Tables 1, 2, Supplementary Table 4 and Supplementary Figure 1). Significantly enriched GO terms represent an over-enrichment of genes with similar functions or within the same network for a given input gene list (Khatri and Draghici, 2005; Reimand et al., 2007). A similar result was reported by Gong et al. (2010) in that drought-susceptible S. lycopersicum cv. M82 exhibited a larger number of DE genes in leaves in response to drought stress than did tolerant tomato lines containing S. pennellii introgressions. Interestingly, although M82 showed a larger number of DE genes, the tolerant lines had a higher number of stress-responsive transcription factors that were differentially expressed. This result is similar to our findings that susceptible NIL163 had a more complex pattern of transcriptional regulation, and tolerant NIL175 exhibited a more targeted transcriptional response to root chilling stress. Research by Liu et al. (2012) on chilling stress in tomato leaves also showed a more focused transcriptional response in tolerant wild S. habrochaites LA1777 and in LA3696, a tomato line containing an S. habrochaites introgression in comparison to susceptible S. lycopersicum LA4204 after whole plant exposure to 4°C.

Of the 34,725 annotated genes in the S. lycopersicum cv. Heinz 1706 reference genome ITAG v2.40 (SGN, solgenomics.net), only 20 genes were differentially transcribed between our paired NILs at the non-stressed 0 h control (Supplementary Table 3). The very small number of DE genes in the non-stressed control suggests that there was little effect of the S. habrochaites chromosome 9 introgression on global transcription under non-stress conditions in our experiment. Interestingly, the majority of the 20 DE genes at 0 h were located on chromosome 9, although outside of the introgressed region. Previous work in Arabidopsis thaliana and Thellungiella salsuginea indicates that tolerant genotypes tend to exhibit increased levels of stress-inducible and stress-responsive proteins even in the absence of stress (Pang et al., 2010; Wendelboe-Nelson and Morris, 2012; Kosová et al., 2015). This suggests that the S. habrochaites introgression present in NIL175 may lead to a priming of the abiotic stress response pathways involved in root chilling tolerance, and that many of these response genes appear to be located on chromosome 9. Since there was an absence of enriched GO terms for DE genes between the two NILs at 0 h, it is not possible to interpret whether specific functional gene groups or networks were preferentially targeted. This evidence, along with the results discussed under “GO terms significantly enriched in NIL175 and NIL163 in response to root chilling” (above) suggests that while the chromosome 9 introgression from S. habrochaites may prime NIL175 to respond to root chilling stress in a tolerant manner in our experiment, it did not fundamentally alter transcription of specific functional gene networks, or genome-wide transcription levels.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.