Co-Transcriptomic Analysis of the Maize–Western Corn Rootworm Interaction

The Western corn rootworm (WCR; Diabrotica virgifera virgifera) is an economically important belowground pest of maize. Belowground feeding by WCR is damaging because it weakens the roots system, diminishes nutrient uptake, and creates entry points for fungal and bacterial pathogens and increases lodging, all of which can significantly suppress maize yields. Previously, it was demonstrated that belowground herbivory can trigger plant defense responses in the roots and the shoots, thereby impacting intraplant communication. Although several aspects of maize-WCR interactions have been reported, co-transcriptomic remodeling in the plant and insect are yet to be explored. We used a maize genotype, Mp708, that is resistant to a large guild of herbivore pests to study the underlying plant defense signaling network between below and aboveground tissues. We also evaluated WCR compensatory transcriptome responses. Using RNA-seq, we profiled the transcriptome of roots and leaves that interacted with WCR infestation up to 5 days post infestation (dpi). Our results suggest that Mp708 shoots and roots had elevated constitutive and WCR-feeding induced expression of genes related to jasmonic acid and ethylene pathways, respectively, before and after WCR feeding for 1 and 5 days. Similarly, extended feeding by WCR for 5 days in Mp708 roots suppressed many genes involved in the benzoxazinoid pathway, which is a major group of indole-derived secondary metabolites that provides resistance to several insect pests in maize. Furthermore, extended feeding by WCR on Mp708 roots revealed several genes that were downregulated in WCR, which include genes related to proteolysis, neuropeptide signaling pathway, defense response, drug catabolic process, and hormone metabolic process. These findings indicate a dynamic transcriptomic dialog between WCR and WCR-infested maize plants.


Introduction
Western corn rootworm (WCR, Diabrotica virgifera virgifera LeConte) is a damaging pest of maize (Zea mays L.) across North America, Central and Eastern Europe [1,2]. WCR, at their larval stage, will start to damage maize roots and, more specifically, the root hairs and cortical tissues to create tunnels to reach the primary root [3][4][5][6][7]. Roots, with their root hairs, are the part of the plant that provides water and nutrient uptake and plant anchorage to soil [8]. WCR adults will feed on above ground tissues, which will impact plant photosynthesis. WCRs have developed resistance to most insecticides classes [9][10][11][12][13], crop rotation [14], and commercial maize hybrids [15][16][17][18][19] because of their ability to rapidly evolve, showing the necessity to understand maize defense mechanisms at various levels for future maize WCR resistance engineering.
Several factors contribute to plant resistance to chewing insects. For example, maize resistance to WCR can be influenced by root architecture [20], biomechanical strength [21], Plants 2022, 11, 2335 2 of 19 or biochemical composition [22]. Among direct defense, plants can adapt their morphology by increasing density of trichomes or wax composition, which can create a physical barrier [22][23][24][25]. Maize inbred lines that have been developed, such as Mp496, Mp704, and Mp708, that have a broad-based resistance to insect pests [26,27]. Both Mp704 and Mp708 have Mp496 in their background, and Mp708 was the result of the crossing of Mp704 and Tx601, which is an inbred maize line susceptible to fall armyworm (Spodoptera frugiperda Smith) [28].
Phytohormones also play a major role in modulating plant defenses in response to herbivore feeding [29]. For example, jasmonic acid (JA) acts as a major phytohormone in providing enhanced resistance to Spodoptera exigua (beet armyworm) in maize and tomato [30]. The crosstalk between JA and salicylic acid (SA) in plants in response to insect herbivory is well documented [31]. For example, chewing insects feeding on host plants activate the SA pathway, which in turn suppresses the JA-dependent defenses [32]. In the caterpillar-resistant Mp708 maize genotype, both JA and ethylene (ET) pathways were required for the accumulation of Maize Insect Resistance 1-Cysteine Protease (Mir1-CP), a key defensive protein in defense against insect pests in maize [33,34]. However, the ET pathway, uncoupled from the JA pathway, can regulate the mir1 expression in providing enhanced defense against the phloem sap-sucking aphids [35]. Furthermore, in response to herbivory, plants activate several defense-related proteins, for example, proteases, protease inhibitors or peroxidases, which could negatively impact insect growth. Resistance to both above and belowground feeding caterpillars have been demonstrated in Mp708 plants [35][36][37][38][39][40][41]. Among the molecular mechanisms deployed by insect-resistant maize inbred lines Mp708 and Mp704, the protein Mir1-CP, as mentioned above, has been characterized as a key component for resistance to pests in maize [35,36,38,42,43]. Accumulation of Mir1-CP has been identified in maize in response to FAW and WCR infestation [38,39,44] and can damage the peritrophic matrix of caterpillars, which significantly reduces caterpillar growth [36,45]. In Mp708 maize genotype, mir1 was expressed in root tissues at 2, 4 and 7 days after WCR infestation [38]. Belowground feeding by WCR also resulted in mir1 transcript expression in the maize whorl tissues of Mp708, but not in susceptible Tx601 plants, indicating that WCR feeding activated defense mechanisms in Mp708 plants. Here, an RNA-seq approach was used to study the interactions between maize and WCR to investigate maize defense signaling networks from below to aboveground tissues and subsequent WCR compensatory transcriptomic responses.

