Integrated Framework of the Immune-Defense Transcriptional Signatures in the Arabidopsis Shoot Apical Meristem

The growing tips of plants grow sterile; therefore, disease-free plants can be generated from them. How plants safeguard growing apices from pathogen infection is still a mystery. The shoot apical meristem (SAM) is one of the three stem cells niches that give rise to the above ground plant organs. This is very well explored; however, how signaling networks orchestrate immune responses against pathogen infections in the SAM remains unclear. To reconstruct a transcriptional framework of the differentially expressed genes (DEGs) pertaining to various SAM cellular populations, we acquired large-scale transcriptome datasets from the public repository Gene Expression Omnibus (GEO). We identify here distinct sets of genes for various SAM cellular populations that are enriched in immune functions, such as immune defense, pathogen infection, biotic stress, and response to salicylic acid and jasmonic acid and their biosynthetic pathways in the SAM. We further linked those immune genes to their respective proteins and identify interactions among them by mapping a transcriptome-guided SAM-interactome. Furthermore, we compared stem-cells regulated transcriptome with innate immune responses in plants showing transcriptional separation among their DEGs in Arabidopsis. Besides unleashing a repertoire of immune-related genes in the SAM, our analysis provides a SAM-interactome that will help the community in designing functional experiments to study the specific defense dynamics of the SAM-cellular populations. Moreover, our study promotes the essence of large-scale omics data re-analysis, allowing a fresh look at the SAM-cellular transcriptome repurposing data-sets for new questions.


Introduction
Unlike animals, plants are sessile and cannot defend themselves from potential danger through the mechanism of fight and flight. Being directly exposed to a biotic and abiotic environment, they perceive plethora of signals that impact their physiological wellbeing. They often re-arrange their internal networks in response to external cues in order to adjust their developmental priorities. Moreover, they regularly replace body organs and tissues, such as leaves, flowers, and cambium, as a token of commitment to their overall Bauplan (body plan) [1]. In plants, the pluripotent stem cells constantly The immune-related functions of the SAM are relatively less characterized as compared to other plant parts. The main objective of this study was to collect a vast array of transcriptomes that deal with various SAM cellular populations. We further aimed at subjecting them to enrichment analysis for immune functions. Another objective of this study was to investigate the transcriptional basis of CLV3p-triggered immunity in the SAM and compare it with flg22-based typical innate immune responses in plants. We further aimed at mapping a transcriptome guided SAM interactome to better comprehend immune-related modules in the SAM cellular populations. The overall purpose of this study was hence to provide a systems biology framework of the immune-related genes and their regulation in 10 different SAM cellular populations. Emphasizing the essence of large-scale data re-usability, our analysis offers repertoire of genes that can be functionally dissected to unearth the molecular basis of defense mechanisms in the SAM.

Acquisition of the Arabidopsis Transcriptome Datasets Pertaining to SAM Cellular Populations and PAMPs-Treated Mesophyll Cells
The Arabidopsis SAM cellular populations are comparatively well studied for aspects of plant growth and development. However, immune responses of the various SAM cellular populations are not well studied. To reconstruct a transcriptional framework of the SAM cellular populations, we acquired 40 sets (raw data) of gene expression profiles from Gene Expression Omnibus (GEO) pertaining to various SAM cell populations, such as CLV3p-expressing stem cells belonging to the CZ, FILP-expressing cells belonging to the organ primordia/PZ, WUS-expressing cells belonging to the RZ, AtML1-expressing cell belonging to the epidermal cell layer, and AtHGD4 cells belonging to the sub-epidermal cell layer, HMG-expressing cells of epidermal cell types; likewise, KAN-and LAS-expressing cells belong to the peripheral zone cell types (see Table 1 for details). Moreover, the transcriptome of the cells belongs to vasculature cell types, such as shoot xylem cells/AtHB8 gene expressing cells and shoot phloem cells/ S17-gene-expressing cells, are also included in the analyses. To get a specific insight into the transcriptional regulation of stem cell populations (the CLV3p-expressing cells), we also included a non-functional (mutant) CLV3n-expressing cells to be compared with functional (wild type) CLV3p-expressing stem cells (Table 1, Table S1). Furthermore, we included a flg22-treated mesophyll cell gene expression profile in our study in order to find potential regulatory overlap and contrast it with CLV3p-mediated transcriptional signatures in the SAM [5]. The transcriptional profiles of all these cell populations (Table 1) were based on cell sorting and protoplasting procedures. To eliminate the effects protoplasting (treatment of the cells with cell wall degrading enzymes [12]) have on the transcription of the SAM-cellular populations and the flg22-treated mesophyll cells, we further included microarray datasets that deal with gene expression changes in cellular populations upon the application of protoplasting cocktails (Table 1).

