Regulation of Cell Type-Specific Immunity Networks in Arabidopsis Roots

Root cell types possess distinct immunity gene networks that are linked to cell identity networks, as revealed by cell type-specific RNA-seq and a paired motif enrichment tool for promoter analyses.


INTRODUCTION
Plant roots are essential for plant health and development. In addition to anchoring plants, roots take up nutrients and water and provide protection from soil-based microbes. Conducting these different tasks is especially challenging under changing environments and acute stress conditions. Roots have evolved complex tissues comprising a diversity of cell types with different functions. Organized in concentric layers, Arabidopsis (Arabidopsis thaliana) roots consist of an outermost epidermis, followed by cortex, endodermis, pericycle, and root vascular tissue containing xylem and phloem cells (Dolan et al., 1993;Brady et al., 2007). This organization is implemented by the stem cell niche in the very root tip where cell fate is determined, and cell types maintain their given identity throughout their lifetime (van den Berg et al., 1995;Sabatini et al., 2003;Wendrich et al., 2017). Cell types nevertheless possess some plasticity, as exemplified for founder cells that originate from xylem-pole pericycle cell files that initiate lateral root formation Du and Scheres, 2018). Recent developments in single-cell transcriptomics have helped to further characterize cell types and to define root development but still cannot provide deep transcriptomic profiles (Birnbaum, 2018;Denyer et al., 2019;Jean-Baptiste et al., 2019;Ryu et al., 2019;Zhang et al., 2019;Rich-Griffin et al., 2020). Studies of cell type-specific transcriptomics based on fluorescence-activated cell sorting (FACS), in turn, have significantly advanced our knowledge of the individuality of cell type function in regulating root integrity under changing environments (Birnbaum et al., 2003(Birnbaum et al., , 2005Dinneny et al., 2008;Gifford et al., 2008Gifford et al., , 2013Bargmann et al., 2013;Geng et al., 2013;Walker et al., 2017). While these studies revealed the importance of a coordinated regulation of cell type-specific gene networks to master root development and secure overall root functionality (e.g., growth) under abiotic stress, the function of root cell types in regulating root immunity remains elusive.
The surface of leaves and roots is the habitat of complex microbiomes consisting of 10 6 to 10 9 microbes (per cm 2 of leaf area or g of soil, respectively); among them are pathogens with lifestyles ranging from biotrophy (life cycle completed on living cells) to necrotrophy (killing cells to complete the life cycle; Berendsen et al., 2012;Bulgarelli et al., 2013;Lareen et al., 2016). Root diseases represent a major threat to crop production, and enhancing root resistance against pathogens by improving processes regulating pattern-triggered immunity (PTI) is an altogether untapped approach to sustain food security (Gewin, 2010;Alexandratos and Bruinsma, 2012;Delgado-Baquerizo et al., 2020;Panth et al., 2020). Plasma membrane-localized pattern recognition receptors (PRRs), to date primarily characterized in leaves, recognize microbe-associated molecular patterns (MAMPs) as nonself molecules from microbes to induce PTI (Jones and Dangl, 2006;Boller and Felix, 2009;Cook et al., 2015). Recognition of the archetypal MAMP flg22 (the active epitope of bacterial flagellin) by the PRR FLAGELLIN-SENSITIVE2 (FLS2) activates PTI responses, including the rapid production of reactive oxygen species (ROS), MITOGEN-ACTIVATED PROTEIN KINASE (MAPK) phosphorylation, and the induction of immunity genes to restrict pathogen infection Gómez-Gómez et al., 1999;Asai et al., 2002;Zipfel et al., 2004). Similarly, root cells can recognize flg22 via FLS2 to trigger effective PTI (Millet et al., 2010;Jacobs et al., 2011;Beck et al., 2014;Wyrsch et al., 2015;Poncini et al., 2017;Stringlis et al., 2018). In addition to MAMPs, plants produce damage-associated molecular patterns (DAMPs) in response to pathogens, which are recognized by PRRs as well. Pep1, one of the best studied DAMPs produced in Arabidopsis, is encoded by PROPEP1 and recognized by the plasma membranelocalized PEP RECEPTOR1 (PEPR1) and PEPR2, triggering similar PTI responses as flg22 (Huffaker et al., 2006;Krol et al., 2010;Yamaguchi et al., 2010;Flury et al., 2013). The PEPR1/2 and FLS2 pathways share common signaling components such as MAPKs (Schulze et al., 2010;Liu et al., 2013;Yamada et al., 2016) but retain certain key differences that might be at least partially explained by additional activities of PEPRs. Qi et al. (2010) identified a unique guanylyl cyclase activity for PEPR1 mediating apoplastic Ca 21 influx upon Pep recognition. PEPR1/2-triggered immune signaling was further shown to maintain PTI in plants impaired in MAMP perception and signaling (Tintor et al., 2013;Yamada et al., 2016). Thus, there are clear differences and interdependencies between flg22 and Pep1-induced PTI in Arabidopsis. Nevertheless, the gene networks and underlying regulatory patterns defining DAMP and MAMP-mediated PTI in roots are currently unknown.
Motivated by recent findings suggesting distinct competences of different root cell types in launching PTI (Wyrsch et al., 2015), we wanted to know if flg22 and Pep1 trigger different transcriptional networks in three Arabidopsis root cell types, epidermis, cortex, and pericycle, and if so, whether it would be possible to identify distinct cell type-specific regulatory patterns. Our study demonstrated that very distinct immunity gene networks are activated in the three cell types. Considering that homomeric or heteromeric tandems of transcription factors (TFs) are often sufficient to determine regulatory specificity in eukaryotic cells (Halfon et al., 2000;Vandepoele et al., 2006;Junion et al., 2012;Ezer et al., 2014), we conducted combinatorial TF binding motif analyses to explain the regulatory patterns of cell type-specific gene networks. More specifically, by developing a statistical test for the enrichment of paired TF motifs that accounted for a multiplicity of TF binding sites, we were able to explain cell type-specific differences of Pep1-and flg22-elicited immune networks by specific TF motif combinations. Moreover, our study suggested the importance of cell identity in determining cell type-specific immunity networks. We discuss the significance of such a regulatory connection in specifying cell type functionality and, thus, in securing root integrity under conditions of environmental stress.