The Maize Transcriptomic Response Varies from Tissues and Time Points
To characterize the impact of WCR feeding on maize Mp708 genotype, both leaves and roots tissues were collected at 1 and 5 days post infestation (dpi) and before feeding (control; 0 dpi). Among the 39,498 annotated protein-coding genes, 35,032 (88.8%) were expressed in at least one condition. On average, 31,570 genes were expressed in each condition, ranging from 30,282 genes (Leaves: L; 0 dpi) to 32,723 (Roots: R; 0 dpi). A principal component analysis (PCA) of the 35,032 genes expressed in at least one condition was performed ( Figure 1). PC1 accounting for 25.2% of the variation separated the transcriptomes by tissues (Figure 1), and PC2 accounting for 6.5% of the variance indicated a separation between 5 dpi time points and other time conditions in both leaves and roots ( Figure 1). As expected, the transcriptional profile of roots was different from that of leaves.  35,032 genes expressed in at least one of the conditions. Conditions are represented with colors (Mp708 leaves: red, Mp708 roots: blue) and shapes (circle: 0 dpi, triangle: 1 dpi and square: 5 dpi). Each shape indicates a replicate.

Differentially Expressed Genes Vary within Tissues and Have Minimal Overlap
Differentially expressed genes (DEGs) (|log2 (FC)| ≥ 1 and q value ≤ 5%) across comparisons and time points are given in Supplemental Table S1. A total of 17,383 DEGs were found for the comparisons: 0 vs. 1 dpi leaves, 0 vs. 5 dpi leaves, 1 vs. 5 dpi leaves, 0 vs. 1 dpi roots, 0 vs. 5 dpi roots, 1 vs. 5 dpi roots, 0 dpi roots vs. 0 dpi leaves, 1 dpi roots vs. 1 dpi leaves and 5 dpi roots vs. 5 dpi leaves. The least numbers of DEGs (<500) were found at 1 dpi leaves and roots (Figure 2A,B). However, more DEGs were found in leaves as compared than roots, indicating WCR herbivory of roots induced a huge transcriptional response in the leaves (Figure 2A,B). At 5 dpi, a substantial number of genes were downregulated in both leaves and roots (Figure 2A,B). A similar observation was made for the comparison of DEGs for 1 vs. 5 dpi leaves and roots ( Figure 2C). These results were suggestive of significant transcriptional remodeling of the leaf (aboveground) and root (belowground) tissue transcriptomes between 1 and 5 dpi.  35,032 genes expressed in at least one of the conditions. Conditions are represented with colors (Mp708 leaves: red, Mp708 roots: blue) and shapes (circle: 0 dpi, triangle: 1 dpi and square: 5 dpi). Each shape indicates a replicate.

Differentially Expressed Genes Vary within Tissues and Have Minimal Overlap
Differentially expressed genes (DEGs) (|log 2 (FC)| ≥ 1 and q value ≤ 5%) across comparisons and time points are given in Supplemental Table S1. A total of 17,383 DEGs were found for the comparisons: 0 vs. 1 dpi leaves, 0 vs. 5 dpi leaves, 1 vs. 5 dpi leaves, 0 vs. 1 dpi roots, 0 vs. 5 dpi roots, 1 vs. 5 dpi roots, 0 dpi roots vs. 0 dpi leaves, 1 dpi roots vs. 1 dpi leaves and 5 dpi roots vs. 5 dpi leaves. The least numbers of DEGs (<500) were found at 1 dpi leaves and roots (Figure 2A,B). However, more DEGs were found in leaves as compared than roots, indicating WCR herbivory of roots induced a huge transcriptional response in the leaves (Figure 2A,B). At 5 dpi, a substantial number of genes were downregulated in both leaves and roots (Figure 2A,B). A similar observation was made for the comparison of DEGs for 1 vs. 5 dpi leaves and roots ( Figure 2C). These results were suggestive of significant transcriptional remodeling of the leaf (aboveground) and root (belowground) tissue transcriptomes between 1 and 5 dpi. The overlap of DEGs for each comparison was next analyzed for each tissue ( Figure  3). A total of 434 DEGs and 1326 DEGs were upregulated in leaves at 1 and 5 dpi, respectively. Of these, 235 DEGs were shared in common between the two time points and represent 54.1% of the total number of DEGs at 1 dpi and 17.7% of the total number of DEGs at 5dpi ( Figure 3A). A gene function enrichment of the 235 upregulated DEGs identified molecular functions related to: DNA-binding transcription factor activity, DNA binding or calcium ion binding (Supplemental Table S2). One hundred fifty-three genes were downregulated at 1 and 5 dpi (leaves), which represent 57.5% of the DEGs at 1 dpi (n = 266 genes) and 6.4% of the DEGs at 5 dpi (n = 2380) ( Figure 3A). The gene function enrichment analysis characterized functions as: oxidoreductase activity, heme binding or peroxidase activity (Supplemental Table S2). The overlap of DEGs for each comparison was next analyzed for each tissue ( Figure 3). A total of 434 DEGs and 1326 DEGs were upregulated in leaves at 1 and 5 dpi, respectively. Of these, 235 DEGs were shared in common between the two time points and represent 54.1% of the total number of DEGs at 1 dpi and 17.7% of the total number of DEGs at 5dpi ( Figure 3A). A gene function enrichment of the 235 upregulated DEGs identified molecular functions related to: DNA-binding transcription factor activity, DNA binding or calcium ion binding (Supplemental Table S2). One hundred fifty-three genes were downregulated at 1 and 5 dpi (leaves), which represent 57.5% of the DEGs at 1 dpi (n = 266 genes) and 6.4% of the DEGs at 5 dpi (n = 2380) ( Figure 3A). The gene function enrichment analysis characterized functions as: oxidoreductase activity, heme binding or peroxidase activity (Supplemental Table S2).
For the root samples, 39 and 31 genes overlapped between the root conditions 1 and 5 dpi, upregulated or downregulated, respectively ( Figure 3B). The low number of genes did not allow us to perform a detailed gene function enrichment analysis. However, we examined the functions of the DEGs. Among the 39 DEGs, gene functions were related to cytochrome P450, SAUR-like auxin-responsive protein family or terpenoid cyclases/protein prenyltransferases superfamily protein (Supplemental Table S1). Among the 31 DEGs downregulated between 1 and 5 dpi, functions were related to jasmonate-zimdomain protein 11 or transcription factors (TFs) such as ERF, AP2 or G2-like ( Figure 3B) (Supplemental Table S1). We also monitored the overlap of DEGs between tissues at 1 or 5 dpi ( Figure 3C,D). At 1 dpi, a low number of up or downregulated genes overlapped between the two tissues. Twenty genes were upregulated in leaves and downregulated in roots ( Figure 3C). At 5 dpi, a higher number of genes overlapped for upregulated (448 genes) and downregulated (715) genes between root and leaf tissues ( Figure 3D). For the root samples, 39 and 31 genes overlapped between the root conditions 1 and 5 dpi, upregulated or downregulated, respectively ( Figure 3B). The low number of genes did not allow us to perform a detailed gene function enrichment analysis. However, we examined the functions of the DEGs. Among the 39 DEGs, gene functions were related to cytochrome P450, SAUR-like auxin-responsive protein family or terpenoid cyclases/protein prenyltransferases superfamily protein (Supplemental Table S1). Among the 31 DEGs downregulated between 1 and 5 dpi, functions were related to jasmonate-zim-domain protein 11 or transcription factors (TFs) such as ERF, AP2 or G2-like ( Figure 3B) (Supplemental Table S1). We also monitored the overlap of DEGs between tissues at 1 or 5 dpi ( Figure 3C,D). At 1 dpi, a low number of up or downregulated genes overlapped between the two tissues. Twenty genes were upregulated in leaves and downregulated in roots ( Figure 3C). At 5 dpi, a higher number of genes overlapped for upregulated (448 genes) and downregulated (715) genes between root and leaf tissues ( Figure 3D).

