Sulfate starvation has a major impact on leaf and root growth in early stages of tomato development
As previously reported for other tomato cultivars [59–64], we found that the growth of Solanum lycopersicum cv. Moneymaker seedlings was severely impacted by sulfate starvation. As shown in Fig. 1A, lack of external sulfate has a clear negative effect on the growth of the aerial tissue, as compared with the full nutrient condition. This effect is noticeable from the third week after sowing. Multifactorial analysis of leaf fresh weight and total leaf area showed that there is a strong interaction between sulfate availability and time factors (p-value <0.001, two-way ANOVA, Fig. 1B and Figure S1), indicating that the growth response of leaves to sulfate availability depends on the plant developmental stage. Accordingly, no significant differences between control and sulfate-starved plants were observed at 2 weeks after sowing (p-value=0.81), while a strong dependence of external sulfate on leaf fresh weight and area was observed 3 and 4 weeks after sowing (p-value <0.001) (Fig. 1B and Figure S1). The same interaction between sulfate and time was observed for root weight, with differences between control and sulfate-starved plants being significant only after 3 and 4 weeks of sowing (p-value <0.001, Fig. 1C). These results indicate that whole-plant growth is severely affected by the lack of external sulfate from 3 weeks after sowing.
Consistent with leaf weight and leaf area measurements over time, there is a steady accumulation of sulfate in leaves, that is abolished in sulfate-starved plants in weeks 3 and 4 (Fig. 1D). In contrast, the increase in root weight over time of plants grown in control condition is not explained by a concurrent increase in root sulfate levels (Fig. 1E). Moreover, no major differences in sulfate accumulation are observed in the roots of full nutrient or sulfate-starved plants (Fig. 1E). Together, these results indicate that additional factors besides sulfate content might be controlling root growth response to sulfate starvation.
Sulfate starvation promotes major changes in leaf and root transcriptomes in a time-dependent manner
In order to gain insights into the molecular mechanisms underlying leaf and root growth repression in response to sulfate starvation, we performed a transcriptomic analysis at a temporal and organ scale, using RNA-seq. To perform this analysis, total RNA from roots and leaves was extracted at 2, 3 and 4 weeks after sowing, considering three independent biological replicates. The RNA-Seq data was pseudo-aligned to the ITAG3.20 transcriptome annotation using kallisto , and normalized expression data in transcripts per million (TPM) was obtained (Table S1).
Principal component analysis (PCA) of the RNA-seq data showed that samples are well separated along the first component, which explains more than 50% of the variability in leaf and root samples (Fig. 2). Interestingly, full nutrient and sulfate-starved samples from leaves are grouped together at 2 weeks after sowing, indicating similar transcript expression profiles at this time point (Fig. 2). In contrast, control and sulfate-starved samples from root tissue at 2 weeks are well separated, suggesting that sulfate starvation alters the plant transcriptome in an organ- and time-dependent fashion, with an earlier response in root tissue.
Accordingly, no genes were found differentially expressed between control and sulfate-starved leaves at two weeks after sowing, while 438 genes were differentially expressed in roots in this time point (log2 Fold Change>1 and q-value<0.05) (Figure S2A and Table S2). In addition, 55% of these 438 genes are subsequently regulated in leaves (241 genes, Figure S2B), which point out an orchestrated whole organism response to sulfate deprivation over developmental time.
A considerable number of genes were differentially expressed by sulfate starvation in roots and leaves: 7024 genes at 3 weeks after sowing and 6163 at 4 weeks after sowing, with most genes being differentially expressed in leaves (Figure S2A and Table S2). This major reprogramming of the plant transcriptome correlates with the observed differences in leaf and root growth of control and sulfate-deprived plants between weeks 3 and 4, suggesting changes in gene expression are partly responsible for the observed phenotypes (Fig. 1A and B).
In order to identify genes that are differentially expressed in response to sulfate starvation, time or a sulfate starvation-time interaction, we performed a multifactorial analysis of RNA-seq data in each tissue using sleuth . In the case of roots, 6052 genes were significantly affected by time, 3755 genes by sulfate starvation and 714 by the interaction of both factors (q-value <0.05, and log2FC >1) (Fig. 3A). In the case of leaves, 6084 genes were affected by sulfate starvation, 7062 genes by time and 4571 genes by the interaction of both factors (q-value <0.05, and log2FC >1) (Fig. 3B). In contrast to root tissue, these results indicate that the response to sulfate starvation in leaves has a strong dependence on the plant age. In fact, about 90% of sulfate-responsive genes in leaves are significantly modulated by time (Fig. 3B).
To validate the results of the RNA-seq analysis with an alternative RNA quantification methodology, we selected 10 sulfate-regulated genes in root and leaves involved in different molecular and biological functions including sulfate metabolism, transport and transcriptional regulation, as well as a group of genes with unknown functions (Table S3). The expression patterns obtained by qPCR analysis showed a significant positive correlation with RNA-seq data (P <0.0001; 0.89 Pearson correlation value; Figure S3), indicating that most of genes analyzed showed a similar expression pattern than the ones identified by RNA-seq analysis.
Comparative analysis of the tomato and Arabidopsis sulfate starvation-responsive transcriptome
In a previous study, we identified 2046 genes regulated by sulfate availability in Arabidopsis thaliana using an integrative meta-analysis of transcriptomic data . This meta-analysis included samples of roots and leaves at different stages of development (from seedling to juvenile plants). In the case of tomato plants, multifactorial analysis identified 7589 sulfate-responsive genes in roots and leaves. In order to compare the sulfate starvation responsive transcriptome of both Arabidopsis and tomato, we first performed an orthology analysis of sulfate responsive genes. The identification of orthologous groups allows cross-referencing of genes from multiple species . We used the PLAZA 4.0 database  to identify the ortholog group to which each sulfate-regulated gene belongs. Thus, we identified 4262 orthologous gene families out the 7589 sulfate-responsive genes in tomato and 1338 orthologous gene families out the 2046 sulfate-responsive genes in Arabidopsis (Table S4). We found that 70% of the Arabidopsis orthologous gene families that are associated with sulfate-responsive genes were also present in tomato plants (Fig. 4A). This result indicates that the RNA-seq analysis in tomato plants captured most of the sulfate-responsive gene families previously detected in Arabidopsis.
Over-representation analysis of biological functions shared by both species revealed several GO terms that are expected to be conserved in the sulfate response as “cellular oxidant detoxification”, “sulfate transmembrane transport” or “sulfur utilization” (Table S5). Interestingly, we also found a significant enrichment (adjusted p-value <0.05) in GO terms that have not been functionally analyzed in previous studies in the context of the sulfate starvation response such as those related to cell wall or cytokinin transport (Table S5). On the other hand, specific GO terms of the Brassicaceae family such as GSL biosynthesis are among the most over-represented biological processes in the case of sulfate-responsive genes in Arabidopsis (Fig. 4B and Table S5). The intersection between the lists of orthologue genes of both species also detected sulfate-responsive processes that are found exclusively in tomato plants. The most over-represented biological process in this group of genes was "cellular response to phosphate starvation" (Fig. 4C and Table S5), suggesting that sulfate deficiency also affects internal phosphate levels in tomato plants. Specifically, we found 52 genes in this biological process including SPX genes , phosphatases and genes encoding enzymes involved in phosphate mobilization by membrane lipid remodeling such as UDP‐sulfoquinovose synthase [70, 71] (Figure S4). In order to determine the response of this set of phosphate-related genes to sulfate starvation in tomato plants, we computed the average fold of change between sulfate starved and control samples in both tissues and we then performed a hierarchical clustering analysis. We found that 60% of the genes involved in the phosphate starvation response (30 out 52 genes) showed lower expression levels in sulfate-starved plants at 3 and 4 weeks after sowing in both tissues (Cluster 3, Figure S4), suggesting that sulfate starvation alters internal phosphate levels. To address this hypothesis, we determined the internal phosphate levels at 3 weeks after sowing in both tissues. The results shown in Fig. 4D, indicate that sulfate-starved plants exhibit significantly higher phosphate levels than control plants in both tissues, with a higher phosphate accumulation in leaves.
Temporal dynamics of gene expression of tomato roots and leaves are altered by sulfate starvation
In order to get further insights about the temporal patterns of expression of the genes and associated biological processes affected by sulfate starvation, we computed Pearson correlation indexes for each pair of sulfate and time-regulated genes in roots and leaves. We performed a hierarchical clustering analysis using Dynamic TreeCut , and identified six co-expression clusters for genes regulated by sulfate and time in roots and leaves (Fig. 5A, Fig. 5B, Figure S5, Figure S6). We found that the majority of genes (60% in roots and 72% in leaves) are contained in Clusters 1 and 2 (Fig. 5A and Fig. 5B). We thus decided to analyze these clusters in more detail.
As shown in Fig. 5C and Fig. 5D, the expression profiles of Clusters 1 and 2 in roots and leaves show similar trends. Genes in Cluster 1 are repressed over time, with an earlier decrease in mRNA levels under sulfate starvation conditions, between weeks 2 and 3. On the other hand, the expression of genes in Cluster 2 slightly increases between weeks 3 and 4 in control conditions, while an important increase in gene expression is shown between weeks 2 and 3 under sulfate starvation conditions (Fig. 5C and Fig. 5D). Although the expression patterns of the genes contained in Clusters 1 and 2 are similar between roots and leaves, the identity of these sulfate-responsive genes differs between these organs (Table S6 and Table S7). As such, we found different enriched biological processes associated with Clusters 1 and 2 in roots and leaves. In roots, Cluster 1 is enriched in genes associated to growth, such as “epidermal cell division” and “plant-type secondary cell wall biogenesis”, as well as processes related to redox activity, such as “hydrogen peroxide catabolic process” and “cellular oxidant detoxification” (Fig. 5C and Table S8). In leaves, we found an enrichment of genes involved in photosynthesis, microtubule-based movement and several GO terms related with cell division such as “mitotic cell cycle checkpoint” or “mitotic spindle organization” (Table S9). GO terms related to photosynthesis are especially abundant in this cluster, so we analyzed in more detail the distribution of these genes across different categories of primary metabolism using Mapman4 annotation framework . In the case of photosynthesis, most genes are involved in light reactions, followed by the Calvin-Benson cycle and sucrose/starch metabolism (Figure S7A). Cluster 1 genes are also involved in other aspects of primary metabolism such as nitrogen assimilation and amino acid biosynthesis as well as fatty acid biosynthesis (Figure S7A). A similar result was obtained when all differentially expressed genes by sulfate starvation in leaves were analyzed using Mapman4 (Figure S7B). On the other hand, genes related to photosynthesis were under-represented in roots (Figure S7C), while other metabolic processes such as cell wall, lipid or amino acid metabolism showed a similar representation in roots and leaves (Figure S7B and Figure 7C). In the case of Cluster 2, and consistent with the marked induction in gene expression triggered by sulfate starvation shown in this cluster, we found an overrepresentation of genes involved in biological processes related to the cellular response to S starvation such as APR genes (Solyc02g032860, Solyc02g080640 and Solyc03g031620), RESPONSE TO LOW SULFUR (LSU) genes (Solyc03g096760, Solyc03g096770, Solyc03g096780) or sulfate transporters (Solyc05g054740, Solyc04g054730, Solyc04g072760, Solyc12g056930) (Fig. 5C and Table S8). Cluster 2 in leaves is enriched in GO terms related with the response to stress and protein phosphorylation (Fig. 5D and Table S9). Interestingly, we did not find an overrepresentation of S-starvation-related processes in Cluster 2 in leaves, as observed in roots. In fact, classical genes involved in the sulfate starvation response such as SULTR, ATPS, APR or LSU belong to Cluster 6 (Figure S6), suggesting an organ-specific expression of these genes in response to external S during development.
In summary, sulfate starvation promotes significant changes in the temporal expression patterns of a large number of genes in roots and leaves, of which many are associated with functions related to cellular growth, photosynthesis and sulfate uptake and metabolism. Moreover, the comparative analysis of major gene co-expression clusters between roots and leaves indicates an organ-specific response of the tomato plants to sulfate starvation.
Identification of transcription factors involved in the sulfate starvation response in tomato
Our results point at major transcriptomic changes occurring in tomato roots and leaves in response to sulfate starvation. As described in Arabidopsis, some of these changes may be attributed to transcriptional regulation by TFs. Although some TFs controlling sulfate responses have been identified in Arabidopsis , to date there is no information about sulfate-responsive TFs and their role in responses to sulfate starvation in tomato. In order to identify key TFs controlling the response to sulfate starvation of tomato, we constructed gene regulatory networks for roots and leaves using co-expression data and regulatory information on TF-DNA interaction obtained from PlantRegMap  .
In the case of roots, we obtained a gene regulatory network consisting of 24 TFs and 172 targets (Figure S8 and Table S10). TFs of the root regulatory network belong to 17 different families (Table S11), with MYB and bZIP TFs being the most abundant families (4 and 3 members, respectively). In the same way, we constructed a TF-target network using the identified sulfate-regulated genes in leaves which was composed by 45 TFs and 566 targets (Figure S9 and Table S12). Specifically, ERF, TCP and WRKY were the most abundant TF families in the leaves network with 7 and 5 members, respectively (Table S11). Notably, the most abundant TF families in the regulatory network of leaves were not shared with the root network (Table S11), suggesting organ-specific regulation of the transcriptional response to sulfate starvation.
We next focused our analysis on TFs whose predicted target genes have overrepresented biological functions (adjusted p-value <0.05), in order to correlate TF regulation to changes in functional processes. As shown in Fig. 6A and Fig. 6B, we found 3 TFs in roots and 10 TFs in leaves that potentially control genes involved in the sulfate starvation response, phytohormone biosynthesis and signaling, senescence, among other biological processes. Most connected factors in terms of outdegree were a MYB TF (Solyc11g071300) in roots and TCP 21 TF (Solyc03g006800) in leaves, showing 40 and 71 targets respectively (Fig. 6A and Fig. 6B). Notably, MYB target genes are involved in S utilization, while TCP21 ones are involved in photosynthesis (Fig. 6A and Fig. 6B, Table S13).
As mentioned above, our results discovered important changes in gene expression in response to sulfate starvation, suggesting the involvement of different sets of TFs acting in both roots and leaves. However, we found one TF shared by both root and leaf gene regulatory networks which is early induced in roots but has a later response in leaves (Fig. 6C). This TF is member of the EIL family and corresponds to ETHYLENE INSENSITIVE 3-LIKE 3 (EIL3, Solyc01g006650). The closest Arabidopsis gene in terms of amino acid identity is SLIM1, which is a key regulatory factor of the sulfate starvation response . In addition, GO enrichment analysis showed that the identified tomato EIL3 TF might control the expression of a set of genes related to sulfate starvation, including classical sulfate-responsive genes such as LSUs as well as genes encoding enzymes involved in the reductive steps of sulfate assimilation (Table S14). Furthermore, we found that 60% of the tomato TF target genes are shared with targets predicted for Arabidopsis SLIM1, according to PlantRegMap, including sulfate-response genes such as LSUs and SDIs (Table S14). Overall, we identified several candidate TFs to mediate the transcriptome reprogramming that occurs during sulfate starvation response in roots and leaves, highlighting the possible key regulatory role of the EIL3 TF in tomato.