RESULTS
flg22 and Pep1 Activate Root Immunity through Partially Nonoverlapping Signaling Pathways Treating Arabidopsis roots with the immunity elicitor flg22 or Pep1 induces PTI responses (e.g., ROS burst, MAPK phosphorylation, induction of PTI marker genes, and eventually inhibited plant growth; Supplemental Figures 1A to 1F). flg22 and Pep1 have been shown to act through overlapping pathways (Krol et al., 2010;Yamaguchi and Huffaker, 2011;Tintor et al., 2013). We previously demonstrated that the beneficial root endophyte Serendipita indica (formerly Piriformospora indica) suppresses PTI to facilitate root colonization and that flg22 treatment of roots inhibits S. indica colonization (Jacobs et al., 2011). In colonized Arabidopsis roots, this fungus inhibits MAPK phosphorylation, PTI marker gene induction, and growth inhibition after flg22 (Supplemental Figures  1B, 1C, and 1E) but not Pep1 treatment (Supplemental Figures 1A to 1C and 1E). Consistent with an effective Pep1-induced immunity, S. indica showed improved root colonization of the Pep1 receptor mutant pepr1 pepr2 (Supplemental Figure 1F). These data suggest that flg22 and Pep1 recruit, at least partially, different signaling pathways to activate PTI in roots.
To further explore if PTI can be activated across different root zones, we treated Arabidopsis lines (Poncini et al., 2017) expressing PTI marker gene promoters fused to nucleus-localized mVENUS (MYB DOMAIN PROTEIN51, pMYB51:NLS-3xmVENUS; PEROXIDASE5, pPER5:NLS-3xmVENUS) with flg22 and Pep1. The analyses were conducted to exclude gene induction by other stresses (e.g., wounding). In contrast to recent reports using the same PTI marker lines (Zhou et al., 2020), both elicitors induced all markers in the root apical meristem (RAM), transition zone (TZ), elongation zone, and differentiation zone of roots grown on ATS medium and to a much lesser degree on half-strength Murashige and Skoog (MS) medium (except for pPER5:NLS-3xmVENUS in RAM/TZ by Pep1; Figure 1A; Supplemental Figures 2A and 2B). Consistently, root growth inhibition was stronger in flg22or Pep1-treated plants grown on ATS medium (Supplemental Figures 2C and 2D). These findings indicate some extent of PTI suppression, likely because the MES-based buffer system commonly used in MS medium (but not in ATS medium) interfered with the well-known induction of pH changes in response to MAMP perception . As a result of this MS medium-based PTI quenching effect, all subsequent experiments were done with plants grown on ATS medium.

Root Cell Types Differ in Their Immunity Gene Networks
Considering the diverse functions of root cell types in root development (Birnbaum et al., 2003;Brady et al., 2007;Gifford et al., 2013) and abiotic stress signaling (Dinneny et al., 2008;Geng et al., 2013), we explored to what extent flg22 and Pep1 affected gene networks in different root cell types. For our studies, we used Arabidopsis lines specifically expressing GFP in epidermis (atrichoblast, pGL2:GFP), cortex (pCORTEX:GFP), or pericycle (xylem pole, E3754; Masucci et al., 1996;Brady et al., 2007;Gifford et al., 2008;Bargmann et al., 2013;Lin et al., 2015) and treated the roots of ;15,000 seedlings (per biological repeat) per line with either flg22, Pep1, or mock (Figure 1B;Supplemental Figures 3A and 3B). We selected these cell types due to the importance of the . In order to make consistent comparisons, the top 128 upregulated genes from each cell type were used to calculate GO enrichment. Tone of color indicates significance degree [log 10 (P adj )]. (E) and (F) The top 5 nonredundant, most significantly enriched GO terms in genes specifically induced by Pep1 in the epidermis (E) or cortex (F). In order to make consistent comparisons, the top 365 upregulated genes from each cell type were used to calculate GO enrichment. Tone of color indicates significance degree [log 10 (P adj )]. (G) and (H) The top 5 nonredundant, most significantly enriched GO terms in genes specifically suppressed by Pep1 in the epidermis (G) or cortex (H). In order to make consistent comparisons, the top 337 downregulated genes from each cell type were used to calculate GO enrichment. Tone of color indicates significance degree [log 10 (P adj )].
(I) Venn diagrams show the overlap of flg22and Pep1-responsive DEG sets aggregated across cell types (top) and the split by cell types for 194 flg22specific (bottom left) and 2392 Pep1-specific (bottom right) DEGs per cell type. epidermis and cortex (as outer, environment-facing cell layers) in protecting the root against pathogen invasion, while pericycle cells (the outermost cell layer of the inner root tissue and intimately associated with the vasculature) regulate lateral root formation and transport processes and are highly responsive to immune elicitors (Takano et al., 2002;Parizot et al., 2012;Wyrsch et al., 2015;Ross-Elliott et al., 2017). Importantly, previous studies have demonstrated PEPR1/2 and FLS2 function in these cell types and that flg22 and Pep1 reach pericycle cells within minutes after treatment (Wyrsch et al., 2015;Ortiz-Morea et al., 2016). flg22or Pep1-elicited roots were gently washed to remove flg22 and Pep1 before generating protoplasts from the roots. For each Arabidopsis line and treatment, ;20,000 GFPexpressing protoplasts were extracted using FACS ( Figure 1B; Supplemental Figures 3A and 3B). Importantly, the cell typespecific expression patterns of the marker genes were unchanged in immunity-activated roots ( Figure 1C), ensuring uniformity of our isolated cell populations. Similar to cell type transcriptome analyses of abiotic stress networks (Dinneny et al. , 2008;Geng et al., 2013), we treated whole roots rather than protoplasts to capture the tissue context of, and intercellular communication between, cell types. In addition, cells were analyzed at 2 h after elicitor treatment to capture early transcriptional changes defining effective PTI and to exclude gene network crosstalk resulting from growth inhibition as a later PTI response (Zipfel et al., 2004).
RNA-seq of FACS-isolated cells ( Figure 1B) resulted in ;315 million reads (;11.6 million read pairs per library) that were uniquely mapped to gene features in the Arabidopsis genome (Supplemental Table 1; Supplemental Figures 4A to 4F). Principal component (PC) analysis revealed that 82% of the variation was contained within the first three principal components, where PC1 (62% of the variation) separates between cell identity and PC2 (16% of variation) between treatments ( Figure 1D; Supplemental Figures 5A and 5B). Cell type marker genes were significantly expressed in the respective cell type populations, indicating that FACS was efficient for cell type isolation (Supplemental Figures 5C  to 5E).
First, we identified differentially expressed genes (DEGs) in each cell type after flg22 or Pep1 treatment. In total, 3276 unique DEGs responded to one or both elicitors in at least one cell type. Consistent with a recent study (Poncini et al., 2017), Pep1 treatment elicited markedly more DEGs (3082) in roots than flg22  (A) to (C) The top 5 nonredundant, most significantly enriched GO terms in cell type-specific identity genes. The 512 most significant upregulated genes from each cell type were used for consistent comparisons across cell types. Tone of color indicates significance degree [log 10 (P adj )].
(D) to (F) Venn diagrams showing the overlap between cell type-specific identity genes and cell type-specific, flg22or Pep1-reponsive immunity genes.
To determine if the cell type specificity of flg22and Pep1responsive gene networks reflects specific functions, we conducted Gene Ontology (GO) analyses. In the epidermis and cortex, flg22and Pep1-induced genes were enriched in immunityassociated terms (Figures 2C to 2F). In total, 21% of epidermisspecific (e.g., ). The GO term analysis revealed functional specificity; notably, immune and hormone responses were more pronounced in flg22or Pep1-treated epidermal cells (e.g., immune system process, hormone metabolic process; Figures 2C and 2E). In turn, flg22and Pep1-induced genes in cortex cells were significantly enriched for terms associated with transport processes (e.g., establishment of localization, oligopeptide transport, organic anion transport; Figures 2D and 2F; Supplemental Data Set 3). Due to the small number of genes (140 genes across all cell types), we were unable to conduct GO analyses for flg22-suppressed genes. However, GO terms associated with Pep1-repressed genes were enriched in terms associated with developmental processes in the epidermis (e.g., root morphogenesis, root system development; Figure 2G) and in flavonoid metabolism and growth hormone synthesis in the cortex (e.g., positive regulation of flavonoid, brassinosteroid biosynthetic process; Figure 2H).