Genes Were Co-Expressed in Tissues and/or Time Points
A weighted gene co-expression network analysis (WGCNA) of the 17,383 non-redundant DEGs was performed, and genes were sorted into 10 modules: module 1 (M1) through M10 ( Figure 4) followed by gene function enrichment analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) (Supplemental Table S3). These modules could be separated by time and tissues.
A weighted gene co-expression network analysis (WGCNA) of the 17,383 non-redundant DEGs was performed, and genes were sorted into 10 modules: module 1 (M1) through M10 ( Figure 4) followed by gene function enrichment analysis using Kyoto Encyclopedia of Genes and Genomes (KEGG) (Supplemental Table S3). These modules could be separated by time and tissues.   Figure 4). In M1, genes were highly expressed at 0 and 1 dpi, with lower expression at 5 dpi; however, M8 (2910 genes) had the opposite gene expression profile with higher expression of genes at 5 dpi ( Figure 4). The 6273 genes included in M1 had functions related to protein localization, post-transcriptional gene silencing by RNA and positive regulation of phosphorylation. For M8, the 2910 genes were involved in regulation of the developmental process, regulation of RNA metabolite process, and cellular response to phosphate starvation processes were enriched (Supplemental Table S3). Genes that were part of M3 (481 genes) were expressed in roots at 0 and 1 dpi only ( Figure 4). Gene function enrichment showed functions for M3 related to ribosome, metabolic pathways, biosynthesis of secondary metabolites or phenylpropanoid biosynthesis, all with potential roles in defense. Genes found in M4 (168 genes) were expressed only in 1 and 5 dpi leaves and were associated with regulation of cellular process, response to auxin, or regulation of nitrogen compound metabolic process ( Figure 4). M5 (3612 genes) and M6 (3052 genes) were composed by genes expressed in leaves with descendant or ascendant expression profiles along the time for M5 and M6, respectively ( Figure 4). Genes in M5 were enriched in functions related to biosynthesis of secondary metabolites, metabolic pathways, fatty acid biosynthesis/metabolism/elongation or porphyrin and chlorophyll metabolism, indicating downregulation of these important leaf developmental processes under WCR-induced stress. M6 DEGs were associated with the regulation of cellular process, acyl-CoA metabolic process, glycerol ether metabolic process (Supplemental Table S3). KEGG enrichment analysis revealed that genes part of M6 were enriched in functions related to circadian rhythm, photosynthesis, metabolic pathways, phenylalanine metabolism or starch and sucrose metabolism, and fatty acid biosynthesis/metabolism/elongation, porphyrin, and chlorophyll metabolism (Supplemental Table S3).