Normalization and Analysis of the SAM Cellular Transcriptomes and Functional Categorization of the DEGs (Differentially Expressed Genes)
To maximally reduce batch and source effects in the analysis of the acquired SAM-based transcriptomes, we re-normalized all datasets of all 3 SAM populations (CLV3, FIL, WUS) acquired from the GEO accession GSE13596 and 7 SAM-cellular populations under the accession GSE28109. Each set was normalized and further analyzed for differential gene expression (DEGs) with R package limma. Taking each individual cellular population as test and the rest of the SAM-cellular populations as background, we identified DEGs in all 10 SAM-cellular populations (Table S1). After excluding 301 protoplasting induced genes, and based on the criteria of adjusted p-values and fold changes, we identified DEGs for SAM cell populations as follows; 192 for CLV3p-expressing cells, 43 for WUS-expressing cells, 143 for SILP-expressing cells, 85 for at ML1-expressing cells, 48 for LAS-expressing cells, 120 genes for KAN1-expressing cells, 188 genes for HMG-expressing cells, and 128 genes for HDG4 cells (Table S1). It is noteworthy to mention that the number of DEGs of AtHB8 and S17 cellular populations are manifold more than the rest of the SAM-cells. For instance, the DEGs of S17-expressing cells are 2824 and that of At HB8-expressing vasculature cells are 1862. Except for the vascular cells, most of the SAM cellular populations have distinct gene expression prolife, with minimal overlap among DEGs of the SAM cells (Table S1).  [12] In order to document the quantitative composition of all 10 SAM-cellular transcriptomes with a focus on the function of each DEG, we used Voronoi Tree-maps to visualize cellular gene expression datasets in conjunction with hierarchical classifications based on Kyoto Encyclopedia of genes and Genomes (KEGG) ortholog. Voronoi Tree tessellations group functionally relate genes whose expression patterns are highly correlated across the different experimental conditions. As shown in Figure 1, each functional category is shown by a color-coded polygon. Broadly, the DEGs of the SAM cellular populations participate in processes, such as cellular metabolism, transcription, replication, and DNA repair mechanisms ( Figure 1). Likewise, signal transduction and cell mobility, as well as immune defense-related functions are overrepresented in our analysis. Owing to the higher number of DEGs of the vascular SAM-cellular populations, the number of polygons is several folds greater in these cells in comparison to other SAM cell populations ( Figure 1). Although the plotted Voronoi Tree-map represents many different functional categories, the scope of this study was limited to cataloguing immune-related transcriptional signatures in the SAM. According to our analysis, the DEGs of various SAM-cell populations, such as HMG, LAS, CLV3p, and FILP-expressing cells, participate in plant-pathogen interaction pathways. Likewise, the DEGs of vasculature SAM cell populations are also shown to have a role in plant-pathogen interaction functions ( Figure 2

GO-Based Enrichment Analysis of the SAM-Cellular Transcriptome for Immune-Related Functions
The Voronoi Tree-maps can efficiently visualize transcriptomes with respect to their regulated pathways, and we already found that several of the SAM cellular populations DEGs ( Figure 2) participate in plant-pathogen interaction functions ( Figure 2). However, we further adopted an inclusive approach to focus on SAM-based immune genes and subjected the DEGs of individual cellular populations to GO-based enrichment analysis. Based on a previous study for the SAM cellular Arabidopsis transcriptome, where shoot and root derived microarrays were assessed for relative closeness in term of their transcriptome profile [12], we categorized SAM cell populations into four clades (sectors); LAS-, HGD4-, and CLV3p-expressing cells as Sector I; AtML-and HMG-expressing cells as Sector II; KAN-, FIL-, and WUS-expressing SAM cells as Sector III; and S17-and HB8-expressing cells as Sector IV in the SAM region. The total list of DEGs for each class is investigated for enriched biological processes using Fisher's exact test, no correction, and p-value limit of 0.05. The keywords "immune defense, pathogen infection, biotic stress, salicylic acid, and jasmonic acid" were used to collect immune-related GO terms for SAM DEGs.
According to our analysis (all immune-related DEGs passing the log2FC threshold in Table S4), 19 DEGs (Table 2) of the CLV3p-expressing (stem) cells participate in various immune functions, such as regulation of defense response by Exocyst subunit Exo70 family protein, Cathepsin B-like protease 3, Beta-D-glucopyranosyl abscisate beta-glucosidase, and Tryptophan aminotransferase-related protein 2. Likewise, Lipoxygenase 3, chloroplastic is a DEG of the stem cells that participates in jasmonic acid biosynthesis. Moreover, DEGs of CLV3p-expressing SAM cells that participate in innate immune responses are E3 ubiquitin-protein ligase ATL31 and CLAVATA 3,as well as NIMIN1. Similarly, the E3-Ubiquitin ligase PUB23 participates in immune effector functions and is expressed by CLV3p-expressing cell. The HGD4 cellular populations express nine DEGs that participate in immune functions, such as defense response to bacteria and fungi, innate immune functions, and jasmonate biosynthesis, as well as defense response to other organisms ( Table 2, details in Table S4). Similarly, three of the LAS cells DEGs are enriched in GO terms, such as defense to fungi and defense to bacteria, as well as the jasmonate biosynthesis. We further analyzed the DEGs that belong to SAM-AtML1-and HMG-expressing cellular populations for immune-related process ( Table 3, details in Table S4) and found that five of the AtML1 cells express genes, such as Protein PAT1 homolog, Duplicated homeodomain-like superfamily protein, Syntaxin-43, TIFY domain/Divergent CCT motif family protein, and Porphobilinogen deaminase, participate in immune-related functions. These DEGs are enriched in GO terms, such as innate immune response, response to external biotic stimulus, and regulation of defense response, as well as response to jasmonic acid (Table 3, details in Table S4). Likewise, the SAM-HMG-expressing cellular population expresses 25 genes which are enriched in immune functions, such as response to external biotic stimulus, response to jasmonic acid, regulation of defense response, defense response to bacterium, regulation of innate immune response, defense response to fungi, and defense response to other organisms (Table 3, details in Table S4). Prominent genes regulating these processes are: Chitinase family protein, Linoleate 9S-lipoxygenase 1, Receptor-like protein 12, WRKY transcription factor 17, Leucine-rich repeat (LRR) family protein, Pathogenesis-related thaumatin superfamily protein, Enhanced Disease Resistance 4 M(EDS4), Respiratory burst oxidase homolog protein D, DCD (Development and Cell Death) domain protein, and so forth ( Table 3). The third clade of the SAM cellular populations is comprised of FIL, KAN1, and WUS expressing cells. As per our analysis, the FIL-expressing SAM cellular population express 23 DEGs that participate in diverse immune processes, such as regulation of immune system processes, defense response to bacteria, and fungi as well as other organisms (Table 4, details in Table S4). Likewise, response to jasmonic acids and jasmonic acid biosynthetic process are among the overrepresented functions. NIMIN1, 4-coumarate-CoA ligase-like 5, Glutaredoxin-C9, Pathogenesis-related thaumatin superfamily protein, pectinesterase inhibitor 12, Serine/threonine-protein kinase OXI1, and WRKY transcription factor 48, as well as NPR5 and NPR6, are among the genes that may regulate immune response in FIL-expressing cells in the SAM of Arabidopsis (Table 4, details in Table S4). The KAN1-expressing cellular population regulate the expression of 15 genes that also participate in immune functions. Pathogenesis-related thaumatin superfamily protein, Defensin-like protein 13, CWINV1, Lipoxygenase 2, chloroplastic, Beta carbonic anhydrase 1, chloroplastic, L-type lectin-domain containing receptor kinase IV.3, and Powdery mildew resistance protein containing protein are DEGs that regulate defense to fungi and defense to bacteria, as well as response to external biotic stimuli ( Table 3, Table S4). As per our analysis, the WUS-expressing cells of the SAM only regulate 5 genes that participate in immune-related processes, such as regulation of the immune effector process, defense response to biotic stimuli, and defense response to fungi (Table 4, details in Table S4). The fourth clade of the SAM cellular populations is comprised of S17-and HB8-expressing cells populations. Being in a more differentiated zone (vasculature), they own many times more DEGs than their counterpart in the central zone. They have many genes that participate in immune functions (Table S2). Overall, there are many genes in the SAM cellular populations that participate in immune functions and are differentially regulated without the SAM being exposed to pathogen infections.

Transcriptome Guided SAM Interactome and Insights Into the Proteins of Immune-Related DEGs of the SAM-Cellular Populations
As mentioned above, the SAM is the niche of stem cells population (CLV3p-expressing cells) in the central zone that gives rise to cells of the peripheral zone which turn into differentiated cells of organ primordia. The cells of the various SAM layers (Reference [11,12]: Table 1) share boundaries through transmembrane proteins or cells of one layer may diffuse into another layer/s. We mapped an interactome (protein-protein interactions: acquired from Arabidopsis Interactome Database) by taking the DEGs of all 10-SAM cellular populations into account. Cellular interactomes are usually generic; however, proteome or transcriptome datasets can be used to define a specific context in a general interactome. We visualized the Arabidopsis Interactome by using the system biology tool "cytoscape" and mapped SAM-transcriptomes ( Figure S1) onto the generic Arabidopsis Interactome (Figure 2A: left panel). We then extracted the transcriptome guided SAM interactome out of the generic interactome and further mapped it into four sub-sectoral networks (Figure 2A: right panel). Each sub-network (sectors) represents cellular populations based upon their closeness [12] in terms of their transcriptional behavior (Tables 2-4; Sector I: LAS-, HGD4-, and CLV3p-expressing cells; Sector II: AtML-and HMG-expressing cells; Sector III: KAN-, FIL-, and WUS-expressing SAM cell; and Sector IV: S17-and HB8-expressing cells). We analyzed the SAM-Interactome to get further insight into the network behavior of immune-enriched genes among SAM-cellular populations.
As per our network analysis, Sector IV has high intra-sectoral interactions, which is contrary to the intra-sectoral behaviors of the nodes of Sector I, II, and III. However, the nodes of all fours sectors are engaged in high intersectoral interactions (Figure 2A: right panel). It is noteworthy to mention that the first three sectors (Sector I, II, and III) are engaged to a higher degree of interaction with nodes of Sector IV in comparison to the interactions within the first three sectors. We further analyzed the network behavior of proteins, the genes of which are enriched in immune functions (Tables 2-4) and are also differentially expressed in the SAM cellular populations (Table S1). We used the cytoscape plugin network-analyzer and found the degree of interaction for densely connected immune nodes (hubs) of the network ( Figure 2B). As per our analysis, only three nodes of Sector I have a degree more than or equal to two, Glycine-rich RNA-binding Protein 7 interacts with 5 other proteins in the network. Likewise, GRP8 and CCA1 have three and two interacting partners in SAM-interactome, respectively. The proteins of the immune genes of the Sector II are sparsely connected and only one of their proteins, Enhanced Disease Resistance 4 (EDS4), has just two interacting partners. Contrary to Sector II, the immune proteins of Sector III are densely connected. For instance, non-specific serine/threonine protein kinase has 9 interacting partners in the SAM-interactomes, whereas Glutaredoxin-C9 and WRKY 40 each have three interactors ( Figure 2B). The immune proteins of Sector IV are highly connected and have degree ranges from 6 up to 14 for some of the proteins. Besides degree of the immune-related proteins of the SAM-interactome, we also analyzed the topological behavior of the SAM-Interactome and found features of a typical regulatory network in our transcriptome guided SAM-Interactome network ( Figure S1).

Transcriptional Basis of the Arabidopsis CLV3p-Triggered Genes in the SAM and its Comparison with flg22-Mediated Gene Expression in the Mesophyll Cells
There is only one report in the literature that attributes the naturally immunized status of the SAM to CLV3p-FLS2 mediated immunity [8]. It is generally believed that FLS2, in its interaction, is only sensitive to flg22, which is an epitope of the bacterial protein flagellin and has been widely used in studies of innate immunity in plants. To compare the transcriptional basis of CLV3p-mediated transcriptional regulation in the SAM and flg22-regulated genes in a differentiated mesophyll cells, we retrieved transcriptome profiles of "CLV3p vs CLV3n" of the SAM-stem cell population and "flg22 vs Mock" of the Arabidopsis mesophyll cells. We split the DEGs regulated by both conditions into Up-Regulated Genes (URGs) and Down-Regulated Genes (DRGs). In order to identify a transcriptional overlap between these two conditions, we plotted the DEGs of both experiments in a Venn diagram (Figure 3). In our analysis, we did not find any reasonable overlap between the regulated genes of these distinct cellular populations. Four genes which are up-regulated by CLV3p are also up-regulated by flg22, whereas no overlap was observed among the down regulated genes of both conditions ( Figure 3). However, seven genes which are up-regulated by flg22 in mesophyll cells are shown to be down-regulated by the CLV3p in the SAM-stem cells population. Among the CLV3p and flg22 four mutually upregulated genes, two of them are cysteine/rich transmembrane domain proteins, while the other two belong to the dehalogenase and cytochrome p450 protein families (Table S3). We did not find mutually downregulated genes caused by the peptides flg22 and CLV3p. The genes which are reciprocally regulated in these two conditions belong to; protein families of sucrose transporters, LOGs, and E3-Ligase, genes related to the maintenance of the shoot apical meristem, toxin extrusion, and nitrogen metabolism (Table S3). Apparently, these genes have no direct involvement in the immune processes. We further analyzed the top 20 up and down differentially regulated genes by both flg22 and CLV3p in their respective cell types. Although there is no overlap among the top 20 up/down regulated genes mediated by both these proteinaceous ligands (flg22 and CLV3p), there is a sharp contrast in the magnitude of their regulation; the repression on gene expression that is caused by flg22 is quite mild in comparison to the stark reduction imposed by CLV3p on the SAM stem cell genes. Likewise, the induction of flg22 is higher than CLV3p in up-regulating genes in their respective cell types (Figure 3). The CLV3p up-regulated DEGs are mostly involved in SAM signaling networks, such as CLV3p, WUS, LOGs, and CKX. In contrast, the up-regulated DEGs of flg22 are typical immune-related genes, such as WRKY 17 and WRKY22, Peroxidase, lipoxygenase and Chitinase family proteins, Elicitors and PAMP-induced peptides, and so forth (Figure 3: lower panel). From these results, we conclude that, although CLV3p and flg22 (two different peptides in terms of sequence and structure) interact with the very same plant innate immunity receptor (FLS2), their transcriptional outputs are very different.
CLV3p up-regulated DEGs are mostly involved in SAM signaling networks, such as CLV3p, WUS, LOGs, and CKX. In contrast, the up-regulated DEGs of flg22 are typical immune-related genes, such as WRKY 17and WRKY22, Peroxidase, lipoxygenase and Chitinase family proteins, Elicitors and PAMP-induced peptides, and so forth (Figure 3: lower panel). From these results, we conclude that, although CLV3p and flg22 (two different peptides in terms of sequence and structure) interact with the very same plant innate immunity receptor (FLS2), their transcriptional outputs are very different.

Discussion
The sterile nature of plant apexes is exploited to produce disease-free plants, which is an old horticultural practice. The exact mechanisms why plant tips grow aseptic is still a terra incognita for the field of plant immunity. The Arabidopsis SAM stem cell niche is an aggregation of approximately 500 cells located at the growing tip of each shoot [9,11,12]. The CZ of the SAM harbors stem cells, the progeny of which are displaced into the adjacent PZ, where they proliferate before differentiating. Despite much progress in The column information deals with the gene IDs, logFC in term of heatmap, and the annotation of the given genes. The scale bar beneath the columns denotes the magnitude of both up-regulation (red) and down-regulation (green) for the given genes.

Discussion
The sterile nature of plant apexes is exploited to produce disease-free plants, which is an old horticultural practice. The exact mechanisms why plant tips grow aseptic is still a terra incognita for the field of plant immunity. The Arabidopsis SAM stem cell niche is an aggregation of approximately 500 cells located at the growing tip of each shoot [9,11,12]. The CZ of the SAM harbors stem cells, the progeny of which are displaced into the adjacent PZ, where they proliferate before differentiating. Despite much progress in understanding the dynamics of immune responses of plants in general, the defense mechanism of the SAM stem cell niche is hardly assessed [5][6][7][8]. The compact nature of cellular populations in the form of layers in the SAM makes it less amenable to perform typical plant-pathogen interaction experiments specifically in the niche where stem cells reside. However, the advent of single cell genomics, advanced microscopy, and sorting techniques, together with sequencing technologies, provide ease in assessing the SAM and RAM cellular populations for detailed analysis. Profiting from these technological advancements, the dynamics of growth and development regulated by stem cells activity in these niches are intensively investigated these days.
Giving momentum to the essence of large-scale data re-usability, as well as re-producibility, we applied a systems-based deconvolution approach to catalogue immune-related genes belonging to various SAM-cellular populations. We acquired raw data pertaining the transcriptomes of 10-SAM cellular populations [12] from the public repository GEO and assessed DEGs ( Figure S1) testing each population against all other cellular populations. Besides finding DEGs across various experimental conditions (cell types and cell treatments), we applied double enrichment analysis: Voronoi Tree-maps as general visualization of DEGs based on KEGG pathways ( Figure 1) and GO-based enrichment of the DEGs of the SAM cellular populations for immune functions (Tables 2-4). In a previous study, Yadav et al. [12] compared DEGs in terms of specific cell populations as group to see the impact of outliers in distinct cell layers (L1, L2, and L3). However, we are broadly interested in genes contained in the SAM cellular populations so as to document their immune functions and motivate experimentalists to further decipher immune functions applying functional genomics approaches.
SAM cellular populations follow a conservative strategy by regulating only a limited number of genes that address their functional requirements. Being committed to functions of meristematic activities, their specific gene expression profile is way too small in comparison to the DEGs in differentiated tissues; the latter express 1000s of genes under given circumstances. We found the response to fungi, response to bacterium, and response to jasmonic acid and jasmonate biosynthesis as statistically overrepresented immune functions in the SAM-cellular populations. CLV3p, HMG, and FILP-expressing cells are among the SAM cellular populations that express immune-related genes (Tables 2-4). Both CLV3p [5] and HMG [12] cells have previously been shown to express genes that are implicated in pathogen defense. Likewise, the cells of S17 and HB8 populations express several genes that contribute to immunity against pathogen infection in plants. As compared to other SAM cell types, we found ( Table S2) that 50% of the DEGs are co-regulated by both S17-and HB8-expressing cells of the SAM. Looking into the number of genes that are enriched in immune functions of the SAM, we conclude that, without being exposed to a pathogen, SAM already has a pre-meditated immunity, whereby various layers (cell populations) express immune genes as a preemptive measure in safeguarding the integrity of the SAM.
To get deeper insight into the function of immune-enriched genes, we linked SAM-DEGs to their encoded proteins. We contextualized a generic interactome AI (Arabidopsis Interactome) into SAM-interactome by mapping DEGs of all SAM-cellular populations onto the AI. The four defined sectors (Figure 2A) of the SAM-interactome have individual cellular populations that share similarities in their gene expression profile [12]. Nodes of the sub-networks have poor intra-sectoral interactions (edges), as well as sparse intersectoral interaction, among the nodes of Sector I, II, and III have been observed ( Figure 2). However, the vast majority of the cross-sectoral interactions happen between the nodes of Sector IV and first three sectors.
There exists a reasonable skepticism about the validity of high-throughput Y2-H based interactomes owing to false positive/negative interactions and technological discrepancies in detecting interactions between soluble proteins and insoluble ones in the cell. Despite this valid skepticism, the computational approaches applied on large cellular interactomes in detecting functional modules endow the experimentalists with a handle to dissect the biological complexity with more targeted approaches. A pertinent question concerning the mapped SAM-interactome ( Figure 2) is: How do interactions between proteins of different cellular populations interact within the SAM's sectors and layers? It remains to be established how proteins of various cellular populations and sectors can interact as a network inside the SAM. There are three plausible explanations to address this apprehension: the movement of the cellular populations inside the SAM; the movements of the proteins, such as WUS, CLV3p across the cells; and the occurrence of transmembrane protein domains at the cellular interfaces in the SAM. Under these circumstances, proteins of different cells may interact as a network in the SAM despite belonging to different cellular population. However, there might be other modes of protein-protein interactions in the SAM-cellular populations, and that margin of error in interactions (artefacts of methods that generate interactomes) cannot be ruled out completely. We identified hub nodes in the interactomes which participate in defense functions and have sectoral and cross-sectoral interactions in the SAM-Interactome. These hub nodes ( Figure 2) merit further investigation to fully unleash their role in SAM-mediated immunity in plants.
An interesting mechanism of bacterial cleansing in the SAM was demonstrated to involve the binding of CLV3p (stem cells expressed gene) to FLS2 and the subsequent activation of Mitogen-Activating Protein Kinase (MAPK) that culminates in the expression of PAMP-Triggered-Immunity (PTI) marker genes in Arabidopsis [5,17]. Contrary to flg22-FLS2 mediated immunity, which involves the inhibition of seedlings growth, the CLV3p-FLS2-triggered immunity operates without growth penalties [6]. This mechanism of stem-cell-triggered immunity was refuted by other reports, demonstrating that CLV3p is blind to FLS2 in eliciting immune responses and that CLV3p does not induce the expression of innate immunity genes via the receptor FLS2 [18,19]. Taking these reports into account, we compared a large scale-transcriptional response mediated by flg22 and CLV3p in their respective cellular populations (Figure 3). Our results clearly demonstrate that there is no reasonable overlap in the DEGs of CLV3p-expressing stem cells population and flg22-treated mesophyll cells. Despite binding the same receptor, the structurally different peptides flg22 and CLV3p regulate different sets of genes in Arabidopsis.
In conclusion, we deconvoluted a set of immune-related genes from large-scale transcriptome data belonging to different SAM-cell populations. The identified transcriptional framework will help the plant research community in dissecting mechanisms of defense in the SAM. Our systems biology analysis infers that not only a single cell population of the SAM is committed to immune functions; rather, many cell types, in their respective layers, express genes that potentially participate in immune processes in the SAM. Once the mechanism of SAM-mediated immunity is fully elucidated, this will provide unique opportunities to incorporate pathogen resistance conferring pathways in crop plants. One of the candidate genes that can create impact in a breeding program is CLV3p, the expression of which is more restricted to the stem cell zone of the SAM. Our systems biology analysis suggests that CLV3p-expressing cells also induce many defense-related genes, which are different from those induced by the flg22-triggered-FLS2-mediated pathway in plants. Whereas the activation of FLS2 based immune pathway has a growth penalty, the CLV3p mediated immunity does not inhibit plant growth and has already been tested in Arabidopsis [6,7]. Thus, the incorporation of such a mechanism in crop plants will have dual benefits: breeding for pathogen resistance and no growth arrest due to active defense responses. It is noteworthy to mention that the introduction of stem cell specific pathways in differentiated plant tissues might have pleiotropic consequences. Therefore, expressing such pathways under inducible promoters would offer a need-dependent expression of the transgene that might reduce the off-target effects of the transgene which sometimes prove more deleterious than the anticipated positive effects on yield and protection. In order to fully translate the benefits of stem cell triggered immunity in crop plants, it is direly required to first establish assays which should easily allow plant-pathogen interaction experiments in the SAM of various crop plants. However, the modern genome editing techniques, state-of-the-art microscopy and imaging analysis, as well the non-invasive plant phenotyping approaches, should allow us to capitalize on the translational benefits of SAM-triggered immunity in crop plants.

Source of Transcriptomes (Gene Expression Data/Profiles): Gene Expression Omnibus (GEO)
We acquired microarrays datasets with accession numbers GSE13596, GSE16472, and GSE28109 from GEO, which is a public depository of transcriptomics/genomics studies obtained with microarrays and sequence-based experiments. The database provides the data, along with the publication if available, and description of the study and samples (all metadata describing the experiments). The data can be directly downloaded and analyzed independently or with the help of tools provided by the database [23,24].
Specifically, GSE13596 and GSE28109 provide gene expression profiles of A.thaliana shoot apical stem cell populations [11] and shoot apical meristem [12], respectively; GSE16472 contains, instead, the gene expression data for the response of Arabidopsis mesophyll cells after bacterial flagellin (flg22) treatment. Raw data from Affymetrix Arabidopsis ATH1 Genome Array is downloaded as CEL files from aforementioned GEO studies. Information on samples in each data can be found in Table 1. The collected data is arranged into three different sets for further analysis. The first set is composed of 10 SAM populations with 3 populations (CLV3, FIL, WUS) being derived by the analysis of GSE13596 and 7 populations from GSE28109. The second set comprises the gene expression data of CLV3p-expressing and lacking stem cells in GSE13596 and the third set is built with the samples from GSE28109 comparing flg22-treated mesophyll cells in comparison to the mock.

R Packages
We normalized the acquired microarray datasets with Limma, which is a Bioconductor tool widely used for the analysis of gene expression data. It compares different conditions provided as experiment design and employs linear models to determine differential expression [25]. We also used Oligo, which is also a Bioconductor tool for the preprocessing and visualization of the oligonucleotide microarrays from Affymetrix and NimbleGen genes. It is composed of various algorithms, e.g., RMA, being commonly used as a tool for background correction and normalization of the microarray data [26][27][28][29].

Voronoi Tree-Maps
To visualize complex microarray datasets, we used Voronoi Tree-maps for the investigation of functional categories with the use of KEGG pathways and for visualization of functions prevailing in a given transcriptome with corresponding quantitative data. The resulting transcriptome map is a composition of polygons of different sizes and colors that represent the functional categories representing highly expressed genes in the submitted data [30].

The Gene Ontology Resource
To assess the analyzed DEGs for immune functions, we profited from the Gene Ontology resource, which is a comprehensive database of gene functions for diverse set of organisms. It provides "gene ontologies" annotated in terms of molecular function, biological process, and cellular component. Further analysis of gene ontologies, e.g., GO enrichment analysis, in a given gene list can be also performed using the database [31][32][33].

CCSB Interactome Database
The Interactome database from the Center for Cancer Systems Biology (CCSB) (http://interactome. dfci.harvard.edu/) includes interactomes from Homo sapiens, Arabidopsis thaliana, Saccharomyces cerevisiae, and Caenorhabditis elegans. We acquired protein-protein interaction datasets from CCSB and mapped transcriptomes datasets to Arabidopsis interactome in order to create a transcriptomes-guided SM-interactome.

Cytoscape 3.7.1 and NetworkAnalyzer
We used Cytoscape, which is a network visualization and analysis tool that allows the assessment of protein-protein interaction networks and signaling pathways in its core and further allows such networks for other analysis owing to several dedicated plugins it provides [34]. We used NetworkAnalyzer for analyzing SAM-interactomes in term of topological and functional analysis. It generates graphs of different topological parameters, such as node degree, number of neighbors, path length, shortest path length, and number of connected nodes [35].

Data Preprocessing, Analysis, and Filtering
Background correction and normalization of each set is performed separately using the RMA algorithm from the R package "oligo" [26][27][28][29]. A soft intensity-based filtering is applied on the normalized data: lowly expressed genes are filtered out by a threshold of 4 for the median intensities of genes as per the instructions of Klaus et al. [36]. Differential gene expression analysis is performed with R package limma. For each set, contrasts are created (control vs. test), and samples are annotated accordingly. After acquiring a linear model fit for each set, empirical Bayes statistics are implemented to get differential expressed genes (DEGs). The differential expressed genes are filtered against p-value < 0.05 and |log2FC| >1. For transcriptome sets, 301 genes reflecting the protoplasting effect are subtracted [11] from the SAM-cellular population gene expression profile.

Pathway Enrichment Analysis
Functional categories for differentially expressed genes in each SAM population is performed employing Proteomaps [30]. For each population, DEGs and the log2FC values are used to generate the Voronoi Tree-maps and third level functional categories that correspond to individual pathway maps in line with KEGG hierarchy levels [37], which were chosen for representation of the differentially expressed genes.

Gene Ontology Enrichment Analysis for Immune-Related Processes
For GO enrichment analysis, SAM populations are separated into two classes. Class I includes 8 SAM populations (FIL, WUS, CLV3, KAN, LAS1, HMG, HDG4, and ATML1: further divided into three categories based on closeness in their expression profile [12]) and class II includes 2 populations (ATHB8 and S17). The total list of DEGs for each class is investigated for enriched biological processes. Fisher's exact test, no correction, and p-value limit of 0.05 is used for the enrichment analysis. The complete list of enriched biological processes is then checked for the keywords "immune", "defense", "pathogen", "biotic", "salicylic acid", and "jasmonic acid". The terms resulted from the search are collected as immune-related GO terms, and the SAM DEGs that take part in them are listed for each GO term.

Network Construction, Visualization, and Analysis
The Arabidopsis thaliana interactome is downloaded from CCSB Plant Interactome Database (http: //interactome.dfci.harvard.edu/A_thaliana/index.php?page=download) and visualized in cytoscape 3.7.1 [34]. The DEGs from all 10 SAM populations are gathered together and then searched in the A.thaliana interactome. Selected nodes and their edges are used to construct the SAM