The experimental platform schematically illustrated in Supplemental Figure 1 , was designed to compare three Arabidopsis genotypes, WT and two mutants that are blocked in autophagy (i.e., atg7 and atg9), all in the Col-0 ecotype background. These genetic stocks were grown in normal conditions (i.e., (+N)-condition) and in conditions that were designed to induce autophagy, by depriving nitrogen from the growth medium (i.e., (-N)-condition). Three days after this environmental induction of autophagy, seedlings were harvested in triplicate for lipidomics and transcriptomics analyses. The harvested seedlings were separated into the root and aerial tissues, which were separately analyzed in order to evaluate if autophagy is distinct between these two organs. Initially a small segment of the roots was stained with monodansylcadaverine and examined by fluorescence microscopy to visualize autophagosomes (Contento et al., 2005), which confirmed that the autophagy process was induced in the WT plants by the (-N)-treatment, whereas this process was inhibited in the atg7 and atg9 mutants (Hanaoka et al., 2002; Thompson et al., 2005). Regardless of the seedling genotype, all seedlings showed an obvious yellowing in the autophagy-inducing condition, which began to appear approximately 48-hours after the (-N) treatment.

Sequencing of the RNA and statistical analysis of the read-counts from each of the resulting tissue samples identified the differentially expressed genes (DEGs) whose expression was affected by the genetic and/or environmental alteration in the autophagic state of the seedlings. The statistical significance of these comparisons was evaluated at two levels of significance, q<0.05 and q<0.1 ( Table 1 ; Supplemental Table 1 ). The entire list of DEGs identified by these comparisons are provided in Supplemental Table 2 . Because the evaluations using a statistical significance level of q<0.1 generated more comprehensive lists of DEG that ultimately provided higher confidence in the GO enrichment analysis, all subsequent evaluations of DEGs were conducted using this latter statistical significance level.

These comparisons indicate that the environmental condition that stimulated autophagy (i.e., (-N)-treatment) induced major changes in gene expression among all genotypes and organs. The atg7 mutant showed the largest number of DEGs in the (-N)/(+N) comparison (~8,300 DEGs in leaves and ~5,200 DEGs in roots), as compared to the numbers observed in the WT genotype (~6,700 DEGs in leaves and ~3,400 DEGs in roots). In contrast, the atg9 mutant showed the lowest number of DEGs in this comparison (~3,200 DEGs in leaves and ~3,700 in DEGs in roots). In all these comparisons, the number of up-regulated DEGs were 10-40% greater than the number of down-regulated DEGs.

The conclusions based on the (-N)/(+N) induced changes in the transcriptomes are reinforced when we compared the number of DEGs observed in the atg7 vs. WT and atg9 vs. WT comparisons in the (-N) condition. Namely, the atg7 mutation has the larger effect in altering the transcriptome (i.e., ~5000 and ~1200 DEGs in leaves and roots, respectively), than the atg9 mutation (i.e., 350 and<100 DEGs in leaves and roots, respectively). These rather large differences in the altered transcriptomes between the two mutants and WT, which are induced by the (-N) treatment, contrast to the situation we found when these comparisons were made in the (+N) condition. Specifically, in the (+N) condition there were only 24 and ~250 DEGs in leaves and roots, respectively, between atg7 and WT plants, and for the atg9 mutant these respective numbers were 34 and 58 DEGs. In most of these comparisons (with the exception of the atg9 mutant), the leaf tissue is considerably more responsive than the root tissue to changes in the (-N)-treatment and/or genotype.

There is very small overlap among the lists of DEGs that respond to the (-N)-treatment between the atg7 mutant and the atg9 mutant. Namely, of the ~5000 DEGs that are detected in the atg7 (-N) vs WT (-N) leaf comparison, only 312 also occur in the atg9 (-N) vs WT (-N) leaf comparison ( Figure 1A ). The same conclusion is drawn from the comparison of root transcriptomes; namely, of the ~1200 DEGs in the atg7 (-N) vs WT (-N) comparison, only 43 also occurred in the atg9 (-N) vs. WT (-N) root comparison ( Figure 1A ). However, when we considered the smaller number of DEGs that are detected in the atg9 (-N) vs WT (-N) comparison, a large majority of these (between 50% and 90%) are also DEGs in the atg7 (-N) vs WT (-N) comparisons of both roots and leaves ( Figure 2A ). It should be noted that a small number of these latter DEGs (<10) switched from being induced or suppressed in one mutant to being suppressed or induced, respectively, when one compares their expression in the other mutant.