flg22-Regulated Genes Are Largely Encompassed within a More Diverse Pep1 Gene Network
We next determined the commonalities of flg22 and Pep1 responses between treatments and across cell types ( Figure 2I). Aggregating across cell types, 78% (690 of 884 DEGs) of flg22responsive genes were also regulated by Pep1, whereas 22% (194 of 884 DEGs) were specific to flg22 ( Figure 2I). Of the 194 flg22specific DEGs, 89% (174 DEGs) were only expressed in one cell type (DEGs: 98 in epidermis, 49 in cortex, 27 in pericycle; Figure 2I; Supplemental Data Set 4). The majority of Pep1-responsive genes were specific to Pep1 (2392 of 3082 DEGs; 77%). These Pep1specific genes were also largely cell type-specific, with 1016 DEGs expressed in the cortex, 583 in epidermal cells, and 152 in pericycle cells ( Figure 2I; Supplemental Data Set 4). GO analyses of these elicitor-specific DEGs revealed that flg22-responsive genes in epidermal and cortical cells were enriched in GO terms defining transport processes (e.g., organic acid transport, amino acid transport, nitrate transport), whereas Pep1-enriched terms were associated with hormone metabolism and signaling (e.g., hormone metabolic process, salicylic acid biosynthetic process) in these two cell types (Supplemental Data Set 5). Overall, our results show that both flg22 and Pep1 activate distinct gene networks in each cell type and that the flg22 response is largely encompassed within a much more diverse Pep1 response.