Time Points: M2, M7, M9, M10
M2 comprised 261 genes expressed at 0 and 1 dpi in both tissues ( Figure 4). Biological process analysis showed gene enrichment related to translation synthesis, βglucan metabolism/biosynthesis process, spindle checkpoint, NADPH regeneration for M2 (Supplemental Table S3). Genes that are part of M7 (222 genes) and M9 (397 genes) were upregulated in both leaves and roots at 5 dpi ( Figure 4). Genes comprising these modules have functions related to cell communication and defense mechanisms. M7 was composed of genes with functions related to the regulation of catalytic activity, regulation of protein metabolic process, regulation of cell communication, regulation of phosphorus metabolic process, and regulation of response to stimulus, which indicate that both tissues respond to WCR attack by expressing genes involved in defense mechanisms. Genes in M9 had functions related to the regulation of RNA metabolic process, regulation of nitrogen compound, biosynthetic process, phosphorus metabolic process, programmed cell death, sucrose metabolic process. M10 had genes upregulated in uninfested roots and leaves at 1 dpi, with functions related to nucleoside phosphate biosynthetic process, monovalent inorganic cation transport or ATP biosynthetic process.

Modulations in Maize Hormone Pathways during WCR Infestation
The transcriptomic activity of genes associated with hormones involved in plant defenses [46] was investigated ( Figure 5).

Ethylene (ET)
There are five majors steps for the synthesis of ET [47]. Nineteen genes involved in ET pathway were differentially expressed between at least two conditions and could be divided into four gene families. SAMS (S-adenosyl-L-methionine synthase) genes were downregulated within the leaf contrasts but upregulated in 1 and 5 dpi roots compared to leaves ( Figure 5). ACS (ACC synthase) genes were downregulated for all the comparisons, except for one gene Zm00001d039487, which was upregulated in two contrasts. Among the ACO (ACC oxidase) genes, six of the eight DEGs were downregulated in uninfested roots compared to uninfested leaves, and within leaf comparisons ( Figure 5). Significant changes in ACO gene abundances were seen in the 1 and 5 dpi root comparisons. These differences were not observed in 1 dpi versus 0 dpi root comparisons. When compared to leaves, six and two ACO genes were down and upregulated, respectively, in the 1 dpi comparisons, and five and two down-and upregulated, respectively, in the 5 dpi comparisons ( Figure 5). In Arabidopsis thaliana, EIN2 is required for ethylene signaling [48]. Here, four genes of the EIN family were DEGs and generally downregulated across all contrasts, except for Zm00001d003451, which was significantly enriched in roots at 5 dpi compared to its control ( Figure 5).

Ethylene (ET)
There are five majors steps for the synthesis of ET [47]. Nineteen genes involved in ET pathway were differentially expressed between at least two conditions and could be divided into four gene families. SAMS (S-adenosyl-L-methionine synthase) genes were downregulated within the leaf contrasts but upregulated in 1 and 5 dpi roots compared to leaves ( Figure 5). ACS (ACC synthase) genes were downregulated for all the comparisons, except for one gene Zm00001d039487, which was upregulated in two contrasts. Among the ACO (ACC oxidase) genes, six of the eight DEGs were downregulated in uninfested roots compared to uninfested leaves, and within leaf comparisons ( Figure 5). Significant changes in ACO gene abundances were seen in the 1 and 5 dpi root comparisons. These differences were not observed in 1 dpi versus 0 dpi root comparisons. When compared to leaves, six and two ACO genes were down and upregulated, respectively, in the

Jasmonic Acid (JA)
α-linolenic acid (18:3) is the starting point for JA biosynthesis and acts as a substrate for lipoxygenases (LOX). The LOX gene family in maize contains 13 loci (ZmLOX1-13) [49]. Further, LOX enzymes were subdivided into two groups, namely 9-lipoxygenases and 13-lipoxygenases, depending on where they oxygenate α-linolenic acid. There are seven 9-lipoxygenases (ZmLOX1, 2, 3, 4, 5, 6, 12) and six 13-lipoxygenases (ZmLOX7, 8,9,10,11,13) in maize. 13-lipoxygenases catalyze the first step to the production of JA, while products of 9-lipoxygenases can still have defensive functions against insect herbivory [50]. Products of LOX catalysis are progressively converted by allene oxide synthase (AOS), allene oxide cyclase (AOC), and oxophytodienoic acid reductase (OPR), ultimately resulting in JA production [51,52]. Genes encoding AOS, LOX, and JAR1 (which conjugate JA to isoleucine [53]) were downregulated between all comparisons, except for transcripts for three LOX genes that were significantly enriched in the 5 dpi root/leaf comparisons ( Figure 5). Genes encoding OPR were downregulated in leaves across all comparisons, and transcripts of two OPR genes that were enriched in the control root/leaf comparisons remained enriched in the 1 and 5 dpi root versus leaves contrast. Notably, two other OPR genes were upregulated at 5 days in roots relative to the expression in 1-day roots. Maize OPR7 and OPR8 function in the conversion of OPDA to JA [54]. Additionally, it was shown that the Mp708 plants had elevated constitutive levels of JA [40,55]. Our results also demonstrate that the Mp708 roots had constitutive expression of OPR7 and OPR8, and these levels remained elevated and unchanged throughout the 5-day infestation of WCR (Supplemental Table S1). These results align with the previous findings that the Mp708 plants have constitutive accumulation of JA.