Functional insights from these DEG lists were gained by conducting GO enrichment analysis via the sequential use of the annotation pipelines, DAVID (Huang et al., 2009; d. W. Huang and Lempicki, 2009) and REVIGO (Supek et al., 2011). These analyses identified 177 GO enrichment terms that are >5-fold enriched among the down-regulated and upregulated DEG lists. These terms were further categorized as being either common between the two tissues (i.e., roots and leaves) or they are uniquely associated with either the root or leaf tissue. Supplemental Figure 2 identifies the 15 GO terms that are common to both leaves and roots, and these include “secondary metabolic process” and “autophagy” terms in atg7 vs. WT comparison, and the “lipid binding” term in atg9 vs. WT comparison. Figures 1B , 2B show the GO terms that differed among the two tissues (i.e., leaves and roots). These GO terms were further categorized into five higher level functional annotations, i.e., metabolism, DNA replication/cell cycle, signaling, photosynthesis or lipid metabolism, and these higher order functionalities are color coded in Figures 1B , 2B . These analyses indicate that the atg7 and atg9 mutations affect different down-stream processes, and in each mutant the processes that are altered in leaves are different from those that are altered in roots. Specifically, in the atg7 mutant the upregulated DEGs in leaves are enriched with GO terms associated with metabolism and signaling functionalities, whereas the down-regulated DEGs in this organ are associated with DNA replication/cell cycle functionalities ( Figure 1B ). The GO terms that are categorized in the signaling functionality identify processes that are specifically associated with autophagy (e.g., pre-autophagosomal structure and autophagosome) and with a combination of stresses that are possibly mediated by changes in the oxidative status of the root tissue (e.g., hypoxia and anoxia, glutathione and ascorbate metabolism, and responses associated with cell death, senescence, stress granules, and herbicides). Similar analyses of the GO terms identified by the DEGs in the atg9 mutant indicate that metabolism functionalities are affected in the leaves, whereas in the roots the predominant functionality that is affected is signaling ( Figure 2B ).

By combining the comparisons of the transcriptomes of the atg7 and atg9 mutants to those of WT plants that were grown in either the normal, (+N)-condition or the autophagy-inducing, (-N)-condition, we were able to distinguish transcriptomic changes that occurred because of the change in the N-status of the plants from those that are associated with the induction of autophagy. This identification depended on the fact that the atg7 or atg9 mutations block autophagy, and therefore we reasoned that the changes in the transcriptome that occurred in response to the (-N) or (+N) treatment of either the atg7 or atg9 mutants that were different from those observed for the WT plants are associated with the process of autophagy. This strategy is schematically illustrated in the expression ratio plots shown in Supplemental Figure 3A . Specifically, the plot illustrates the expression ratio of a hypothetical gene in the (-N)/(+N) comparison in WT plants (green datapoints), and in autophagy mutant plants (orange datapoints). Genes that show an expression ratio where the green and orange datapoints are indistinguishable from each other, but show altered expression due to the (-N) treatment, were identified as responding to the (-N)-treatment of the plants but not altered by autophagy ( Supplemental Table 4C ). All other combinations of gene expression ratios, where the green and orange datapoints are statistically distinguishable from each other (q<0.1), are changing their expression levels in response to the autophagic state of the plants. These autophagy-responding genes can be further classified into three categories: a) those genes that showed enhanced differential expression in the autophagy mutant compared to the WT; b) those genes that showed reduced differential expression in the autophagy mutant compared to the WT; and c) those genes that showed opposite differential expression in the autophagy mutant than in the WT.

Using this integrated comparison strategy, we identified 863 DEGs in leaves and 239 DEGs in roots whose expression is responding to the atg7-defect in autophagy. In parallel we identified 14 DEGs in leaves and 8 DEGs in roots whose expression is responding to the atg9-defect in autophagy ( Figure 3 ; Supplemental Figure 3B ). Thus, by integrating these comparisons we were able to distinguish the smaller subset of autophagy-responding DEGs from the larger number of DEGs that are responding to the change in the N-status of the plants ( Table 1 ). The functionalities of the genes that are responding to the (-N)-treatment ( Supplemental Figures 4 – 6 ) or to the autophagy defects ( Figure 3 ) were explored by conducting GO enrichment analysis.

