Data Set
To generate a comprehensive data resource combining experiments involving heat, drought and salinity stress on tomato cultivars, we searched the NCBI Sequence Read Archive (SRA) for “Solanum lycopersicum AND heat”, “Solanum lycopersicum AND drought”, “Solanum lycopersicum AND salt”, respectively, and filtered the resulting entries according to the following criteria: reads contain total mRNA, samples were taken from non-mutant/wildtype plants, growth, treatment, and sampling conditions were sufficiently described to determine the genotype and tissue, plants were not subjected to a combination of multiple stresses, for heat stress: The heat stress temperature was at least 32°C and plants were grown in (semi-)controlled conditions to ensure consistent exposure to the set temperature. 23 data sources fulfilled the criteria, to which we added 3 previously unpublished heat stress datasets from 3 different experiments: a comparison of different leaf ages and temperature levels performed at the IPB Halle, a time series of seedling grown on media including day and night cycles performed at the IGZ Großbeeren, and an experiment on ovaries in the greenhouse performed at the Hebrew University of Jerusalem (HUJI); growth and sampling conditions for these datasets are described in the methods section. An overview of the used data sources can be found in Table 1.
Table 1
Data Sources used for the meta analysis. Journal Publications were frequently not linked to the data set on NCBI, they were manually added.
Data Source ID
|
NCBI BioProject Accession
|
GEO Accession
|
Reference
|
DS-1
|
PRJNA657834
|
GSE156402
|
-
|
DS-2
|
PRJNA484882
|
-
|
(Zhou et al. 2019)
|
DS-3
|
PRJNA376497
|
-
|
(Liu et al. 2017)
|
DS-4
|
PRJNA644865
|
-
|
(Diouf et al. 2020)
|
DS-5
|
PRJNA869132
|
-
|
-
|
DS-6
|
PRJNA635375
|
GSE151277
|
(Q. Wang et al. 2021)
|
DS-7
|
PRJNA800740
|
-
|
-
|
DS-8
|
PRJNA745152
|
-
|
(Nie et al. 2023)
|
DS-9
|
PRJNA746070
|
-
|
-
|
DS-10
|
being submitted to NCBI
|
(data from HUJI)
|
this publication
|
DS-11
|
being submitted to NCBI
|
(data from IPB)
|
this publication
|
DS-12
|
PRJNA624892
|
GSE148530
|
(T. Zhao et al. 2020, 2021)
|
DS-13
|
being submitted to NCBI
|
(data from IGZ)
|
this publication
|
DS-14
|
PRJNA560638
|
-
|
(Wen et al. 2019)
|
DS-15
|
PRJNA603594
|
-
|
(Almeida, Perez-Fons, and Fraser 2021)
|
DS-16
|
PRJNA562700
|
-
|
(Batista et al. 2020)
|
DS-17
|
PRJNA639840
|
GSE152620
|
(Lopez-Delacalle et al. 2021)
|
DS-18
|
PRJNA753098
|
-
|
-
|
DS-19
|
PRJNA769714
|
GSE185583
|
(Jansma et al. 2022)
|
DS-20
|
PRJNA730730
|
GSE174607
|
-
|
DS-21
|
PRJNA319398
|
GSE80556
|
(Keller et al. 2017)
|
DS-22
|
PRJNA559982
|
-
|
(Kashyap et al. 2020)
|
DS-23
|
PRJNA476376
|
-
|
(Keshishian et al. 2018)
|
DS-24
|
PRJNA888477
|
-
|
(B. Wang et al. 2023)
|
DS-25
|
PRJNA282940
|
GSE68500
|
(Fragkostefanakis et al. 2016)
|
DS-26
|
PRJEB42497
|
-
|
(Mellidou et al. 2021)
|
The data set contains 338 samples: 92 from heat stressed plants, 72 from drought, 37 from salt, and 137 from respective controls. They stem from 8 different tissues and 22 different genotypes. A complete list of the used samples and their annotated metadata can be found in Table S1 and the interactive shiny app accompanying this manuscript. After comparing stress treated replicates with their respective controls and calculating fold-changes and p-values for differential expression, we end up with 88 contrasts, 39 of which pertain to heat, 36 to drought, and 13 to salt stress (see Figure Fig. 1A). Data is dominated by leaf contrasts (58) followed by seedling and fruit with 10 contrasts each. Only a minority of data sources represent anther, ovaries, pollen, seed, and root (Fig. 1A). In the following sections we will refer to the dataset of normalized log-transformed transcripts per million (TPM) as the logTPM data, and to the fold changes to respective controls as fold change data.
Data source is the strongest factor determining both gene expression and transcriptomic stress response
Analysing the variance components of the logTPM data revealed that the largest known source of variance was the experiment/data source, accounting for 28.6% of the variance in gene expression. An additional 11.9% of the gene expression in our data set was determined by the genotype, 7.7% by tissue and only 3.7% by the treatment (control, heat, drought, or salt stress). The remaining 48.1% of variance could not be explained by these factors. In the first two dimensions of a principal component analysis (PCA), representing 47.0% of the variance, most samples from the pollen, ovary and anther tissues placed distantly to the other samples, even from the same genotype, reflecting their different expression profile compared to vegetative tissues (Fig. 2A, this effect has also been observed in Arabidopsis (Klepikova et al. 2016) although not in pollen and ovary tissue specifically). Samples from seedling, root, fruit, and seed co-localised or formed a sub-cluster within the leaf samples. In contrast, we observed almost no separation by treatment, and no following principal component (PC) could be attributed to a general separation of the different treatments (Figure S1 and shiny app).
The same holds for the fold-change data, where the experiment/data source was still the largest known source of variation (27.6%), followed by treatment (10.9%), tissue (8.0%), and genotype (4.9%). The remaining 48.6% of data variance was not explained by the model. In the first two components of the PCA (representing 22.0% of the variance), the tissues formed tighter clusters but neither tissue nor treatments separated clearly (Fig. 2B). Accordingly, none of the following PCs could be associated with the treatment and tissue factor (Figure S2 and shiny app).
These results indicate that the large majority of responses observed in any individual heat, drought or salt stress experiment is specific to the given experimental setup.
Heat stress triggers surprisingly diverse transcriptomic responses
In general, heat stress led to extensive changes in gene expression, affecting more genes than drought stress but fewer than salt stress. The number of differentially expressed genes (DEG, using an FDR threshold of < 0.05) between treatment and control in heat varied between 2187 (contrast DS-10-1, equivalent to 6.3% of all quantified transcripts) and 14376 (contrast DS-11-5-1, equivalent to 41.44% of all quantified transcripts) with a median number of 8788 changes (Fig. 3A). This compares to a median of 4229 and 11347 DEGs in drought and salt stress, respectively. No correlation with the length of treatment or treatment intensity could be observed for the heat stress contrasts (Fig. 3A).
To describe the gene expression response on the functional level, we translated the DEG results into up- and down-regulation of functional groups (gene ontology terms; GO terms). GO term analysis for each contrast revealed enrichment of 151 to 1128 GO terms among the upregulated genes (median 656) and 437 to 1888 GO terms among the downregulated genes (median 1158; adjusted p-value < 0.05). Again, no correlation with the length of treatment or treatment intensity could be observed (Fig. 3B). The GO term “response to heat” (GO:0009408) was enriched in the upregulated genes of 26 out of 28 contrasts so the criteria for compiling the dataset were successful in selecting experiments which caused heat stress. Also consistently enriched were the terms “protein folding” (GO:0006457), “response to high light intensity” (GO:0009644; both in 28 out of 28 contrasts), “post-embryonic development” (GO:0009791) and “reproductive process” (GO:0022414; each in 27 contrasts) and “response to oxidative stress” (GO:0006979), “response to ethanol” (GO:0045471) and “cellular catabolic process” (GO:0044248) in 26 out of 28 contrasts.
To establish the size of the total potential heat stress response, we selected all genes significantly affected in at least one of the examined heat stress contrasts (p-value < 0.05, adjusted for all comparisons). Surprisingly, this yielded 26,133 genes, more than two thirds of the 34,688 protein coding genes in the genome. The genes cover a variety of biological processes and mechanisms, of which some but not all are directly linked to heat stress (Table S2). The majority of the affected genes (16,788 of 26,133) were significantly upregulated in some contrasts but downregulated in others.
Examining the data on drought and salt stress experiments the same way, we obtained analogous results and found considerable overlaps with the heat stress gene set as well as the enriched GO terms, representing potential targets downstream of the crosstalk between the stresses (Fig. 3C and D).
These results confirm that a majority of tomato stress responses is activated in interaction with additional factors, such as the tissue, the genotype, the duration or the intensity of the stress. Therefore, extraction of core stress-specific responses is not straightforward. Thus in the next step, we searched for a strategy to identify the most likely elements of such a core stress response and evaluated the number of independent experiments and statistics needed to provide a reasonable estimate.
Combining experiments selects more meaningful candidate genes than ordering by p-value or π score
When faced with a large number of differentially expressed genes, one often has to prioritize and filter out genes for follow-up experimental functional characterisation. One simple strategy is ranking DEGs by p-value and choosing from among the genes with the lowest p-value. But statistical significance does not always reflect biological impact. To alleviate this, some methods also take into account the magnitude of the fold change such as the π significance score developed in (Xiao et al. 2012) which multiplies the log10-transformed p-value with the absolute fold change value. In terms of meta-analyses, gene candidate selection relates to identification of responses that are most likely to be confirmed in multiple independent experiments.
To explore these strategies for this dataset, we calculated the number of genes selected for a given p-value and π score cutoff for each contrast in our dataset (Fig. 4A 1–2). Expectedly, in both cases the number of potential candidates quickly decreased with decreasing p (or π) value cutoff irrespective of whether heat, drought or salt stress was examined. As a benchmark for how representative the selected genes are for the phenomenon of interest, we then determined the proportion of those candidates that are annotated with the primary GO term for this stress type or any of its children terms (GO:0009408 response to heat, GO:0009414 response to water deprivation, GO:0009651 response to salt stress). For comparison, we did the same with an nonspecific term which was comparably enriched in the union set (GO:0009203 ribonucleoside triphosphate catabolic process, GO:0042454 ribonucleoside catabolic process, GO:0098542 defense response to other organism), Fig. 4A 4–5. For the drought stress and salt stress treatments, the proportion of genes representative of the stress response did not change much, showing that this strategy was not very useful to identify the characteristic “core response” to that stress. The enrichment with the nonspecific term did not change either. In contrast, for heat stress, both strategies increased the proportion of representative genes. However, the confidence intervals of these curves are quite large, indicating that this observation was true for some but not consistent across all experiments. Overall, lowering the p-value cutoff reduced likely false positives at least in some heat stress experiments but was not a generally reliable strategy and while including the fold-change via a π score smoothed out the curves, it did not improve the results by much.
As an alternative approach for meta-analyses, we tested intersecting DEGs from multiple experiments by selecting only those genes which are significantly over- or under-expressed in at least 80% of the examined contrasts of the respective stress type (Fig. 4A 3). As before, the number of selected genes quickly decreased as more contrasts were added with most genes being eliminated in the first 5–10 contrasts. However, unlike the previously described strategies, the selected genes became more and more enriched in the respective stress-describing GO term for all three examined stresses with relatively small confidence intervals, meaning the specificity in the gene sets was increasing with more data (Fig. 4A 6). This suggests that carrying out a treatment repeatedly in different laboratories under slightly varying conditions facilitates the identification of a more representative core response, addressing the limitations associated with a single-experiment analysis, which were especially visible in our drought and salt stress datasets.
Defining the core heat stress response as differentially abundant genes in ≥ 80% of our 28 contrasts, we find 46 up- and 11 downregulated genes which are responsive consistently across the different tissues, genotypes, and data sources (Table 2, Fig. 4B). The biological processes enriched in these upregulated genes are much more specific: response to heat (GO:0009408), protein folding (GO:0006457), response to oxidative stress (GO:0006979), and response to light intensity (GO:0009642) are the top hits (Table S3), the last probably because there is substantial cross-talk between heat and light regulation (Franklin 2009; Legris et al. 2017).
About one third of the upregulated genes encode heat shock proteins and other chaperones (e.g. HSP21 and BAG6-like*; see Table 2 for full gene description and ID), on top of that there are genes encoding regulatory proteins and transcription factors (SlC3H68, VA1*), enzymes (GGCT*, APX*), and transmembrane transporters (COPT*, NPF1.2*), as well as some hypothetical and previously uncharacterized proteins (labeled as unknown; BLASTp results for them in Table S4). The downregulated genes include those for fatty acid desaturases and enzymes important for cell wall biosynthesis (FAD*, UGE4*), transporters (SIVGT2, MGTRNIPA*) and regulatory proteins (CBP*, RSMG*). Twelve of the upregulated genes have a predicted binding site for heat shock factors (Fig. 4B). Eight of them encode heat shock proteins, two are unknown proteins and the last two are ZFP* and GST*.
There are few overlaps between the genes in this core set and the analogously derived core sets for drought and salt stress (Fig. 4C). The up- and down-regulated genes for each stress are clearly separated and the stresses largely distinct. Of the three stress types, drought and salt stress seem to share more of their transcriptional response with each other than with heat stress.
Table 2
Heat Stress Core genes (transcript IDs and gene descriptions from the ITAG4.1 annotation). Αbbreviations were chosen to more easily refer to genes within this manuscript. Where possible, existing symbols from UniProt or Solgenomics were used. Otherwise abbreviations were chosen based on the gene description and marked with an asterisk throughout (*) to show they do not represent confirmed gene names.
Effect of HS
|
Transcript ID
|
Gene Description
|
Abbreviation (Source)
|
UPREGULATED
|
Solyc01g095320.4.1
|
Bag family molecular chaperone regulator 6-like protein
|
BAG6-like*
|
UPREGULATED
|
Solyc01g102960.3.1
|
22.0 kDa class IV heat shock protein
|
HSP22-1*
|
UPREGULATED
|
Solyc02g080130.5.1
|
Chaperone protein dnaJ 50
|
DnaJ50*
|
UPREGULATED
|
Solyc03g082420.3.1
|
Small heat shock protein, chloroplastic
|
HSP21 (UniProt)
|
UPREGULATED
|
Solyc03g113930.3.1
|
22.0 kDa class IV heat shock protein
|
HSP22-2*
|
UPREGULATED
|
Solyc04g011440.4.1
|
ethylene-responsive heat shock protein cognate 70
|
ER21 (Solgenomics)
|
UPREGULATED
|
Solyc04g014480.3.1
|
15.7 kDa heat shock protein, peroxisomal
|
HSP15.7*
|
UPREGULATED
|
Solyc05g014280.4.1
|
small heat shock protein 1
|
VIS1 (Solgenomics)
|
UPREGULATED
|
Solyc06g011370.4.1
|
Chaperone protein ClpB
|
ClpB*
|
UPREGULATED
|
Solyc06g036290.3.1
|
heat shock protein 90
|
HSP90 (Solgenomics)
|
UPREGULATED
|
Solyc06g076520.1.1
|
Class I heat shock protein
|
HSPI-1*
|
UPREGULATED
|
Solyc06g076560.2.1
|
Class I heat shock protein
|
HSPI-2*
|
UPREGULATED
|
Solyc08g062340.3.1
|
17.9 kDa class II heat shock protein
|
HSP17.9*
|
UPREGULATED
|
Solyc08g078700.3.1
|
Mitochondrial small heat shock protein
|
MTSHP (Solgenomics)
|
UPREGULATED
|
Solyc09g015000.4.1
|
Class I heat shock protein
|
HSPI-3*
|
UPREGULATED
|
Solyc11g020040.2.1
|
heat shock protein 70
|
HSP70 (Solgenomics)
|
UPREGULATED
|
Solyc11g066100.2.1
|
Heat shock protein 70
|
HSP70-2*
|
UPREGULATED
|
Solyc11g071830.2.1
|
DnaJ protein like
|
DnaJ-like*
|
UPREGULATED
|
Solyc01g104740.3.1
|
Multiprotein-bridging factor 1c
|
MBF1c-ER24 (Solgenomics)
|
UPREGULATED
|
Solyc02g030450.3.1
|
Protein ROOT HAIR DEFECTIVE 3 homolog
|
RHD3*
|
UPREGULATED
|
Solyc02g079020.4.1
|
B3 domain-containing transcription repressor VA.1
|
VA1*
|
UPREGULATED
|
Solyc03g063760.4.1
|
DNA-binding protein SMUBP-2
|
SMUBP-2*
|
UPREGULATED
|
Solyc03g083190.3.1
|
gamma-irradiation and mitomycin c induced.1
|
GMI1*
|
UPREGULATED
|
Solyc03g117440.3.1
|
RNA helicase DEAD14
|
SIDEAD14 (Solgenomics)
|
UPREGULATED
|
Solyc04g054290.3.1
|
Pentatricopeptide repeat-containing protein
|
PCMP*
|
UPREGULATED
|
Solyc05g055010.4.1
|
RNA-binding protein 42
|
RBP42*
|
UPREGULATED
|
Solyc06g076670.3.1
|
RNA-binding (RRM/RBD/RNP motifs) family protein
|
RNBP*
|
UPREGULATED
|
Solyc08g005420.3.1
|
cell division cycle 20.2, cofactor of APC complex-like
|
CDC20.2*
|
UPREGULATED
|
Solyc08g067960.3.1
|
CHY-type/CTCHY-type/RING-type Zinc finger protein
|
ZFP*
|
UPREGULATED
|
Solyc10g006610.5.1
|
Auxin-responsive GH3 family protein
|
AUXGH3*
|
UPREGULATED
|
Solyc11g066830.2.1
|
Zinc finger transcription factor 68
|
SlC3H68 (Solgenomics)
|
DOWNREGULATED
|
Solyc06g082710.5.1
|
KH domain-containing protein
|
KHP*
|
DOWNREGULATED
|
Solyc07g006830.4.1
|
ribosomal RNA small subunit methyltransferase G
|
RSMG*
|
DOWNREGULATED
|
Solyc12g055930.3.1
|
Calmodulin-binding protein
|
CBP*
|
UPREGULATED
|
Solyc01g065840.4.1
|
F-box/RNI-like superfamily protein
|
F-box-like*
|
UPREGULATED
|
Solyc02g062130.5.1
|
Ferredoxin-NADP reductase, chloroplastic
|
FNR*
|
UPREGULATED
|
Solyc03g113590.3.1
|
5-formyltetrahydrofolate cycloligase
|
AF516366 (Solgenomics)
|
UPREGULATED
|
Solyc08g006150.3.1
|
Gamma-glutamylcyclotransferase
|
GGCT*
|
UPREGULATED
|
Solyc09g007270.3.1
|
Ascorbate peroxidase
|
APX*
|
UPREGULATED
|
Solyc10g084400.2.1
|
Glutathione S-transferase
|
GST*
|
DOWNREGULATED
|
Solyc01g102950.3.1
|
Lycopene beta/epsilon cyclase protein
|
LCYB*
|
DOWNREGULATED
|
Solyc04g080850.3.1
|
Thioredoxin
|
TXN*
|
DOWNREGULATED
|
Solyc08g063090.2.1
|
Fatty acid desaturase
|
FAD*
|
DOWNREGULATED
|
Solyc09g091510.3.1
|
chalcone synthase 1
|
CHS1 (Solgenomics)
|
DOWNREGULATED
|
Solyc10g011810.4.1
|
Delta(8)-fatty-acid desaturase
|
D8FAD*
|
DOWNREGULATED
|
Solyc10g012410.3.1
|
UDP-glucose 4-epimerase 4
|
UGE4*
|
UPREGULATED
|
Solyc01g101070.3.1
|
ABC transporter-like
|
ABCTR-like*
|
UPREGULATED
|
Solyc08g066940.4.1
|
Protein NRT1/ PTR FAMILY 1.2
|
NPF1.2*
|
UPREGULATED
|
Solyc09g011410.4.1
|
Copper ion transmembrane transporter
|
COPT*
|
DOWNREGULATED
|
Solyc01g103760.3.1
|
magnesium transporter NIPA (DUF803)
|
MGTRNIPA*
|
DOWNREGULATED
|
Solyc03g096950.4.1
|
Vacuolar glucose transporter 2
|
SIVGT2 (Solgenomics)
|
UPREGULATED
|
Solyc03g097960.1.1
|
Cysteine/Histidine-rich .1
|
unknown-1*
|
UPREGULATED
|
Solyc06g150138.1.1
|
Unknown protein
|
unknown-2*
|
UPREGULATED
|
Solyc07g065410.1.1
|
hypothetical protein
|
unknown-3*
|
UPREGULATED
|
Solyc09g008970.1.1
|
plant/protein
|
unknown-4*
|
UPREGULATED
|
Solyc10g081480.2.1
|
maltase-glucoamylase, intestinal protein
|
unknown-5*
|
UPREGULATED
|
Solyc11g070060.2.1
|
hypothetical protein
|
unknown-6*
|
Tissue specific responses
After establishing the core heat stress response with 57 genes in total, we hypothesized that genes in sub-functionalized gene families escape the intersect criterion but cannot be considered noise. To investigate this, we employed Orthogonal Partial Least Squares-Discriminant Analysis (OPLS-DA) with the aim of identifying genes whose expression profiles most effectively discriminate between different tissues under heat stress. However, the analysis failed to identify any orthogonal components with significant predictive capability, indicating that the data set does not contain enough data to determine the tissue specific variance component of the heat response.
Potential for breeding and the heat stress core response
The core genes involved in the tomato heat stress response are spread across all 12 chromosomes, with notable clustering in certain regions, particularly towards the end of chromosome 2 and the start of chromosome 9 (Fig. 4E). (Ayenan et al. 2019) reviewed four mapping studies on heat tolerance traits and combined them into 13 meta quantitative loci (from now on referred to as QTLs for simplicity). These traits include heat stress effects on pollen viability, flower and inflorescence development as well as photosynthetic parameters. Only 17 of our identified core heat stress response genes co-localized with these QTLs. Individual examples include: BAG6-like*, putatively a modulator of chaperone activity of HSP70 (Corduan et al. 2009), in MQTL 1.4 (chromosome 1) for the anther length and heat-induced membrane damage (measured as relative electrical conductivity (Wen et al. 2019)) traits, and VA1* and DnaJ50* in the region of MQTL2.2 (chromosome 2) for the same two phenotypes. The DnaJ proteins indeed have been indicated to increase heat tolerance in tomato while overexpressed (Kong, Deng, Wang, et al. 2014) and affect the efficiency of Photosystem II in chilling stress (Kong, Deng, Zhou, et al. 2014). The remaining chromosomal regions enriched with core genes could be valuable targets for future research into heat stress resilience mechanisms in tomato.
Heat Stress Response Genes are Ancient
To verify whether heat stress response is, as we hypothesized, an evolutionarily conserved process, we determined the most recent common ancestor (MRCA) of all orthologs found in the Orthologous Matrix (OMA) database for each of our heat stress core genes. Genes whose orthologs have a more recent MRCA likely diverged more recently in evolutionary history, i.e. are evolutionarily younger and potentially less conserved across diverse species. We performed the same analysis for the drought and salt stress core genes as well as a baseline for comparison of 100 genes randomly sampled from the tomato genome (repeated 5 times, Fig. 5A). In comparison to the baseline, all core responses are enriched in genes which are evolutionarily older and likely more conserved across the tree of life. Out of the three stress responses, the heat stress genes seem to have the deepest evolutionary history with more than half of the heat core stress genes having orthologs outside of the Viridiplantae. Heat shock proteins have their predicted root at many different levels (from LUCA up until the Mesangiospermae) while several regulators and transporters (COPT*, VA1*, NPF1.2*), as well as the unknown genes mostly appear later in evolutionary history.
If the genes involved in heat stress response are conserved and ancient, it could be expected that their functional interactions, i.e. the heat stress response process, is equally conserved across the tree of life. To test this, we analyzed the co-evolution of the heat stress core response genes using the HogProf algorithm for phylogenetic profiling (Moi et al. 2020) first only within the Streptophyta and then across all of life (Fig. 5B). Strong co-evolution between two proteins (whenever one is present, the other is present as well) indicates involvement in the same biological process and a mutual functional dependency of the two proteins. Within the Streptophyta clade, four small clusters form with a perfect co-evolution score: FAD* and D8FAD* are predicted to be paralogs with the duplication event having happened on the level of the Spermatophyta, HSP22-1* and HSP22-2* are paralogs duplicated at the level of lamiids, and ER21 and HSP70-2* paralogs duplicated at the level of Viridiplantae. HSP70 is not marked as a paralog of ER21 or HSP70-2* and HSPI-1* and HSPI-2* map to the same entry on OMA, hence they have identical co-evolution with other genes and perfect co-evolution with each other. Expectedly, these clusters are equally strong across the whole tree of life. However, within the Streptophyta, two larger clusters emerge linking also non-paralogous genes (e.g. D8FAD* to CHS1, GGCT* to VIS1). These clusters only partially hold up when examining the co-evolution across all of life, indicating a part of core stress response genes evolves under a plant clade-specific selective pressure. Interestingly, the uncharacterized genes 4 and 5 are part of such a cluster as well.
Tools for meta analyses
We provide this analysis as a template for researchers interested in performing similar transcriptomic meta analyses; therefore all code for the analyses used here is designed in an easily reusable way and freely available at https://github.com/NAMlab/tomato-rna-meta. We also provide an rnaseq-mapper tool to quantify the transcript abundances of many RNA sequencing samples by mapping to a common reference transcriptome. The rnaseq-mapper is able to process a combination of both local FASTQ files and published sequence files on the NCBI Sequence Read Archive (single- or paired-end each) which it automatically downloads and combines with the local files to produce one coherent dataset. The tool is designed to be scalable also for large datasets by utilizing Nextflow to automatically parallelize tasks across CPU cores on a local machine or nodes in a computing cluster, and it is economical with disk space use by deleting data that is no longer required while running. The tool is available at https://github.com/NAMlab/rnaseq-mapper and includes all required software dependencies in an Apptainer/Singularity container which simplifies installation and execution to one line each.