Salicylic Acid (SA)
We evaluated the genes involved in SA biosynthesis, as SA is often implicated in defense against pests. Biosynthesis of SA can take place via the isochorismate (IC) and/or phenylalanine ammonia-lyase (PAL) pathways [56,57]. In the IC pathway, IC produced in chloroplasts from chorismate by ICS is eventually converted to SA in the cytosol, whereas chorismate is exported from the chloroplasts to the cytosol and is used for the biosynthesis of phenylalanine. Phenylalanine is the primary precursor for the phenylpropanoid pathway and is converted to 4-cinnamic acid by PAL (phenylalanine ammonia lyase). 4-cinnamic acid is subsequently converted to benzoic acid, and ultimately to SA [54]. SA can be methylated by SAMTs to form methyl salicylate or conjugated to sugars by glycosyl transferases to form SA-sugar derivatives that are stored in vacuoles. Glycosidases can convert the SA-glycosyl esters back to free SA. Methylesterases (MES1) can convert SA methyl ester back to free SA and play a role of increasing cellular SA levels. SAMT1 was generally downregulated except in the 5 versus 1 dpi root comparisons. One copy of ICS2 was part of the DEGs and was downregulated in all comparisons ( Figure 5). Four copies of putative methylesterases, MES1, were found, of which three were downregulated in all comparisons, and one Zm00001d042209 appeared to be more enriched in roots as compared to leaves. One PAL gene appeared to be more enriched in roots relative to leaves, whereas all the other copies were downregulated across comparisons ( Figure 5). Similarly, many of the potential SA glycosyltransferases (SAGT1) were downregulated across comparisons, except for two that showed a mild enrichment in the roots. A UGT gene that was part of the DEGs was also downregulated in all comparisons ( Figure 5).

Maize DIMBOA Pathway Was Turned Off after Extended WCR Feeding
Benzoxanoids are secondary metabolites that act as natural pesticides. Their biosynthesis involves nine enzymes to form a linear pathway leading to the storage of DIMBOA as glucoside conjugate. Here, we identified DEGs encoding proteins involved in the synthesis of DIMBOA and DIMBOA-glucoside [58][59][60] (Figure 6). Out of the five genes associated with DIMBOA biosynthesis, transcripts for four were more enriched on roots relative to leaves across all comparisons, and one was downregulated, suggesting greater DIMBOA biosynthetic capacity in roots ( Figure 6). DIMBOA is glycosylated by several glucosyl transferases (GTs) using UDPG as a substrate to produce DIMBOA-ß-D-glucoside and UDP. Among the 14 GTs that can catalyze this reaction [61], eight were differentially expressed. In roots, two genes (Zm0001d019254 and Zm00001d019259) were upregulated and one gene (Zm00001d019250) was downregulated for the comparison 0 vs. 5 dpi and 1 vs. 5 dpi (Figure 6). For the roots, three genes (Zm00001d029250, Zm00001d019254 and Zm00001d034692) were upregulated at 0 vs. 5 dpi and 1 vs. 5 dpi. Only one gene was upregulated at 0 vs. 1 dpi (Zm00001d019251) and downregulated at 1 vs. 5 dpi (Figure 6). glucosyl transferases (GTs) using UDPG as a substrate to produce DIMBOA-ß-D-glucoside and UDP. Among the 14 GTs that can catalyze this reaction [61], eight were differentially expressed. In roots, two genes (Zm0001d019254 and Zm00001d019259) were upregulated and one gene (Zm00001d019250) was downregulated for the comparison 0 vs. 5 dpi and 1 vs. 5 dpi (Figure 6). For the roots, three genes (Zm00001d029250, Zm00001d019254 and Zm00001d034692) were upregulated at 0 vs. 5 dpi and 1 vs. 5 dpi. Only one gene was upregulated at 0 vs. 1 dpi (Zm00001d019251) and downregulated at 1 vs. 5 dpi (Figure 6).

WCR Transcriptome Is Remodeled between 1 and 5 dpi
Complementary to the maize transcriptome, we also investigated the WCR transcriptome before infestation (0 h), and after infestation at 1 and 5 dpi. A PCA analysis of the 19,358 genes expressed in at least one condition showed that samples were separated according to the time course (PC1, 17%), indicating differential regulation of transcription following feeding (Figure 7). In total, 4609 genes were differentially expressed (Supplemental Table S4). At 1 and 5 dpi, a higher number of genes were downregulated ( Figure  8). Furthermore, the number of up-or downregulated genes between conditions 1 and 5 dpi were comparable ( Figure 8A). We investigated the function of the up-and downreg-