Consistent with prior studies, which indicate that the atg7 mutant has a more severe phenotypic defect than the atg9 mutant (Shin et al., 2014), we found a larger number of DEGs that are responding to the atg7-defect in autophagy than the number responding to the atg9-defect in autophagy. The Venn diagrams shown in Supplemental Figure 3B indicate that there is very little overlap among the autophagy-responsive genes in the atg7 and atg9 mutations, and similarly, there is a very small overlap in these autophagy-responsive genes between leaf and root tissues. These findings suggest that in the absence of either ATG7 or ATG9 gene-functions there are distinct down-stream effects, which may also be reflective of the differential expression pattern of these genes between leaves and roots.

Because of the small number of genes whose expression is affected by autophagy due to the atg9 mutation (22 DEGs), it was difficult to further assess functionalities that are associated with this genetic defect. However, such analyses were pursued with the gene list identified from the atg7 mutant ( Figure 3A ). Only a small fraction (~5%) of these autophagy-responsive DEGs occur in both leaves and roots, and the majority (~80%) are unique to leaves. Using the list of leaf and root DEGs as inputs we conducted GO term enrichment analysis to identify a broad overview of the types of processes that are affected by autophagy. These analyses identified 92 GO terms that show statistically significant enrichment (p< 0.1) ( Supplemental Figure 7 ), and Figure 3 focuses on those GO terms that show more than 5-fold enrichment. Comparison of the GO terms recovered with the list of leaf and root DEGs indicates that different processes are being affected by autophagy in the two organs. This difference is most pronounced with the GO Biological Process terms, which are populated by signaling processes in roots, whereas in leaves they are populated by metabolism associated processes. This enrichment of metabolism associated terms is also apparent in the GO Molecular Function category.

Lipidomics data were gathered in parallel to the described transcriptomics data. These analyses identified 156 specific lipid species, which included LysoPC, LysoPE, LysoPG, PC, PA, PG, PE, PI, PS, MGDG and DGDG. Upon statistical analyses, we obtained rigorous data on the abundance changes for 97 lipid molecules. PCA analyses of the leaf lipidomics data indicates that the major driver that separates the lipidomes of these samples is associated with the nitrogen treatment of the seedlings ( Figure 4A , Supplemental Figure 8 ). Namely, PC1 that accounts for ~54% of the variance in the data, separates samples based upon whether they were grown on (+N) or (-N) media. Separation of the sample by genotype is captured by PC2, and accounts for only 7% of the variance in the data.

Figure 4B illustrate that there are only a limited number of lipids species whose abundances are significantly affected by the atg7 mutation (3 and 14 specific lipids in the (+N) and (-N) conditions, respectively), and these appear to be enriched with lyso-PC, PC, PI and MGDG lipids. As with leaf tissue, the effect of the atg7 mutation on the root lipidome was more evident in the (-N) treatment, with the abundances of 13 specific lipids being altered in the mutant ( Figure 4B ), whereas in the (+N) condition, no lipid species where affected by the atg7 mutation ( Figure 4B ). Comparison of the lipidomes of roots and leaves indicate that the atg7 mutation induced different changes in the two organs; only ~15% of the altered lipids are shared between leaves and roots. In contrast to the effect of the atg7 mutation on the leaf and root lipidomes, the atg9 mutation did not affect a detectable change in the lipidomes of these organs, irrespective of whether the plants were grown in the nitrogen replete or nitrogen deficient condition ( Supplemental Figure 8 ).