Immunity Networks Only Partially Overlap with Cell Identity Gene Networks
Our observation that Pep1-repressed genes function in root growth and development ( Figure 2G) is in agreement with the reported root growth-inhibiting effect of PTI (Gómez-Gómez et al., 1999;Jacobs et al., 2011). It is known that the maintenance of the identity of each cell type, which is defined by cell type-specific functions, is essential for overall root integrity and root growth, especially under stress (Iyer-Pascuzzi et al., 2011;Geng et al., 2013). To determine if PTI affects cell (type) identity networks, we first defined the cell identity-specific transcriptomes (using our RNA-seq data from mock-treated samples) and identified 950 genes as specifically enriched in epidermis, 512 in cortex, and 1055 in pericycle (Supplemental Data Set 6). These enriched data sets were confirmed to strongly overlap (P < 10 26 , Fisher's exact test) with published cell identity gene sets (Bargmann et al., 2013), and distinct GO terms were associated with each set of identity genes specifying the different functions of each cell type (Figures 3A to 3C; Supplemental Data Set 6). By comparing cell identity with cell type-specific PTI gene sets (by combining flg22 and Pep1 DEGs per cell type), we found that PTI affected 18% (epidermis), 28% (cortex), and 5% (pericycle) of respective cell identity genes (Figures 3D to 3F; Supplemental Data Set 7). Similarly, salt stress or iron deprivation-regulated networks overlapped with cell identity networks, which was found to support cell type-specific responses to abiotic stresses (Dinneny et al., 2008;Iyer-Pascuzzi et al., 2011). We therefore wanted to understand how cell identity and immunity networks are linked. Heatmaps represent P values (log10 Padj ) for enrichment of motif pairs in the promoters of flg22-induced epidermis (A) and cortex genes (B), promoters of Pep1-induced epidermis (C) and cortex genes (D), and promoters of Pep1-suppressed epidermis (E) and cortex genes (F). x and y axes of all heatmaps show the same order of TF binding motifs. Pixel colors indicate significance values (log10 Padj ) for the abundance of given motif pairs in the promoters of identity genes per cell type.

Specific TF Pairing Links Cell Identity with Cell Type-Specific Immunity Networks
To regulate gene networks, TFs exert their activity at their site (cell type) of synthesis but can also move across cell boundaries. For instance, fundamental root developmental processes such as root patterning, cell fate decision, and root growth depend on the mobility of TFs such as SHORTROOT, KNOTTED1, and PHLOEM EARLY DOF (Nakajima et al., 2001;Xu et al., 2011;Clark et al., 2016;Miyashima et al., 2019). We therefore analyzed the presence and abundance of TF binding motifs in the promoters of DEGs to identify cell type-specific gene regulatory patterns and reveal any interdependencies between cell identity and cell type-specific immunity. Based on the clear evidence that combinatorial TF pairing is most crucial in controlling gene expression (Achard et al., 2009;Van de Velde et al., 2014;Lewis et al., 2015), we developed the Paired Motif Enrichment Tool (PMET) method to identify pairs of TF binding motifs within the promoters of our cell identity and cell type-specific DEG sets based on the following criteria, as depicted in Figure 4A: (1) a multiplicity of TF binding motifs in each promoter, and (2) acceptance of limited but (3) rejection of major overlapping of binding sites with (4) high motif specificities ( Figure 4A; see Methods). By tolerating limited motif overlap, our analysis allowed us to consider TF complexes associated with compound binding sites (Rodríguez-Martínez et al., 2017).
In the first step, we identified enriched motif pairs in our cell identity gene networks using a TF binding motif database derived from the work of Franco-Zorrilla et al. (2014). To enable direct comparisons with equal statistical power, we balanced the data sets in that we compared the promoters of all 512 (368 1 144 genes; Figure 3E) cortex identity genes with the promoters of the 512 most significantly upregulated epidermis and pericycle identity genes (Figures 3D to 3F). We found highly significant enrichment of a number of motif pairs within the promoters of 472 epidermis (92%), 442 cortex (86%), and 461 pericycle (90%) cell identity genes (Figures 4B to 4D; Supplemental Data Set 8). Comparing epidermis, cortex, and pericycle, the pattern of enriched motif combinations was distinctive between all cell types. For instance, within the promoters of epidermis identity genes, motifs for four different WRKY TFs (WRKY12/38/45 and to a lesser extent WRKY18) were found to pair uniquely with a wide range of motifs ( Figure 4B). In particular, WRKYs were enriched with binding motifs for AT-HOOK MOTIF CONTAINING NUCLEAR LOCALIZED (AHL), ARABIDOPSIS THALIANA HOMEOBOX (ATHB), and ARABIDOPSIS NAC (ANAC) TFs within the promoters of 40% of all epidermis identity genes (203 of 512 genes). In the cortex gene promoters, we observed specific pairing of ATHB with MYB TF binding motifs, whereas pairing between AHLs and DOF AFFECTING GERMINATION (DAG2) or ZINC FINGER OF ARA-BIDOPSIS6 (ZAT6) was pericycle-specific ( Figures 4C and 4D). We also noted some overlap in enriched motif pairs between epidermis and cortex identity networks, with MYCs and PHY-TOCHROME-INTERACTING FACTORs (PIFs) showing pairing with AHL and ATHB motifs for both cell types ( Figures 4B and 4C).
Interestingly, we observed the distinct pairing of stress-and development-associated TFs in all cell types. In the epidermis, for instance, AHL and ATHB TFs, which function in growth and development (Matsushita et al., 2007;Hur et al., 2015;Miao et al., 2018), paired with WRKYs, a large family of Arabidopsis TFs (>70 members) with regulatory functions in plant innate immunity and abiotic stress responses (Pandey and Somssich, 2009;Rushton et al., 2010), whereas in the cortex, AHLs and ATHBs paired with plant immune/jasmonate-responsive MYC2/3/4 TFs (Fernández-Calvo et al., 2011;Schweizer et al., 2013). Based on this apparent connection between developmental and stress networks within identity genes, we analyzed cell type-specific flg22and Pep1responsive promoter sets to test whether this connection was also present following immune activation. We had to exclude pericycle data (Figures 2A and 2B) from these analyses due to the low number of flg22or Pep1-responsive genes. To make the enrichment scores across treated cell types comparable, we again equalized gene set sizes (see Methods).
In summary, the analyses again revealed the interaction of stress-and development-associated TFs as a consistent theme, irrespective of flg22, Pep1, or cell type. Interestingly, in the promoters of Pep1-repressed genes, which were enriched in developmental GO terms ( Figures 2G and 2H), we detected a higher abundance of developmental TF motif pairing (with a complete absence of WRKY motifs). While Pep1-suppressed genes in the epidermis showed an enriched pairing of MYC2/3/4 and PIF3/4 with AHL12/20/25 and YAB1 motifs, the promoters of Pep1suppressed cortical genes revealed ATHB12/15/51 motifs paired with MYB46/52/111. Our paired-motif enrichment analyses thus revealed a pairing of stress-and development-associated TFs in regulating immunity networks in a cell type-specific manner. Overall, for flg22, our promoter analyses uncovered an enrichment of highly specific motif pairs in the promoters of 90% of flg22induced epidermis genes and 86% of cortex genes ( Figures 5A  and 5C). For Pep1, we identified paired-motif enrichment in the promoters of 58 and 60% of epidermis-or cortex-induced genes, respectively ( Figures 5B and 5D) as well as 57 and 50% of epidermis-or cortex-suppressed genes, respectively (Figures 5E and 5F; Supplemental Data Set 10). Despite some overlap, the promoters of DEGs specifically regulated by Pep1 and flg22 generally In (A) to (I), native promoters with or without deletions in predicted TF binding motifs were analyzed in two independent plant lines stably expressing the Promoter:LhG4 > pOp6:YFP transactivation system. Please see the results for both lines per construct in Supplemental Figure 8. SCRI Renaissance 2200 (SR2200) was used to stain root cells. showed clear differences in motif enrichment within and across cell types.

TF Binding Motif Pairs Determine Cell Type-Specific and Elicitor-Responsive Gene Regulation
As a final step, we wanted to confirm the accuracy of the predictions from our promoter motif enrichment tool in interpreting cell type-specific expression patterns. We first identified a set of native promoters from genes that showed cell type-specific expression and enrichment for a motif pair. For the epidermis, WRKY45 and PLANT INTRACELLULAR RAS GROUP-RELATED LRR2 (PIRL2) showed specific expression and predicted KAN-WRKY pairs. AtM10 and BASIC HELIX-LOOP-HELIX92 (bHLH92) were chosen as cortex-specific genes enriched in MYC-WRKY pairs in their promoters. For our analyses we used the pPRO-MOTER:LhG4 > pOp6:YFP transactivation system (Moore et al., 1998;Craft et al., 2005;Costa et al., 2014), where our native promoters were fused to LhG4 (pPROMOTER NATIVE :LhG4), which binds to the Op6 promoter (pOP6) to transactivate yellow fluorescent protein (YFP) expression. For each promoter variant, we identified two independent Arabidopsis lines for further analyses. Excitingly, except for bHLH92 constructs, where we could not identify any lines, we confirmed the predicted cell type-specific expression patterns in planta: pOp6:YFP lines expressing pWRKY45 NATIVE :LhG4 or pPIRL2 NATIVE :LhG4 showed epidermisspecific YFP expression, whereas pAtM10 NATIVE :LhG4 mediated cortex-enriched YFP expression ( Figures 6A, 6D, and 6G; Supplemental Figures 7A to 7F and 8A, 8D, and 8G). Moreover, replacing one type of the predicted motifs in the otherwise unaltered native promoters by a nonfunctional sequence (see Methods for design) was sufficient to abolish YFP expression ( Figures 6B, 6C, 6E, 6F, 6H, and 6I; Supplemental Figures 8B, 8C, 8E, 8F, 8H, and 8I), indicating the accuracy of our prediction tool in identifying motif pairs that determine cell type-specific expression. We further confirmed the predicted function of these motif pairs in cell type-specific regulation by flg22 or Pep1. Again, eliminating one type of motif was sufficient to abolish (or significantly reduce Pep1-induced YFP transactivation in plants expressing pWRKY45 NATIVE_DWRKY ) cell type-specific elicitor inducibility ( Figures 6J to 6L).
In the case of plants expressing pPIRL2 NATIVE_DKAN (lacking KAN motifs in the native PIRL2 promoter), we detected enhanced basal and elicitor-induced, cell type-specific expression of YFP (under mock,flg22,or Pep1;Figures 6D to 6F;Supplemental Figures 8D to 8F). This led us to investigate the roles of the identified motifs as activators or repressors of gene transcription. We therefore generated synthetic promoters consisting of motif sequences linked with short nonfunctional linker sequences (see Methods for design) and, to exclude any other regulatory motifs, lacking any promoter backbone sequence. These synthetic promoters were run by a proximally placed minimal CaMV35S promoter (CaMV35S min ), and plants expressing this minimal CaMV35S promoter alone (pCaMV35S min :LhG4) did not induce YFP expression. To validate this positional effect further, we generated a second set of synthetic promoters with 43 WRKY and 43 MYC motifs. When 43 MYC motifs were placed proximally to 43 WRKY motifs (pWRKY-MYC-Motifs SYNTHETIC :LhG4), the plants showed reduced YFP expression ( Figure 7G), whereas full epidermis-specific YFP expression levels were detected when 43 MYC motifs had a distal position to 43 WRKY motifs (pMYC-WRKY-Motifs SYNTHETIC :LhG4; Figure 7H; Supplemental Figures 9A and 9B). This indicates the importance of motif positions in balancing gene expression and that KAN motifs apparently functioned as repressor elements in our setup.
Taken together, our promoter motif prediction tool identified paired TF binding sites, and we confirmed the importance of the   predicted pairs of TF binding motifs in regulating cell type-specific gene expression. This tool is also suitable for identifying relevant promoter motifs whose function can be studied in synthetic promoters as a helpful step in generating customized or minimal synthetic promoters for the targeted analysis of gene network regulation in plants.

DISCUSSION
The physical inaccessibility of roots impedes the detection and control of soil-borne pathogens and explains the high relevance of root diseases for staple crop production (Gewin, 2010;Alexandratos and Bruinsma, 2012;Panth et al., 2020). Under the projected future global warming, the frequency of root diseases is expected to rise at a global scale (Delgado-Baquerizo et al., 2020). Therefore, new control strategies are needed, which requires a better understanding of the genetic disease resistance potentials and regulatory mechanisms of underlying immune responses in roots. In this study, we analyzed the immunity gene networks of epidermis and cortex cells, which build the outer frontier to the rhizosphere, as well as the pericycle, the "outer frontier" of the inner root vasculature. To obtain a comprehensive picture of root immunity, we analyzed cell type-specific PTI in response to the MAMP flg22, which activates PTI against bacteria, and the DAMP Based on our paired motif enrichment analyses, TFs can be assigned to cell identity or cell type-specific immunity networks, and immunity networks differ between epidermis and cortex cells. Those TFs that can be found in both (identity and immunity) networks might act as core TFs (bracketed) to connect epidermis identity with cell type-specific immunity networks. Among all core TFs, WRKY12/18/38/45 and MYC2/3/4 together with ATHB15/51 showed the strongest links in the epidermis and cortex, respectively. Synthetic promoters only consisting of a set of four motifs or motif combinations and a CaMV35S minimal promoter (but lacking any native promoter backbone sequence) were analyzed in plant lines stably expressing the Promoter:LhG4 > pOp6:YFP transactivation system. Except for WRKY-KANADI motif combinations, two independent plant lines were analyzed per promoter construct by confocal microscopy of axial root sections. Pep1, a plant-derived PTI elicitor that is activated upon the perception of different microbes and defense hormones and, hence, might trigger the full array of all PTI responses against a larger variety of pathogens (Ryan et al., 2007;Lori et al., 2015). By combining RNA-seq data with PMET, a promoter analysis approach, we were able to identify distinct regulatory patterns within immunity networks activated in a cell type-specific manner.
The analytic design of PMET considers the pairing of TF binding sites on promoters, which proved to be sufficient to explain the regulatory patterns of 50 to 90% of DEGs per cell type and treatment. The application of PMET in the analysis of complex transcriptomes is not restricted to plants but works for eukaryotes in general. We found it equally efficient for elucidating network regulation in neuroinflammatory disorders in mouse (Schang et al., 2018). In this study, PMET revealed the close interplay of cell typespecific immunity with identity networks that are essential for tissue organization and development. Consistent with this finding, cell type-specific transcriptomics revealed a close interaction of abiotic stress with cell identity networks. The adaptation to different abiotic stresses includes a highly coordinated and cell typespecific redirection of gene networks to maintain root function (e.g., growth; Dinneny et al., 2008;Gifford et al., 2008;Iyer-Pascuzzi et al., 2011). In addition, the underlying rebooting of cell type-specific signaling under salt stress involved hormones, such as abscisic acid, to readjust root growth and adapt root system architecture (Geng et al., 2013). Besides hormones, RAPID ALKALINIZATION FACTOR (RALF) peptides in interaction with their receptors THESEUS1 and FERONIA regulate lateral root formation (Haruta et al., 2014;Murphy et al., 2016;Gonneau et al., 2018;Jourquin et al., 2020). Interestingly, RALF-FERONIA complexes also control PTI (Stegmann et al., 2017). These studies suggest a close interdependency of developmental and stress networks at the perception (receptor) to signal transduction level (e.g., hormones) in order to adjust overall root development to changing environments. The observed interplay of cell identity and cell type-specific immune gene networks in our study adds regulatory mechanisms at the gene level to support this crosstalk at the plant development-environment interface. Altogether, these findings demonstrate a remarkable complexity of immunity signaling in roots.
Although cell type specificity in immunity gene regulation may require a higher degree of coordination (e.g., numerous regulatory and signaling proteins), it apparently adds to the robustness and flexibility required for a root system to adapt to changing environments. Accordingly, recent studies observed flg22 responsiveness in cell types of different developmental ages, although qualitative differences in the immune competences of root cell types appeared to occur (Beck et al., 2014;Wyrsch et al., 2015;Poncini et al., 2017). Supported by PMET-base analyses, we obtained insights into how this tight coordination of cell typespecific immune responses may be achieved. Our observation that root cell types keep their identity under biotic stress, as indicated by cell type marker expression, principal component analyses, and studying the regulation of cell identity genes under immunity ( Figures 1C, 1D, and D to F3D to 3F), might be most critical in this respect. We found that the pairing of TF motifs for DEGs differed depending on the cell type, and we confirmed the significance of TF pairs for cell type-specific expression of immunity genes in our functional promoter analyses (Figure 6; Supplemental Figures 7 and 8). Moreover, our data suggest a model where certain TF combinations prevailed in specific cell types in a treatment-dependent manner and revealed "core" TFs that link cell identity with cell type immunity networks (Figure 8). For the epidermis-specific networks, WRKY12/18/38/45 appear to act as core TFs (together with ANAC46/55/58 and AHL12/20/ 25) to connect identity and immunity networks by pairing with ATHB, MYC, and PIF TFs of the epidermis identity network and with KAN and YAB TFs of the epidermis immunity network. In the cortex, in turn, MYC2/3/4 and ATHB15/51 (together with PIF3/4/5 and AHL12/20/25) might serve as core TFs that pair with WRKYs to direct cortex-specific immunity as well as with ATHB12 and MYBs to corroborate cortex identity networks. Thus, the cooperation of different TF families in specific combinations underpins the highly cell type-specific immunity networks that we observed.
In addition to defining links between cell identity and immunity networks, our promoter motif prediction tool can disentangle apparently contradictory functions of TFs. We observed a striking difference in the context-dependent association of MYC motifs in Pep1 responses. MYC TFs have been implicated in Pep1mediated signaling. Pep peptides specifically activate the MYC2-dependent branch of jasmonic acid signaling , and MYC2 can act as both an activator and a repressor of jasmonic acid-mediated gene expression (Dombrecht et al., 2007). Consistent with this observation, MYC2/3/4 binding motifs were found to pair with WRKY and ATHB motifs in the promoters of Pep1 upregulated cortex genes ( Figure 5D). By contrast, MYC2/3/4 motifs paired with AHL12/20/25 motifs in the promoters of downregulated genes in the epidermis ( Figure 5E). The analyses thus suggest a potential multifunctionality of specific TF family member combinations, whereby they can act with heterogeneous partners in different cell types to control contrasting transcriptional outputs. In addition, the observed context dependency of regulatory functions may reduce the perceived redundancy among TFs that recognize highly similar sites if investigated in isolation.
By computing the presence of TF binding motif combinations for DEGs in our RNA-seq data, we observed that cell identity and stress-responsive gene networks coexist in each root cell type. As reported for abiotic stress integration (Dinneny et al., 2008;Iyer-Pascuzzi et al., 2011;Geng et al., 2013), our data support the concept that cell identity underpins transcriptional reprogramming leading to cell specificity in response to signal perception, even modulating outputs from strong elicitors such as MAMPs. Linking immunity to cell identity networks would guarantee cell type-specific coordination of immune responses with individual and measured contributions from each cell type. Such a coregulatory model would likely be applicable to root responses to all environmental stresses and add the necessary flexibility to gene regulation (Dinneny et al., 2008;Van de Velde et al., 2014). In the future, it will be interesting to explore to what extent cell identity and core TFs contribute to cell type specificity of immune responses and resistance against pathogens as well as root growth upon the activation of immunity. Our analyses using promoter lines ( Figure 6; Supplemental Figures 7 and 8) suggest the feasibility of such functional dissection studies.

Plant Materials, Growth Conditions, and Treatments
Seeds of Arabidopsis (Arabidopsis thaliana) ecotype Col-0 were obtained from the Nottingham Arabidopsis Stock Center. Seeds of the pepr1-2 pepr2-2 mutant were kindly provided by Y. Yamaguchi (Osaka Prefecture University). Marker lines for the three cell types, pGL2:GFP (epidermis atrichoblast; Lin et al. 2015), pCORTEX:GFP (cortex), and E3754 (xylempole pericycle; Bargmann et al. 2013), were obtained from Miriam Gifford (University of Warwick). PTI marker lines pMYB51:NLS-3xmVENUS and pPER5:NLS-3xmVENUS were provided by Silke Lehmann (University of Warwick; Poncini et al., 2017). Plants were grown on vertical square Petri dishes on ATS medium (Lincoln et al., 1990)  For experiments with Serendipita indica, roots of 9-to 10-d-old Arabidopsis plants were inoculated with 1 mL of a 500,000 chlamydospores mL 21 spore suspension per Petri dish. Control plants were treated with water containing 0.02% Tween 20 (mock). If not stated otherwise, plants were treated on plates with 1 mL per plate with a 1 mM solution of flg22 or Pep1 or with water as a control. For all experiments, flg22 and Pep1 peptides were used as described (Gómez-Gómez et al., 1999;Krol et al., 2010). All data are based on at least three independent biological experiments (if not specified otherwise).

Measurement of Plant Growth Inhibition
For seedling growth inhibition assays, plants were grown on square Petri dishes for 10 d before inoculation with S. indica or mock treatment. After 3 d, plants were supplemented with 1 mM flg22, 1 mM Pep1, or water (control). Eleven days later, plant fresh weights were determined. For root growth inhibition assays, pMYB51:NLS-3xmVENUS and pPER5:NLS-3xmVENUS lines were grown on square Petri dishes for 10 d and treated with 1 mL per plate of 1 mM flg22 or 1 mM Pep1 solution in water. Control plants were treated with water (mock). Four days later, the plates were photographed and root lengths were measured using ImageJ software (https://imagej. net). If not stated otherwise, all biological experiments were done in triplicate and at least 12 plants per line per treatment were evaluated.

Measurement of ROS Burst
Roots of 2-week-old plants were grown on solid ATS medium and treated with 1 mM flg22, 1 mM Pep1, or mock at 3 d after inoculation with S. indica or mock treatment. For ROS burst quantification, roots were cut into 1-cmlong pieces (10 mg per assay) and transferred to a luminol-based assay as described (Gómez-Gómez et al., 1999). Data were analyzed by Student's t test.

Gene Expression Analysis Using qRT-PCR
For gene expression analyses of whole roots, root material was harvested 2 and 24 h after flg22, Pep1, or control treatment. TRIzol (Invitrogen) was used to extract total RNA. After DNase treatment, RNA was reverse transcribed into cDNA using a qScript cDNA synthesis kit (Quanta Biosciences), and 10 ng of cDNA was used as a template in qRT-PCR using SYBR Green JumpStart Taq ReadyMix (Sigma-Aldrich) and a Stratagene Mx3005P Real-time PCR Detection System (Agilent Technologies) following the manufacturer's recommended protocol. The 2 2DCt method (Schmittgen and Livak, 2008) was used to determine the differential expression of marker genes for immunity activation, and the housekeeping genes UBIQUITIN5 (UBQ5, AT3G62250) and ELONGATION FACTOR1a (EF1a, AT5G60390; for primer sequences, see Supplemental Table 5) were used for normalization. Data were analyzed using Student's t test.

Quantification of S. indica Colonization by qRT-PCR
Genomic DNA was isolated from roots using a Plant DNeasy Kit (Qiagen). A total of 40 ng of genomic DNA served as the template in qRT-PCR using SYBR Green JumpStart Taq ReadyMix (Sigma-Aldrich) in a Stratagene Mx3005P Real-time PCR Detection System (Agilent Technologies) following the manufacturer's recommended protocol. The 2 2DCt method (Schmittgen and Livak, 2008) was used to determine the extent of fungal colonization by subtracting the raw cycle threshold values of S. indica internal transcribed spacer from those of Arabidopsis UBQ5 (for primer sequences, see Supplemental Table 5). Differences in S. indica colonization were determined by Student's t test.

FACS
For FACS experiments, plants were grown for 12 d on square ATS plates and treated with 1 mL per plate of 1 mM solutions of flg22 or Pep1 peptide or water as a control for 1 h. Taking into account the time required for protoplast generation and cell sorting (together ;1 h), the status of all sampled cells was 2 h after flg22 or Pep1 treatment. Briefly, whole roots were cut into pieces and incubated in protoplast solution (1.5% cellulase R10 [Duchefa Biochemie], 1.2% cellulase RS [Duchefa Biochemie], 0.2% macerozyme R10 [Duchefa Biochemie], and 0.12% pectinase [Sigma-Aldrich] in 600 mM mannitol, 2 mM MES hydrate, 10 mM KCl, 2 mM CaCl 2 , 2 mM MgCl 2 , and 0.1% BSA, pH 5.7; Walker et al., 2017) for 45 min. Protoplasts were filtered through 70-mm followed by 40-mm cell strainers, centrifuged at 300g for 3 min, resuspended in protoplast solution lacking cell wall-degrading enzymes, and subjected to FACS. Three independent biological experiments were performed for each marker line. GFP-expressing protoplasts were collected using a BD Influx cell sorter (BD Biosciences) following previously published protocols (Birnbaum et al., 2003;Gifford et al., 2008;Grønlund et al., 2012). The cell sorter was equipped with a 100-mm nozzle, and BD FACSFlow (BD Biosciences) was used as the sheath fluid. BD Accudrop Fluorescent Beads (BD Biosciences) were used prior to each experiment to optimize sorting settings. A pressure of 20 p.s.i. (sheath) and 21 to 21.5 p.s.i. (sample) was applied during the experiments. Drop frequency was set to 39.2 kHz, and event rate was generally kept at <4000 events s 21 . GFP-expressing protoplasts were identified using a 488-nm argon laser, plotting the outcome of a 580/30 bandpass filter versus a 530/ 40 bandpass filter, to differentiate between green fluorescence and autofluorescence. Different cell populations were collected for microscopy in preexperiments to determine the presence of GFP-expressing protoplasts. As previously reported (Grønlund et al., 2012), these protoplasts were present in the high 530-nm/low 580-nm population. Sorting gates were set conservatively in subsequent experiments based on these observations (Supplemental Figure 10). For RNA extraction, GFP-expressing protoplasts were sorted into Qiagen RLT lysis buffer containing 1% (v/v) b-mercaptoethanol, mixed, and immediately frozen at 280°C. At least 10,000 GFP-expressing protoplasts were sorted per experiment and treatment condition. Sorting times were kept below 20 min.

RNA Isolation, RNA-Seq Library Construction, and Sequencing
Total RNA was extracted from the samples using a Qiagen RNeasy Plant Mini Kit including on-column DNase treatment with a Qiagen DNase kit. A 6000 Pico Kit (Agilent Technologies) was used to check quantity and quality of the RNA on a Bioanalyzer 2100 (Agilent Technologies). Preparation of amplified cDNA from total RNA and RNA-seq library construction were performed using the Ovation RNA-seq System V2 and Ovation Ultralow Library Systems Kits (NuGEN Technologies), respectively, following standard protocols. Sequencing was performed by the High-Throughput Genomics Group at the Wellcome Trust Centre for Human Genetics on an Illumina HiSeq2500 System.

RNA-Seq Quality Control and Read Mapping
For each sample, read quality was evaluated using FastQC software (Andrews, 2010). The paired-end libraries (2 3 100-bp reads) were mapped to the Arabidopsis TAIR10 genome using STAR (default parameters; Dobin et al., 2013). The reads mapping to exons were counted using LiBiNorm (Dyer et al., 2019; settings: -f bam -s no -i Parent -t mRNA) using an Araport 11 annotation GTF file (Ensembl release 39). On average, 43% of reads uniquely mapped to exons (full details are given in Supplemental Table 1). The quality of read mapping was assessed using the Integrative Genomics Viewer (Robinson et al., 2011). The quality of replicates was assessed by plotting read counts of samples against one another and assessing the dispersion and presence of any artifacts between samples. Due to preferential amplification in some samples, reads corresponding to rRNA and ribosomal proteins had to be removed for subsequent analyses (Supplemental Data Set 11). The mitochondrial and plastid chromosomes were also removed, as this work focused on nucleus-encoded genes. Principal component analysis was calculated using the R function prcomp and visualized using ggplot2 implemented in R. Fragments per kilobase per million values were calculated for genes remaining after filtering using exonic gene lengths from the Araport 11 annotation GTF file (the same GTF file used for HTSeq-count; Supplemental Data Set 12).

Differential Gene Expression and Functional Analysis
Within each cell type, genes that were differentially expressed in response to each treatment (compared with mock) were identified using DESeq2 (Love et al., 2014). All counts (minus filtered genes described above) were normalized using DESeq2 (default parameters). The model for differential expression included the replicate information in the model matrix (;batch 1 condition) in order to account for batch effects. Genes were defined as differentially expressed if the adjusted P value was less than 0.05. We included all statistically significant genes and did not use a threshold on the fold change to (1) avoid the use of an arbitrary threshold, especially as thresholds do not necessarily correlate with translation and/or biological function, and (2) maximize statistical power for downstream analyses such as GO term or motif enrichment analyses. In order to define cell identity genes, first we identified genes that were differentially expressed between all possible pairs of cell type mock-treated samples using DESeq2 (again accounting for batches in the model matrix; P < 0.05). We then defined cell identity genes as those that were significantly upregulated in one cell type compared with both other cell types. For example, to define epidermis identity genes, we first identified genes differentially upregulated compared with the pericycle and cortex separately, then we took the intersection of these two lists of genes to be the epidermis identity genes. The significance of the overlap between cell identity genes and published cellspecific gene sets (Bargmann et al., 2013;Supplemental Data Set 13) was tested using Fisher's exact test. From the published data set, epidermisspecific genes were taken to be the union of those labeled "trichoblast" and "atrichoblast," cortex genes were those labeled as "cortex," and pericycle genes were labeled as the union of "phloem-pole pericycle" and "xylempole pericycle" genes (Supplemental Data Set 14). The fit of the DESeq2 model to our data was tested by plotting replicates against one another and overlaying the read counts of DEGs (Supplemental Figure 4). Subset analysis to determine cell type-and treatment-exclusive genes was performed in R using built-in set functions and the VennDiagram (Chen and Boutros, 2011) package. Proportional visualizations of three-set Venn diagrams were created using eulerAPE (Micallef and Rodgers, 2014), nonproportional Venn diagrams were created using the Venn function from the R packaged gplots, and the six-set Venn diagram was created using the interactive Venn tool by Bardou et al. (2014).
GO enrichment analysis was performed using the R package GOStats (Falcon and Gentleman, 2007) with an additional Benjamini-Hochberg multiple testing correction applied. To make the P values across cell types directly comparable, we equalized gene set sizes by taking the top K genes from the larger data set where K is the size of the smaller data set for flg22 upregulated genes (cortex set: 128 genes; Figure 2A), Pep1 upregulated genes (epidermis set: 365 genes; Figure 2B), and Pep1 downregulated genes (epidermis set: 337 genes; Figure 2B). For the cell identity genes, we took the top 512 genes ( Figure 3E). For the figures, the top 5 nonredundant terms were shown, where redundant terms were defined as terms that corresponded to exactly the same sets of genes. These were not shown in the plots in order to maximize the scope of the top GO terms.

PMET
A total of 113 motifs (in the form of letter-probability matrices) were obtained from microarray studies performed by Franco-Zorrilla et al. (2014). These motifs characterize the target sequence specificity of 63 plant TFs representing 25 families. To conduct accurate comparisons, we equalized gene set sizes by taking all of the cortex-specific flg22 upregulated genes (128 genes; from Figure 2A) and the top 120 most significantly upregulated epidermis-specific genes. For Pep1, we included all Pep1-induced/suppressed epidermis genes (365 and 337 genes, respectively; Figure 2B) and the equivalent top Pep1-induced/suppressed genes in the cortex. Promoter regions corresponding to 1000 bp upstream from the transcription start site were collected from the TAIR database (Lamesch et al., 2012; www.arabidopsis.org) for all nuclear genes (including transposable element genes and excluding genes on the mitochondrial or plastid chromosomes) in the Arabidopsis genome. For each motif and each promoter, the sequence was scanned for occurrences of the motif using FIMO (Grant et al., 2011), which assigns a probability score to each potential hit. In order to determine the number of hits to consider and to compute an overall score for motif presence in a promoter, we computed the geometric mean p of the top k FIMO probability scores for nonoverlapping hits and computed the binomial probability of observing at least k hits of probability p in a 1-kb promoter. The value of k minimizing the binomial probability was taken to indicate the most likely number of binding sites, and the k hits were recorded for subsequent analysis (1 # k # 5). For each motif, the promoters were ordered by increasing binomial probability, and the top n 5 5000 promoters were considered as containing the motif. Homogeneously setting n to the same value for each motif lays the foundation for computing P values for motif pairs that are comparable (see below) and avoids biases brought about by motifs of varying specificity. The parameter n was set to the high value of 5000 for high sensitivity (rather than specificity), as stringency is introduced when the pairing of motifs is considered. The binomial probability of the n-th promoter was recorded for each motif as a threshold. For each pair of motifs and for each promoter containing both motifs, overlaps of recorded motif hits were identified and the total information content (IC) of the overlap (based on motifs) was calculated. Total IC across the overlapping positions (i) of one motif was calculated as where f b,i corresponds to the probability of observing base b at position i, according to previous work by Franco-Zorrilla et. al. (2014), and where B 5 fA, C, G, Tg and I is the set of positions in the overlap. If the IC of the overlap for either motif exceeded 4 (indicating that highly conserved bases are part of the overlap), these hits were removed and the binomial probability was recalculated for the remaining hits. The process of overlap removal is, therefore, only dependent on the IC and not directly dependent on the length of the overlap. If the recalculated scores were still below the recorded motif-specific threshold, the two motifs were considered to be colocalized in the promoter. Finally, gene sets of interest were tested for enrichment of paired motifs using a pairwise hypergeometric test based on the MATLAB function proposed by Meng et al. (2009). Hypergeometric P values were corrected for the number of motif pairs using a stringent Bonferroni correction (calculating the correction for each gene set separately). As a negative control, random gene sets were tested against the same statistical method, resulting in no significant P values. Corrected P < 0.05 values are considered significant. For each comparison of results made between conditions, the gene sets tested were of equal size to make P values comparable. To this end, gene set sizes were equalized by taking the top K genes from the larger gene set, where K is the size of the smaller gene set.

Generation of LhG4/pOp6 Transactivation Lines to Study Promoter Activities
In order to assess the PMET results, we decided to validate two pairs of motifs: WRKY-KAN enrichment in epidermis-specific flg22-responsive genes and MYC-WRKY enrichment in cortex-specific Pep1-responsive genes. We used short native promoters of three candidates, all of which were expressed strongly in a cell type-specific manner. WRKY45 and PIRL2 were chosen as genes to validate WRKY-KAN enrichment and AtM10 was chosen to validate MYC-WRKY enrichment. Native promoter fragments (300 to 550 bp) were chosen such that they contained at least three binding sites for each TF type (WRKY, KAN, or MYC) but were short enough to reduce complexity arising from the presence of other TF binding motifs known to be present based on prior motif results. In the mutant versions of the promoters, WRKY, KAN, or MYC TF binding sites were replaced by the random sequence GAACTT. For synthetic promoters, we stacked TF binding sites in quadruplicate, each separated by 6 bp of a random spacer sequence (GAACTT), in front of the CaMV35S minimal promoter (nucleotides 246 to 11; Bhullar et al., 2003). In all cases, the presence or absence of TF binding sites was confirmed using FIMO, and we made sure that the insertion of random sequences did not result in the prediction of new TF binding sites. Native promoters or their mutated versions as well as synthetic promoter fragments were synthesized (Integrated DNA Technologies) flanked by attb sites for recombination into the Gateway-compatible pBIN1-GW-LhG4 vector (Costa et al., 2014). The resulting pBIN1Promoter-LhG4 "effector" plasmids were transformed into a stable pOp6:YFP "activator" line (Costa et al., 2014) via Agrobacterium tumefaciens-mediated transformation (Chang et al., 1994). Transformants were selected on hygromycin-containing medium, and the presence of the effector constructs was confirmed by PCR. For all assays, plants were grown in vertical square Petri dishes on ATS medium for 7 d and treated with 1 mL per plate of 1 mM flg22 or Pep1 or water (mock). Two hours later, whole roots of ;30 plants per line and treatment were harvested for expression analysis by qRT-PCR (with UBQ5 and EF1a for normalization as described above) or whole seedlings were harvested for confocal laserscanning microscopy. The 2 2DDCt method was used to determine YFP expression in transactivation plants relative to Col-0. To quantify YFP fluorescence in the synthetic promoter lines, maximum projections of the YFP channel from Z-stacks taken from root segments were analyzed using the ImageJ software package (https://imagej.net), and mean gray values for each image were determined. One-way ANOVA with Bonferroni posthoc test was applied to determine significance levels between native and synthetic promoter lines (see the Supplemental File for ANOVA tables).

Confocal Laser-Scanning Microscopy
All images were taken with a confocal laser-scanning microscope (Zeiss LSM 880 or Leica SP5). For all experiments, GFP was excited with a 488-nm laser line and detected between 500 and 545 nm, and YFP and VENUS were excited at 514 nm and detected between 520 and 560 nm. pMYB51:NLS-3xmVENUS, pPER5:NLS-3xmVENUS, pGL2:GFP, pCORTEX:GFP, and E3754 lines were imaged 1 h after flg22, Pep1, or mock treatment. For immunity and cell type-specific marker lines, viable roots were stained in a 10 mg mL 21 propidium iodide solution. Propidium iodide was excited at 561 nm, and fluorescence was detected between 570 and 720 nm. Due to the short time frame (2 to 3 h after treatment) for imaging of LhG4/pOp6 transactivation lines, we used a SCRI Renaissance 2200 (SR2200)/ ClearSee staining method (Kurihara et al., 2015;Musielak et al., 2015), which involves fixation in paraformaldehyde. For this, seedlings were incubated in SR2200 staining solution (Musielak et al., 2015) for 30 min, followed by a washing step with PBS and clearing for at least 2 d in ClearSee solution (Kurihara et al., 2015). SR2200-stained cell walls were imaged using a 405-nm laser line for excitation and a bandwidth between 410 and 480 nm for detection.

Uploading to ePlant
A scalable vector graphics (SVG) image representing our experimental setup was generated using the open-source SVG editing program Inkscape (https://inkscape.org). Expression data summarized as fragments per kilobase per million reads mapped were data based on the Bio-Analytic Resource for Plant Biology, and a web service was created to return expression values for a given gene across all samples in JSON (Javascript object notation) format. An XML file for a new Root Immunity Elicitation ePlant view was created to permit mapping of sample names to specific SVG parts and incorporated into ePlant (Waese et al., 2017). The appropriate group tags in the SVG file were edited to correspond to those in the XML file. The view may be freely accessed with ePlant at http://bar.utoronto. ca/eplant/?ActiveSpecies5Arabidopsis%20thalianaandGenes5AT2G35980 andActiveGene5AT2G35980andActiveView5RootImmunityElicitationView (view for NHL10).

Accession Numbers
The following supplemental data sets have been submitted to the Gene Expression Omnibus and are available at https://www.ncbi.nlm.nih.gov/ geo/query/acc.cgi?acc5GSE112960: paired-end raw .fastq files and raw counts produced using LiBiNorm (Dyer et al., 2019) as .txt files. The data generated and analyzed during this study are available from the corresponding author on reasonable request. PMET has been made available as a web-based tool at http://bar.utoronto.ca/index.html (linked to http://nero.wsbc.warwick.ac.uk/tools), and the PMET source code is available on GitHub at https://github.com/kate-wa/PMET-software.

Supplemental Data
Supplemental Figure 1. The ability of S. indica to suppress flg22 but not Pep1-triggered immune responses in roots.  Table 1. Alignment statistics table based on output from htseq-count.

Supplemental
Supplemental Table 2. Number of DEGs responding to flg22 and Pep1 in three root cell types.
Supplemental Table 3. Genes differentially expressed in all three cell types in response to both flg22 and Pep1.
Supplemental Table 4. Number of DEGs responding specifically to either flg22 or Pep1 in one or more cell types.
Supplemental Data Set 1. Lists of all differentially expressed genes in response to flg22 and Pep1 from DESeq2 output.
Supplemental Data Set 2. Lists of cell type-specific DEGs following flg22 or Pep1 treatment.
Supplemental Data Set 3. Lists of cell type-specific GO terms after flg22 and Pep1 treatment.
Supplemental Data Set 4. Lists of flg22and Pep1-specific DEGs per cell type.
Supplemental Data Set 5. Lists of flg22and Pep1-specific GO terms per cell type.
Supplemental Data Set 6. Lists of cell identity genes and enriched GO terms.
Supplemental Data Set 7. Lists of PTI genes within cell type-specific cell identity genes.
Supplemental Data Set 8. Paired-motif analysis results for cell identity genes.
Supplemental Data Set 9. Paired-motif analysis results for flg22 upregulated genes.
Supplemental Data Set 10. Paired-motif analysis results for Pep1 upand down-regulated genes.
Supplemental Data Set 11. Genes omitted from the analysis.
Supplemental Data Set 12. Fragments per kilobase per million (FPKM) matrix of filtered genes.
Supplemental Data Set 13. Lists of differentially expressed genes between mock-treated cell types.
Supplemental Data Set 14. List of cell-specific genes as published (Bargmann et al., 2013).

ACKNOWLEDGMENTS
We thank Jim Beynon and Murray Grant for helpful comments on the manuscript, Jeanette Selby and Lesley Ward at the Genome Research Facility (University of Warwick, School of Life Sciences) for their help in RNA-seq library preparation, Christina Neumann and Rebekka Schmidt (Institute of Phytopathology, Justus Liebig University Giessen) for help with PTI assays, and Jose Gutierrez-Marcos for sharing the LH4/pOP6 transactivation system. We thank the High-Throughput Genomics Group at the Wellcome Trust Centre for Human Genetics for the generation of sequencing data. This work was supported by the Biotechnology and Biological Sciences Research Council (Grant BB/M017982/1 to P.S. as well as to C. R.-G. through MIBTP), the Deutsche Forschungsgemeinschaft (Grants SCHA1444/3-3 to SCHA1444/5-2 to P.S.), and the Alan Turing Institute (Grant R-WAR-006 to P.E.B.).