WCR Transcriptome Is Remodeled between 1 and 5 dpi
Complementary to the maize transcriptome, we also investigated the WCR transcriptome before infestation (0 h), and after infestation at 1 and 5 dpi. A PCA analysis of the 19,358 genes expressed in at least one condition showed that samples were separated according to the time course (PC1, 17%), indicating differential regulation of transcription following feeding (Figure 7). In total, 4609 genes were differentially expressed (Supplemental Table S4). At 1 and 5 dpi, a higher number of genes were downregulated ( Figure 8). Furthermore, the number of up-or downregulated genes between conditions 1 and 5 dpi were comparable ( Figure 8A). We investigated the function of the up-and downregulated WCR genes for each comparison (Supplemental Table S4). At 1 dpi, the 1297 downregulated genes had functions related to defense response, response to biotic stimulus, response to other organisms, and response to external stimulus. In contrast, the 895 upregulated genes had functions related to the monocarboxylic acid metabolic process, cellular metabolic compound salvage, fatty acid metabolic process, cellular compound organization, or cytoskeleton organization (Supplemental Table S5). At 5 dpi, 2154 downregulated genes had functions linked to proton transmembrane transport, translation, peptide biosynthetic process, organonitrogen compound biosynthetic process, cation transport or defense response. We also functionally characterized the DEGs between conditions 1 and 5 dpi. The 707 downregulated genes at 5 dpi had functions related to defense response, response to stimulus signal transduction, cell communication, while genes upregulated at 5 dpi were associated with functions as DNA replication, DNA metabolic process, post-transcriptional gene silencing by RNA. peptide biosynthetic process, organonitrogen compound biosynthetic process, cation transport or defense response. We also functionally characterized the DEGs between conditions 1 and 5 dpi. The 707 downregulated genes at 5 dpi had functions related to defense response, response to stimulus signal transduction, cell communication, while genes upregulated at 5 dpi were associated with functions as DNA replication, DNA metabolic process, post-transcriptional gene silencing by RNA. Eighteen genes were upregulated at 1 dpi compared to 0 h (FPKM0dpi = 0). These genes were related to functions classified as bromodomain-containing protein, phosphatidylinositol 4-kinase alpha-like, Krueppel-like factor 12, epidermal growth factor receptor substrate 15-like 1, asparagine synthetase domain-containing (Supplemental Table S5). Genes upregulated at 1 dpi compared to 0 h and expressed in control condition (877 genes) were Eighteen genes were upregulated at 1 dpi compared to 0 h (FPKM 0dpi = 0). These genes were related to functions classified as bromodomain-containing protein, phosphatidylinositol 4-kinase alpha-like, Krueppel-like factor 12, epidermal growth factor receptor substrate 15-like 1, asparagine synthetase domain-containing (Supplemental Table S5). Genes upregulated at 1 dpi compared to 0 h and expressed in control condition (877 genes) were related to functions: calphotin-like, endoglucanase-like, chaoptin-like, venom peptide HsVx1-like or cuticle protein LPCP-23-like (Supplemental Table S5). Downregulated genes were associated with functions related to acaloleptin A-like, glycine-rich protein-like, cecropin-B2-like or acidic mammalian chitinase-like.

Discussion
Plant metabolism is impacted in response to biotic stress, and intraplant communication plays a crucial role in this process. We have previously shown that maize plants subjected to leaf herbivory by the European corn borer (ECB; Ostrinia nubilalis) were more resistant to WCR larval feeding on the roots [41]. In addition, the analysis of the maize root transcriptomes of plants infested with ECB indicated an upregulation of genes with functions related to phytohormone or defense in root tissues. Here, we show that the WCR

0 h vs. 5 dpi
Twelve DEGs were found upregulated at 5 dpi, but not expressed at 1 dpi or 0 h (FPKM 0dpi = 0 and FPKM 1dpi = 0). The functions of these genes were linked to: ctenidin-1-like, cuticle protein LPCP-23-like, neuropeptide-like 3, glucose dehydrogenase, neuropeptide-like 3, putative gustatory receptor 28 b (Supplemental Table S5). Genes expressed at 5 dpi but not expressed in control (1759 DEGs) have functions related to: bromodomain-containing protein, neuropeptide-like 3, prisilkin-39-like, phosphatidylinositol 4-kinase alpha-like, elongation of long-chain fatty acids protein, cuticle protein or proline-rich protein (Supplemental Table S5). Functions of downregulated genes (2154 DEGs) were associated with: alpha-crystallin B chain-like, thaumatin-like protein 1 b, tenecin-1-like, venom serine protease inhibitor-like, cathepsin L-like proteinase, putative glucosylceramidase 3, acaloleptin A-like, drosomycinlike or venom allergen 3-like (Supplemental Table S5). In addition, due to the overlap of the DEGs, a large number of genes up-or downregulated were shared between the comparisons: 0 vs. 1 dpi, 0 vs. 5 dpi or 1 vs. 5 dpi (Figure 8B). At 5 dpi, 782 genes were downregulated, whereas 777 genes were downregulated at 5 and 1 dpi when compared with the controls (0 dpi). Further, the gene function enrichment analysis identified functions related to cellular component organization, actin cytoskeleton organization, protein transport, or lipid biosynthetic process ( Figure 8B, Supplemental Table S5). Overlapping genes downregulated or upregulated for the three comparisons indicate that these genes are differentially expressed during the experimental time course. Here, 242 genes were downregulated and 160 genes were upregulated within the three comparisons 0 vs. 1 dpi, 0 vs. 5 dpi and 1 vs. 5 dpi ( Figure 8B). Gene function enrichment analysis showed that the 242 downregulated genes were associated with defense response, innate immune response, reproductive behavior, mating behavior, or visual perception (Supplemental Table S5). The 160 upregulated genes have functions related to DNA replication or DNA metabolic process. Three hundred and forty six genes were downregulated between 1 vs. 5 dpi and 0 vs. 5 dpi, suggesting that these genes were not impacted at 1 dpi ( Figure 8B). These genes were involved in proteolysis, neuropeptide signaling pathway, defense response, drug catabolic process, and hormone metabolic process.