We used MALDI-MSI to query the spatial distribution of specific lipids in the organs from WT and the autophagy mutants. Initial studies evaluated both roots and leaves of these plants grown either in the (+N) or (-N) conditions. Unfortunately, examination of root cross-sections by MSI proved technically difficult, probably due to the low abundance of lipids in this organ, and only a limited set of data could be collected ( Supplemental Figure 10 ). Therefore, a complete set of MS images were collected from the leaves of WT, and atg7 and atg9 mutant plants that were grown either in the (+N) or (-N) conditions. To ensure reproducibility of the data, all images were generated from at least three different biological samples, and Supplemental Figure 9 shows images of an exemplary set of triplicate data of the lipid, LysoPG(18:2). Similar reproducibility was obtained with all the lipid images, and Figure 5 shows representative images for 13 distinct phospholipids. These lipid distribution images were selected to ensure we determined the spatial distribution of a wide range of phospholipids, and warranted by a robust imaging signal, which ensured reliable reproducibility in the collected data. These data indicate that of the lipids that were imaged, there is no apparent differential distribution of each phospholipid molecule among the cell populations of these leaves, irrespective of the autophagic status of the seedling determined by the combination of genotype and nitrogen treatment.

Gene loci encoding functions associated with lipid metabolism in Arabidopsis are categorized in the AraLip database (http://aralip.plantbiology.msu.edu/) (Li-Beisson et al., 2013). AraLip identifies 775 Arabidopsis genes that have either been demonstrated to be involved in lipid metabolism or are hypothesized to be associated with such processes. We examined the transcriptome data to specifically identify those lipid genes whose expression were affected by either the atg7 or atg9 mutation in the context of the nitrogen treatment. These analyses identified 197 lipid genes (i.e., the altered lipid transcriptome) whose expression was significantly affected by either the atg mutations or the nitrogen treatment ( Figure 6 ; Supplemental Table 3 ). The altered lipid-related transcriptome is larger in the atg7 mutant than in the atg9 mutant, and in both mutants the size of this altered transcriptome is increased significantly by the (-N) treatment ( Figure 6 ); specifically, there are 184 lipid-related DEGs in the atg7 mutant, and only 13 such DEGs in the atg9 mutant. This parallels the changes we observed in the lipidome of the atg7 and atg9 mutant seedlings, with hardly any detectable changes in the lipids of the latter mutant ( Figure 4 ; Supplemental Figure 8 ). Moreover, the size of the altered lipid-related transcriptome is larger in the leaves than in the roots [i.e., more DEGs in leaves (146 genes) than roots (75 genes)] ( Figure 6 ). In the case of the atg7 mutant grown in the (-N) condition, where the altered lipid-related transcriptome is sufficiently large to provide statistical rigor, the number of the DEGs that are upregulated is between 2- and 4-fold larger than the number of down regulated genes, in both leaves and roots.

Insights into the lipid-related gene functionalities that were affected by these treatments were gleaned by identifying the lipid pathways pinpointed by the DEGs ( Supplemental Table 3 ). These annotations identified 15 lipid metabolic pathways, which can be classified as either a) biosynthetic (i.e., fatty acid biosynthesis, elongation and metabolism, chloroplastic lipid metabolism, cuticle and wax esters deposition and triacylglycerol biosynthesis); b) catabolic (i.e., fatty acid degradation, triacylglycerol metabolism and lipases), c) membrane lipid reorganization (i.e., lipid transport and binding proteins, and phospholipid metabolism) and d) signaling lipids (i.e., phosphoinositide, oxylipin and sphingolipid metabolism). Because the affected lipid-related transcriptome is larger in the atg7 mutant than the atg9 mutant, we were able to surmise that in the leaves and roots of the former mutant, the down-regulated genes in the (-N) treatment are enriched with lipid biosynthetic pathways and pathways that reorganize membrane lipids In contrast, the up-regulated genes are enriched in processes associated with lipid catabolism and the metabolism of signaling lipids.

Seed stocks of wild-type Arabidopsis (ecotype Col-0) and atg7-2 (Chung et al., 2010) and atg9-4 (SALK_145980) (Floyd et al., 2015) mutant lines were obtained from the Arabidopsis Biological Resource Center (Columbus, OH). Seeds were sterilized with 70% ethanol, followed by a 10-min incubation in 0.1% (v/v) Tween 20 (Thermo Fisher Scientific, Waltham, MA) and 50% (v/v) bleach solution. The seeds were then washed with sterile water, at least three times. Subsequently, the suspended seeds were vernalized by incubating at 4°C in darkness for 2 days.

After vernalization, the seeds were suspended in sterile 0.1% (w/v) agarose (VWR, Radnor, PA), and sown on 3.5 cm x 4 cm autoclaved stainless steel growth mesh (14 Mesh T304 Woven Stainless 0.017” wire diameter, TWP Inc. Berkeley, California), which was laid on ½ strength Murashige and Skoog (MS) solid medium composed of 2.15 g/L Murashige and Skoog Basal Salt Mixture (MilliporeSigma, Burlington, MA), 0.05% (v/v) Murashige and Skoog Vitamin Solution (MilliporeSigma), 1% (w/v) sucrose (Thermo Fisher Scientific), 6g/L Phytoblend Agar (Caisson Labs, Smithfield, UT), and 2mM MES (MilliporeSigma) at pH 5.7. Each 10 cm x 10 cm square Petri dish, containing 4 growth meshes, was placed in a growth room maintained at 22 °C under continuous illumination (50 ± 10 μE m⁻² s⁻¹) for a period of 5 days.

Subsequently, the growth mesh, carrying the germinated seedlings, was sterilely moved onto a sterile 7.5 cm x 8.5 cm stainless steel platform mesh (10 Mesh Woven Stainless 0.025” wire diameter, TWP Inc.), which was in a 11.4 cm × 8.6 cm × 6.4 cm Phytatray dish (MilliporeSigma) that contained sterile liquid medium, composed of ½ strength MS liquid media, which contained 10 mM NH₄NO₃ and 9.4 mM KNO₃ (i.e., the (+N) condition). The volume of the medium was adjusted so that the growth mesh that carried the seeds was in contact with the surface of the medium, and thus as seedlings grew the root system extended into the liquid medium. After 1 day incubation in the +N liquid medium, half the growth meshes from each Phytatray dish were moved into a Phytatray dish that contained nitrogen-deficient liquid medium (-N medium, which contains no nitrogen salts). This medium was composed of 5% (v/v) Murashige and Skoog Basal Salt Micronutrient Solution (MilliporeSigma), 0.05% (v/v) Murashige and Skoog Vitamin Solution, 1.5 mM CaCl₂, 0.75 mM MgSO₄, 0.625 mM KH₂PO₄, 2.5 mM KCl, 2 mM MES and 1% (w/v) sucrose. After an additional 3-day incubation, the seedlings from both the +N and -N media were harvested by cutting the hypocotyls that were extending below the growth mesh, and leaf and root tissues were collected separately.

For MDC staining, Arabidopsis seedlings were stained with 50 mM monodansylcadaverine (MDC) (MilliporeSigma) in phosphate-buffered saline, pH 7.4 for 10 minutes, followed by 3 brief washes with phosphate-buffered saline, as previously described by (Contento et al., 2005). MDC-stained seedlings roots were imaged in the late elongation zone and neighboring cells in the differentiation zone, while root tips and older differentiation region cells were excluded. MDC fluorescence was visualized using a Zeiss Axio Imager.A2 upright microscope (Jean, Germany) using a 40X objective and a UV lamp and 4’,6-diamidino- 2-phenylindole-specific filter (DAPI). The number of motile MDC-stained puncta in all cells visible in the focal plane was quantified and expressed per frame for at least 10 images per sample for 3 biological replicates.

Excess moisture was removed from the collected fresh tissue by blotting the tissue with dry paper towels, then the tissue was transferred to 15 mL centrifuge tubes (Thermo Fisher Scientific), flash frozen with liquid N₂, and stored at -80°C. About 1 g of harvested plant tissue was pulverized under liquid N₂ in an RNase-free mortar and pestle. RNA was extracted with TRIzol Reagents (Thermo Fisher Scientific), followed by treatment with DNase I (QIAGEN, Hilden, Germany) and clean-up with RNeasy Mini Kit (QIAGEN). The quality of RNA was analyzed with an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA). The RNA integrity number (RIN) of all samples was ≥7.8.