Discussion
Plant metabolism is impacted in response to biotic stress, and intraplant communication plays a crucial role in this process. We have previously shown that maize plants subjected to leaf herbivory by the European corn borer (ECB; Ostrinia nubilalis) were more resistant to WCR larval feeding on the roots [41]. In addition, the analysis of the maize root transcriptomes of plants infested with ECB indicated an upregulation of genes with functions related to phytohormone or defense in root tissues. Here, we show that the WCR attack remodels the transcriptomes of both below and aboveground maize tissues. In addition, we highlighted the changes of the WCR transcriptome during the time-course of this study.
A previous study reported no overlap of DEGs between roots and leaves tissues after pest infestation [63]; here, we observed an overlap between up-and downregulated genes between leaves and roots at 5 dpi in WCR-infested plants relative to 0 dpi. Genotype differences and age of WCRs used in our study and that of Erb et al. [63] could account for this disparity. Nevertheless, our results indicated a strong transcriptomic response with common gene activation in below and aboveground tissues with increasing time of WCR feeding on maize Mp708 roots. For both leaves and roots tissues, we found a higher transcriptional response at 5 dpi compared to other time points indicating that response to WCR feeding built up over time.
Maize Mp708 line was defined as a resistant line in response to insect feeding associated with the expression of mir1. Here, mir1 (Zm00001d036542) was expressed in both roots and leaves tissues, but the expression of mir1 was constitutive in the leaves at all the time points (0, 1, and 5 dpi), whereas mir1 was upregulated in the root at 5 dpi compared to other time points (Supplemental Figure S1). This result suggested that maize defense mechanisms were activated after the initiation of WCR feeding and is not an immediate effect. The accumulation of mir1 transcripts was previously linked with JA and ET levels in Mp708 in response to fall armyworms [34]. Here, we showed that genes involved in the ET pathway are upregulated in root tissues, and their expression increased over time, indicating that ET plays an important role in defense response in Mp708 against WCR herbivory. Both ET and JA are well documented for their defensive roles during plant responses to insect herbivory [34,35,47,64,65]. Here, we have shown that genes related to SA pathway were mainly downregulated. SA and JA were previously described with an antagonistic action, but not all the time. In Arabidopsis, it was also demonstrated that both SA and JA act synergistically to promote effector-triggered immunity against Pseudomonas syringae pv. maculicola [66]. Because Mp708 plants had elevated levels of JA [40,55], higher constitutive expression of OPR7, a gene involved in the conversion of OPDA to JA, remained high before and after WCR feeding for 7 days [38]. Similarly, our results also demonstrated elevated constitutive expression of OPR7 and OPR8 genes in Mp708 roots before and after WCR feeding for 5 days (Supplemental Table S1). However, a recent study by Ye et al. (2022) [67], showed that the JA biosynthesis genes were upregulated in maize roots after WCR feeding for 72 h. The constitutive and induced expression of JA biosynthesis genes in the maize genotypes used in our study and in Ye et al. [67] may contribute to the discrepancy between results in maize-WCR interactions. In the future, a metabolomic approach would benefit to validate our transcriptomic results.
In addition to phytohormones, plants also use secondary metabolites as a defense mechanism. Among these metabolites, maize plants store benzoxazinoids in a non-toxic form when they are not injured, but tissue injury results in benzoxazinoids breaking down into toxic compounds [68]. However, these benzoxazinoids can be redirected and controlled by the herbivores and can be used as self-defense for WCR against their natural predators [69]. In roots, we found that the DIMBOA-synthesis genes were downregulated, while genes associated with DIMBOA-glucosyl transferase and DIMBOA-glucoside dehydrogenase were mostly upregulated. Collectively, our results suggest that the expression of genes related to benzoxazinoids in WCR resistant maize plant were altered after WCR herbivory.
The transcriptomic analysis of Mp708 revealed that after 1 day post infestation, expression profiles of genes with functions related to cell wall organization, cellulose metabolic process, or biosynthetic process were unchanged. However, their expression was drastically decreased after 5 dpi, suggesting that it takes at least 5 days post-WCR infestation to start affecting plant tissues and physiology. It is plausible that the newly hatched WCR larvae did not actively feed on the maize roots at the beginning (1 dpi), but as the days progressed (5 dpi), the late instar WCR larvae were able to modulate physiological changes in the plant tissues. Alternatively, it is equally likely that the "cues" released from late instar WCR larvae were able to suppress the maize physiological responses. Similarly, the transcription profile of mir1 revealed no change in the leaf, but the expression in the roots increased significantly at 5 dpi. These data suggest that there is a fine tuning of intraplant defense, providing resources to regions under active herbivory (roots) and not to distal sites (leaves).
Similar to the WCR herbivory affecting plant transcriptomes, WCR transcriptomes were also impacted after feeding for 1 and 5 days on maize Mp708 plants. Changes in the WCR transcriptome were broadly related to their development and counteracting plant defense mechanisms. At 1 dpi, upregulated genes compared to before infestation had functions related to larval growth and contained many uncharacterized genes. A large number of WCR DEGs with unknown functions highlight the need for further annotation of the WCR genome. As mentioned before, plant defense mechanisms were activated in maize roots at 5 dpi. This activation was also reflected in the WCR transcriptome with the suppression of the gene transcription associated with functions related to allergen or defense. Some of these WCR genes that were previously reported to be downregulated in the WCR larval transcriptome in response to two maize toxins Cry34/35Ab1 [70,71] were also downregulated in our study. These genes included zinc ion binding (downregulated at 1 and 5 dpi), lipase activity (downregulated at 1 dpi and upregulated at 5 dpi), lipase activity (upregulated at 5 dpi and 1 dpi) or upregulated cell communication (upregulated 0 dpi vs. 5 dpi) [72][73][74]. These results demonstrate that the plant defenses in a resistant maize genotype can negatively impact all the developmental stages of WCR.