cDNA library construction and sequencing were conducted using a BGISEQ-500 system with 100-cycle paired-ends at BGI Americas (https://www.bgi.com/us/). Data files have been deposited in the NCBI Sequence Read Archive (https://www.ncbi.nlm.nih.gov/sra) under accession number PRJNA769614. The annotation file in the gff3 format and the genome DNA sequence of TAIR10 Release-39 were downloaded from http://ftp.ensemblgenomes.org/pub/plants/release-39/gff3/arabidopsis_thaliana/ and http://ftp.ensemblgenomes.org/pub/plants/release-39/fasta/arabidopsis_thaliana/dna/respectively. The splicing sites are extracted from the annotation file. The reads are aligned using HISAT2 (hierarchical indexing for spliced alignment of transcripts) (V2.1.0) (Kim et al., 2015) with known splice sites, mapped reads are counted by HTSeq-count (V0.9.1) (Anders et al., 2015) for raw counts and by cufflinks (V2.2.1) (Trapnell et al., 2010), for normalized counts.

Lipids were extracted using a modification of a standard protocol (Shiva et al., 2018). The collected plant tissues (leaf or root) were transferred to a 50 mL Teflon-lined screw-capped glass tube (Thermo Fisher Scientific) containing 3 mL preheated isopropanol (Thermo Fisher Scientific) containing 0.01% (v/v) butylated hydroxytoluene (BHT) (MilliporeSigma) and 1 µM 1,2-didecanoyl-sn-glycero-3-phosphocholine (MilliporeSigma) as an internal standard. The tubes were incubated at 75°C for 15 min to quench the action of any lipases.Following the addition of 1.5 mL chloroform and 0.6 mL water the mixture was vigorously shaken at room temperature for 1 h. The clear liquid extract was transferred to another 50 mL tube using glass Pasteur pipettes, and the remnant tissue was further extracted for 30-minutes with another 4 mL chloroform/methanol (2:1) that contained 0.01% BHT. The clear liquid from this second extraction was removed and combined with the initial liquid extract. This chloroform/methanol (2:1) extraction was repeated three times, the last extraction being incubated overnight. The residue tissue remaining after lipid extraction was dried at 105°C, and the dry weight of each sample determined; each leaf tissue sample weighed approximately 20 mg, and each root tissue sampled weighed approximately 10 mg.

All extract aliquots from each biological sample were combined into a single screw capped tube and stored at -80°C under a nitrogen gas atmosphere. The solvent from each extract was removed by evaporation with the aid of a stream of N₂ gas, and the lipid residue was dissolved in 1 mL chloroform and transferred to 2.0 mL clear glass vial with Teflon-lined screw cap (Thermo Fisher Scientific). The solvent was again evaporated with N₂ gas, and the vials were shipped overnight on dry ice to the Kansas Lipidomics Research Center (https://www.k-state.edu/lipid/) for lipidomics analysis.

Mass spectrometric analysis of the lipidome was conducted at the Kansas Lipidomics Research Center using a Xevo TQ-S mass spectrometer (Waters Co., Milford, MA). Individual lipids were identified by direct infusion in positive ion mode with precursor and neutral loss scans (Xiao et al., 2010; Peters et al., 2010; Li et al., 2014), using the scans shown in Supplemental Table 4 . Response factor corrections were applied to the MGDG and DGDG analyses to correct for differences in the response of the mass spectrometer to unsaturated galactolipid species as compared to the saturated internal standards. Phospholipid data did not require such response factor corrections, as the biological phospholipids and the internal standard have similar response factors.

The metabolomics data from this study are available at the NIH Common Fund’s National Metabolomics Data Repository (NMDR) website, the Metabolomics Workbench (https://www.metabolomicsworkbench.org) (Sud et al., 2016), where it has been assigned Study ID ST002252. The data can be accessed directly via its Project DOI: http://dx.doi.org/10.21228/M8KQ6P.

Lipidomics data were gathered from 44 individual samples that represented three or four replicates per +N/-N treatment for each genotype that was evaluated. Collected data were normalized using the upper-quartile normalization method (Bullard et al., 2010). Transcriptomics data were collected from 41 tissue samples that were gathered in parallel to those used for lipidomics analysis. Due to potential contamination of three root samples by leaves– indicated by assessing the abundance of transcripts encoding for 15 genes coding for photosynthetic functions that would be expected to be expressed in leaves and not in roots [i.e., RUBISCO small subunit genes (At1g67090 and At5g38410), photosystem I subunit genes (e.g., AT1G30380 and AT1G52230)] - the root replicate samples (one of each genotype) were removed from further analysis. The remaining 38 samples that represented between three and four replicates of root and leaf tissues of the three genotypes, grown in either (+N) or (-N) conditions were further analyzed. Initially, transcripts which showed undetectable expression levels across all the samples were removed, which provided abundance data for 26093 transcripts in the root samples, and 26356 transcripts in the leaf samples. The transcriptomics data was analyzed by using the DESeq function in R package DESeq2 (Love et al., 2014). Data was normalized by using the median of ratio method (Anders and Huber, 2010) in DESeq2 package.

Three statistical tests were applied to the transcriptomics datasets:

Test 1) comparison between N-treatments for a given genotype, e.g.,

Test 2) comparison between genotypes for a given N-condition (+ or -), e.g.,

Test 3) N-treatment effect tests among the 3 genotypes, e.g.,

Pearson’s correlation was calculated between logcpm, which was defined by log₂[10⁶• (r_i,j + δ)/R_j ], where r_i,j is the expression data for transcript i, δ is the half of the smallest non-zero measurement/count, R_j is the upper quartile for sample j. The log fold change (LFC) and q-values (Nettleton et al., 2006) which were calculated from p-values were used for interpretation of the results.

Principal component analysis (PCA) was conducted using the pcomp function of the factoextra package in R. Volcano plots were generated by plotting the negative logarithm of the p value on the y-axis and the LFC on the x-axis in Excel. Heat maps of gene expression levels were created with Perseus (https://maxquant.net/perseus/) (Tyanova et al., 2016).

Gene lists were subjected to analysis via DAVID (https://david.ncifcrf.gov/home.jsp) (Huang et al., 2009; d. W. Huang and Lempicki, 2009). Function annotation chart was generated by DAVID for each gene list, and terms for Biological Process, Cellular Component and Molecular Function were further reduced using REVIGO (http://revigo.irb.hr/), choosing “Small (0.5)” as the option for reduction (Supek et al., 2011). The fold-enrichment and p-value from DAVID is matched to the reduced lists and plotted.

Leaves were prepared for mass spectrometric imaging as previously described (Klein et al., 2015). The first pair of true leaves were collected and the abaxial side was adhered to Scotch box sealing tape (3M, Maplewood, MN). After 1-2 h drying in a vacuum, the packing tape was folded so that the adaxial surface of the leaves adhered to the tape, and the tape-encased leaf-tape “sandwich” was passed through a rolling mill to ensure adherence of the tape to the surfaces of the leaf. The two adhering tapes were separated, which resulted in fracturing of the leaf through the mesophyll cell interior, resulting in two separated pieces of packing tape, each retaining the abaxial and adaxial halves of the leaf. Each of these pieces of packing tape, with the exposed leaf surface facing upward, were adhered to a microscope slide using double-sided tape. After acquiring optical images, the samples were subjected to matrix deposition by sublimating (Hankin et al., 2007) 2,5-dihydroxybenzoic acid (DHB; 98%; MilliporeSigma) or 1,5-diaminonaphthalene (DAN; 97%; MilliporeSigma) at 140°C for 4 min at a pressure of~50 mtorr.

Mass spectrometry imaging data were collected using a MALDI-linear ion trap-Orbitrap mass spectrometer (MALDI-LTQ-Orbitrap Discovery; Thermo Fisher Scientific). The instrument was modified to incorporate an external 355 nm frequency tripled Nd : YAG laser (UVFQ; Elforlight, Daventry, United Kingdom) and a f=60 mm focus lens (Feenstra et al., 2017). TunePlus and XCalibur (Thermo Fisher Scientific) were used to define imaging parameters and to acquire data, respectively. Leaves were acquired using 100 μm raster step size. Mass spectra were acquired with 10 laser shots per spectrum mostly using an Orbitrap mass analyzer (resolution of 30,000 at m/z 400) for an m/z scan range of 200–1000. Data were obtained in positive and negative ion mode, using DHB and DAN, respectively.

MS images were generated using ImageQuest (Thermo Fisher Scientific) with a mass window of ±5 ppm, and without normalization. The same minimum and maximum plot values were used across all images and compounds. Peak assignments were based on accurate mass measurement, and comparison with online databases (METLIN, https://metlin.scripps.edu). Images were saved from a screenshot of the 2-D image window.