Insect Growth
WCR eggs were obtained from Crop Characteristics, Inc., Farmington, MN, and were maintained as described previously [41]. Briefly, the WCRs used in this study were originally collected from susceptible non-Bt maize fields in Minnesota, United States. This colony was further maintained on non-Bt maize plants at Crop Characteristics Inc. for several generations (>80). The WCR eggs obtained were kept in growth chamber with 14:10 (L:D) h photoperiod at 23 • C for hatching. Newly hatched neonate larvae were used for the experiments.

Plant Growth and WCR Infestations
Mp708 plants were grown in 3.8 cm × 21.0 cm plastic Cone-tainers (Hummert International, Earth City, MO, USA) containing soil mixed with vermiculite and perlite (PRO-MIX BX BIOFUNGICIDE + MYCORRHIZAE, Premier Tech Horticulture Ltd., Olds, AB, Canada) in growth chambers with 14:10 (L:D) h photoperiod, 160 µE m −2 s −1 , 25 • C, and 50-60% relative humidity. All plants for the experiments were used at the V2-V3 developmental stage (~2 weeks) [75]. Maize roots were infested with ten newly hatched neonate WCR larvae [41] and plant tissues (both roots and leaves) were collected at 1 and 5 day post infestation (dpi). For leaf sample collection, the tissues in the yellow-green region of the whorl (where the caterpillars normally feed) were harvested. For the root tissue collection, 10 cm from the root-shoot junction, which forms the middle region of the entire root system, was harvested, as described previously [76]. WCR uninfested plants were used as controls. In addition, WCR used to infest the maize plants were collected at 1 and 5 dpi and before infestation (i.e., 0 dpi, which are the newly hatched larvae that were never exposed to Mp708 roots). In total, three replicates were collected for each experimental condition. For each replicate, three infested or uninfested leaves and root tissues and WCR were collected and flash-frozen in liquid nitrogen.

RNA Extraction and RNA-Seq Libraries Construction and Sequencing
Maize and WCR tissues (80-100 mg) were ground using 2010 Geno/Grinder ® (SPEX SamplePrep, Metuchen, NJ, USA) for 40 s at 1400 strokes min −1 . Total RNA was extracted from the homogenized tissue using Qiagen RNeasy Plant Mini Kit. Extracted total RNA was quantified through Nanodrop 2000c Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). Then, RNA-seq libraries were constructed by using the mRNA-seq standard TruSeq protocol from Illumina. RNA-seq libraries were sequenced to produce 50 bp pairedend reads. Each sample had an average of 20 million reads. The transcriptomics dataset is available under Bioproject: PRJNA781637.

Functional Annotation
The GOBU package was used for enrichment calculations [82]. The full set of maize gene annotation was used as the reference comparison set against down-or upregulated DEGs. The p values were calculated using Fisher's exact test and were corrected for multiple testing with FDR method using the R module called 'p-adjust'.

Conclusions
The transcriptomic analysis of maize above and belowground tissues after WCR infestation indicated that there are changes in the gene expression both at the site of WCR feeding and tissues distal to WCR feeding (i.e., leaves). Plant transcriptomic responses modulated over time and were higher at 5 than at 1 dpi. Based on our data and previous findings, it is highly likely that the resistance mechanism in the Mp708 genotype is a multitrait phenotype, and Mir1-CP could act as a central defensive protein involved in imparting pest resistance. Determining the role of Mir1-CP-dependent protein networks or other defense proteins connected to Mir1-CP that contribute to above-and belowground defense signaling is critical to understand global and tissue-specific defenses in the insect-resistant maize genotype. Additionally, we also observed changes to the WCR transcriptome with the deactivation of defense genes after 5 dpi. Understanding the major maize defense signaling networks and compensatory responses in WCR to combat plant defenses could lead to the development of novel pest management strategies.
Supplementary Materials: The following supporting information can be downloaded at: https: //www.mdpi.com/article/10.3390/plants11182335/s1, Supplemental Table S1: Contrast summary for the genes expressed in maize roots and root tissues at 1 dpi and 5 dpi. Supplemental Table S2: Gene function enrichment output from Venn Diagram. Supplemental Table S3: Gene function enrichment output for each module. Supplemental Table S4: Contrast summary for the genes expressed in WCR tissues at 1 dpi and 5 dpi. Supplemental Table S5: Gene function enrichment for WCR. Supplemental Figure S1: The expression level of mir1 (Zm00001d036542) across each maize conditions. Author Contributions: J.L. conceived and designed the research. L.P. performed the computational analysis. S.B. performed the experiments. G.S. and N.N.V. contributed to method development and reagents. W.P.W. developed the maize genotype used in the study. L.P., G.S. and J.L. wrote the article with contributions from all authors. All authors have read and agreed to the published version of the manuscript.
Funding: This work was supported at different times by the Nebraska Agricultural Experiment Station with funding from the Hatch Act (accession #1007272) through the USDA National Institute of Food and Agriculture and funds from the University of Nebraska-Lincoln to J.L.

Data Availability Statement:
The transcriptomics dataset is available under Bioproject: PRJNA781637.