Seeds of maize inbred lines B73, Zong3, Chang7-2, A188, and Mo17 were sown in 12 cm-diameter plastic pots filled with humus soil and vermiculite (about 7 to 1 ratio). Plants were grown in a greenhouse at ~17°C (night) and ~28°C (day) under a 14 h photoperiod until the fourth leaf had fully expanded (V4, 3-week-old seedlings). A R. padi colony was established from a single first instar nymph collected from a field in the Kunming Institute of Botany, Chinese Academy of Sciences, in Kunming, China, in 2014. The insect colony was maintained in an 80 × 80 × 150 cm cage on sorghum plants at 17–28°C, 14 h day length.

To evaluate the survival rates of R. padi, 20 apterous adults were transferred to the middle of the forth leaves and enclosed in clip cages (3 cm diameter), and after 24 h, the adult aphids were brushed off and 20 neonates were left in the clip cages to evaluate every day in the next 5 days (20 replicated maize plants for each group) (Carena and Glogoza, 2004; Mao and Zeng, 2014). To determine days from birth to first reproduction, life span, and progeny per aphid, one new neonate was placed on the middle of the forth leaf of each plant and the performance parameters were estimated until its death (20 replicated maize plants for each group) (Kettles et al., 2013).

For treatment on maize, aphid feeding and mock treatments were performed as the following: for each aphid-treated plant, 50 adult aphids were gently placed onto the middle of the fourth leaf and enclosed in a clip cage, and the plants of the control group (mock) were fixed with empty clip cages. Three replicates were used for the aphid feeding and control group. In order to minimize the variations caused by circadian rhythm and other differences in environmental conditions, we treated leaves at different time points prior to the collection time. Aphids were quickly removed with a soft paintbrush at the time of harvest, and the regions in the cages were immediately excised and snap-frozen in liquid nitrogen.

50 adult aphids were feeding on plant for 0 and 48 h, five biological replicates were used, and a previously established method was employed (Qi et al., 2016). In brief, samples were extracted with 80% methanol and then analyzed on an Agilent 1200 Rapid Resolution Liquid Chromatography (RRLC) system equipped with a ZORBAX SB-Aq column (2.1 × 100 mm, 1.8 μm) followed by an Agilent 6510 Q-TOF running in positive ionization mode. The column temperature was set at 50°C, and the flow rate was 0.3 mL min⁻¹. All metabolites were identified based on pure standards, except some Bxs: DIMBOA was confirmed by comparing the Bx profiles between maize inbred line H88 and bx1 mutant (GRMZM2G085381); DIMBOA-Glc, DIM₂BOA-Glc, HDMBOA-Glc, and MBOA were confirmed with purified standards, and the other putative Bxs were identified by their molecular masses. The abundance of each metabolite was estimated from the peak area. Peak detection and matching were performed using bioconductor XCMS and CAMERA packages (http://www.bioconductor.org/.) Samples (metabolites after integration) with P < 0.05 were considered to be differentially regulated. The acquired data were introduced to the SIMCA-P software package (v11.5, Umetric, Umea, Sweden) for Principal Component Analysis (PCA). The data was filtered by orthogonal signal correction (OSC) to remove variations from non-correlated factors.

Three biological replicates were used for RNA-seq analyses. Total RNA was extracted from ground leaf samples using TRIzol reagent (Invitrogen). The concentration of total RNA was measured using a NanoDrop 2000 Spectrophotometer (Thermo Scientific). Sequencing was performed at 3.75 G depth on a HiSeq2500-PE125 platform (Illumina). The raw sequence data are deposited in NCBI's Short Read Archive under the BioProject ID (PRJNA338511). The resulting sequences were trimmed based on quality scores and mapped to the maize B73 reference genome sequence V2 and maize working gene set V5a with Tophat2 (Kim et al., 2013) using the following modifications from the default parameters: maximum intron size 100,000; minimum intron size 20; up to two mismatches allowed (Trapnell et al., 2012). The expression levels of genes were estimated using Cufflink and were normalized using the numbers of reads per kb of exon sequence in a gene per million mapped reads. Differential expression analysis was performed using the Cuffdiff (Trapnell et al., 2012) with a cutoff of 2-fold change relative to the control levels.

All AS analyses were based on splice junctions obtained from the BAM files produced by TOPHAT2. For differential splicing analysis, we used Cuffdiff to categorize individual transcripts based on their transcription start site (Sammeth et al., 2008), thus grouping all the transcripts with a common pre-mRNA molecule. Thereafter, we performed Cuffdiff to estimate significant expression differences of each transcript using the Jensen-Shannon divergence metric across the different samples (Trapnell et al., 2012).

For gene ontology (GO) analysis, we obtained the GO terms for each maize gene using the AgriGO web service (Du et al., 2010), and WEGO (Ye et al., 2006) was used for GO functional annotation; Plant MetGenMAP (Joung et al., 2009) was used for pathway enrichment analysis. Sets of genes showing a similar pattern of transcription were assumed to be functionally correlated. Transcription patterns were clustered using MeV 4.9.0 (Saeed et al., 2003) with Euclidean distances and the hierarchical cluster method of complete linkage clustering and Java software. MeV 4.9.0 was also used to create heat-map visualizations of differentially expressed genes. To identify the transcription factors, we referenced a maize transcription factor database retrieved from the Plant Transcription Factor Database 4.0 (Jin et al., 2015, 2017).

To evaluate the levels of resistance to R. padi in different maize inbreed lines, we measured the survival rates of R. padi on the fourth fully expanded leaves of the inbred lines B73, Zong3, Chang7-2, A188, and Mo17 for 6 days (Figure 1A). R. padi showed the highest survival rate on B73 (86.3% survived on day 1), while exhibited the lowest survival rate on Mo17 (69% survived on day 1). The survival rates decreased in B73 and Mo17 and the trend continued in the next 5 days (65.6 and 10.5% aphids survived on B73 and Mo17 respectively on day 6) (Figure 1A). Congruently, the life spans of R. padi on B73 and Mo17 were 22.9 and 18.6 days, respectively, and the new neonates on B73 needed 8.7 days to start to give offspring, while 11 days were needed on Mo17 (Figure 1B). R. padi fecundity on B73 was much higher than on Mo17 (39.7 and 12.9 offspring per aphid on B73 and Mo17, respectively) (Figure 1C). Thus, B73 and Mo17 are a R. padi-susceptible and a-resistant line, respectively. These results were similar to previous findings, where Mo17 exhibited a much greater defense level against the maize leaf aphid R. maidis than did B73 (Meihls et al., 2013; Betsiashvili et al., 2015).

To explore the metabolites in B73 and Mo17 that might account for the observed R. padi resistance, leaf samples were collected 48 h after R. padi infestation. The profiles of Bx and phenolic metabolites were obtained on an UPLC-Q-TOF MS system (Table S1).

Given the importance of Bxs in maize defense against insects (Long et al., 1977; Rostas, 2007; Niemeyer, 2009; Glauser et al., 2011), we quantified the contents of Bxs (Figure 2A, Table S1). In general, Bxs were much higher in Mo17 than in B73. For example, DIMBOA, DIM₂BOA, DIMBOA-Glc, and DIM₂BOA-Glc in Mo17 were at least 2-fold > B73 (Figure 2A, Table S1). Notably, R. padi feeding did not have an obvious impact on the contents of Bxs in either line. Another group of metabolites that we examined were the phenolics. The contents of phenolics also showed large differences between B73 and Mo17. Without R. padi treatment, 10 phenolic compounds, apigenin O-hexoside, apigenin C-hexoside C-pentoside, quercetin O-hexoside, quercetin O-deoxyhexosyl-hexoside, luteolin, luteolin O-hexoside, luteolin O-hexoside O-deoxyhexoside, neochlorogenic acid, isorhamnetin O-hexoside, and isorhamnetin O-rhamnoglycoside were 0.64–3.29-fold higher in B73 than in Mo17, while the content of scytalone was 34.9% lower in B73 (Figure 2B, Table S1). After 48 h of R. padi feeding, most phenolics decreased in B73 (39–77%), but remained constant in Mo17.

We used RNA-seq technology to analyze the genome-wide gene expression responses of B73 and Mo17 to R. padi feeding at 6 and 24 h.

Differentially expressed genes (DEGs) were selected based on significance (P < 0.05) (Table S2): 3,789 and 3,420 DEGs were found in B73 after 6 and 24 h R. padi feeding, and the same treatment resulted in 4,127 and 1,786 DEGs in Mo17. Gene functions were further classified using GO analysis. As shown in Figure S1 (Table S2): 5,657 DEGs in B73 and 4,683 DEGs in Mo17 were assigned to 25 GO terms (Cellular Component, Molecular Function, and Biological Process: 8, 6, and 11 GO terms, respectively). The GO terms in B73 and Mo17 were largely similar except that some from Cellular Component (Extracellular Region Part and Macromolecular Complex) and Biological Process (such as Biological adhesion and Cellular Component Organization) showed significant differences.

Detailed expression profiling, based on the criteria that the genes should be up- or down-regulated more than 1-fold and P < 0.05, indicated that both B73 and Mo17 responded to R. padi infestation, but in a different manner (Table S3). After 6 h of R. padi feeding (Figure 3A, Table S3), 745 (319 up- and 426 down-regulated) and 853 (555 up- and 298 down-regulated) genes showed changes in expression levels in B73 and Mo17, respectively, and among these, only 162 were regulated in both lines. After 24 h of R. padi infestation (Figure 3B, Table S3), still 567 DEGs were uniquely regulated in B73, while Mo17 had specific 221 DEGs; only 80 genes showed altered levels in both B73 and Mo17 at 24 h. Thus, B73 and Mo17 have very different transcriptomic responses to R. padi infestation.

The DEGs were screened with a cluster analysis based on the k-means method using Pearson's correlation distance (TM4 software, http://www.tm4.org/). Each cluster represents genes whose expressions showed similar response patterns to R. padi infestation (Figure 4A, Table S4). In B73, 1,157 genes were classified into three expression cluster groups (Clusters 1–3). B73 Cluster 1, which contained 338 genes, showed a positive slope from 0 to 6 h, followed by a drop in the expression level at 24 h. The most abundant group was the B73 Cluster 2, with 628 genes showing a negative slope from 0 to 6 h, and then returned to the basal levels at the 24 h time point. B73 Cluster 3 had 191 genes whose expression showed a positive slope from 0 to 24 h. In Mo17, 958 transcripts were clustered into two groups (Figure 4A, Table S4). Cluster 1 contained 617 genes whose expression showed a positive slope from 0 to 6 h, followed by a drop in the expression level at 24 h, and Cluster 2 contained 341 genes whose expression showed an initial drop from 0 to 6 h, followed by reversion to the basal levels at 24 h (Figure 4A, Table S4). To reveal the different biological processes in each gene expression cluster, pathway enrichment analysis was performed using MetGenMAP (Joung et al., 2009). Pathways that were enriched among the transcripts in R. padi-treated B73 included (1) many plant defense- and stress response-related genes (Figure 4B, Table S4), such as benzoate degradation II in Cluster 1; (2) proline degradation, tryptophan biosynthesis, and chlorogenic acid biosynthesis I in Cluster 2; (3) photorespiration, Rubisco shunt, glycine biosynthesis, phospholipid desaturation, and glycolipid desaturation, which were mostly related to primary metabolism in Cluster 3. Notably, defense-related salicylate biosynthesis and chlorogenic acid biosynthesis could be found in Cluster 1, and DIMBOA-glucoside degradation was in Cluster 2 (Figure 4B, Table S4).

In order to gain further insight into the regulatory mechanism underlying the differences between the responses of B73 and Mo17 to R. padi on the gene expression level, differentially regulated TFs were identified. It was found that R. padi herbivory specifically regulated 35 and 14 TFs in B73 and Mo17, respectively, and 9 TFs were regulated in common (Figure 5A, Table S5). We found that the 44 TFs distributed in seven families (MYB, WRKY, b HLH, NAC, GRAS, HSF, and Dof) in B73, most of which were down-regulated (Figure 5B, Table S5). In contrast, the 23 TFs in Mo17 distributed in five families (MYB, WRKY, b HLH, NAC, and GRAS), and most of them were up-regulated (Figure 5C, Table S5). However, 7 out of the 9 commonly regulated TFs in B73 and Mo17 were down-regulated in B73 but were up-regulated in Mo17 (Figure 5D, Table S5).

Given the large gene expression differences between B73 and Mo17 in response to R. padi feeding, we expected that these two lines may also show a distinct transcript pattern at individual time points.

We compared the DEGs between B73 and Mo17 under the control condition and after 6 and 24 h of R. padi feeding, and it was found that 2,174 genes were always differentially expressed under all these conditions, but specifically, 1,437, 496, and 748 genes exhibited different expression levels only under control conditions and after 6 and 24 h of R. padi feeding, respectively (Figure 6A, Table S6).

Specifically, (1) the DEGs between B73 0 and Mo17 0 h were clustered to 11 pathways, most of which were associated with primary metabolism, except salicylate biosynthesis and chlorogenic acid biosynthesis I (Figure 6B, Table S6); (2) DEGs between B73 and Mo17 at 6 h were linked to four defense-related pathways, which were phenylpropanoid biosynthesis, chlorogenic acid biosynthesis I, tyrosine biosynthesis, and phenylalanine biosynthesis; (3) DEGs between B73 and Mo17 at 24 h were enriched on five pathways: Proline biosynthesis I, Abscisic acid biosynthesis, Proline biosynthesis III, Arginine biosynthesis IV, and Citrulline biosynthesis (Figure 6B, Table S6). Thus, B73 and Mo17 showed different response to R. padi feeding in gene expression pattern. More detailed transcriptomic comparisons between B73 and Mo17 over time are shown in Table S6.

We next determined the transcript levels of TFs in B73 and Mo17 under control condition, 6 and 24 h after R. padi treatment (Figure 7, Table S7). Despite the large differences between the transcriptomes of B73 and Mo17, only 96 TFs were found to have different expression levels between B73 and Mo17 under the control condition, while 62 and 57 TFs were differentially expressed between these two lines after being treated with aphid feeding for 6 and 24 h, respectively (Figure 7A, Table S7); and among these, 38, 12, and 11 TF genes were specifically different only in the control samples and those induced by 6 and 24 h of R. padi herbivory. Only 28 TFs were differentially expressed between B73 and Mo17 at all time points, and these TFs showed the same patterns of regulation (Figure 7B, Table S7). The entirety of TFs that were differentially expressed between B73 and Mo17 were distributed in WRKY, MYB, NAC, b HLH, C2H2, GRAS, HSF, Dof, and b ZIP (Figure 7C, Table S7). Under all conditions, most TFs of the WRKY, MYB, NAC, b HLH, and C2H2 group had greater expression levels in B73 than in Mo17. Notably, the number of WRKY, MYB, and NAC TFs, whose transcript levels were higher in B73 than in Mo17, decreased after R. padi feeding; for example, in B73 0 h samples, 19 WRKYs had higher expressions than in Mo17, but in B73 6 and B73 24 h the number of TFs that had higher expression in B73 decreased to 9 and 12, respectively (Figure 7C). It is likely that TFs may have played an important role in shaping the differences of B73 and Mo17 transcriptomes under control and R. padi-induced conditions.

In order to understand the possible genetic basis for the different levels of Bxs in B73 and Mo17, the transcriptome data were screened for Bx biosynthesis-related genes that had been annotated (http://www.maizegdb.org/), and the newly reported BX13 and BX14 (Handrick et al. (2016). Most of the Bx biosynthetic genes had different expression levels between B73 and Mo17 (Table S7). Importantly, BX1 and BX13 levels were much higher in Mo17 than in B73 (at least 7.8 and 33 times greater in Mo17 than in B73) in all samples (Figures 8A,B, Table S7). BX1 is very important for Bx biosynthesis (Butron et al., 2010); without the expression of BX13, DIM₂BOA-Glc cannot be detected, suggesting BX13 plays a significant role in Bx biosynthesis (Handrick et al., 2016). As Bx levels, such as DIM₂BOA-Glc and HDMBOA-Glc, were much higher in Mo17 than in B73 (Figure 2A), BX1 and BX13 likely accounted for the high Bx contents in Mo17.

Only four TFs, a MYB (GRMZM2G108959), a GRAS (GRMZM2G015080), a NAC (GRMZM2G179049), and a WRKY (GRMZM2G425430), showed greater expression levels in Mo17 than in B73 under all conditions (Figure 8C, Table S7), suggesting that these four TFs may play a role in regulating the expression of BX1 and BX13, resulting in the high Bx levels in Mo17.

Terpene synthases (TPSs) catalyze the formation of diverse volatile terpenes (Bohlmann et al., 1998; Cheng et al., 2007; Dudareva et al., 2013). To investigate whether the volatile terpene biosynthesis pathways respond to R. padi feeding, we screened for the expression levels of 11 TPS genes in B73 and Mo17 that had been annotated by the Plant MetGenMAP to be 10 sesquiterpene synthases (TPS1 to 10) and one monoterpene synthase (GRMZM2G030583) (Joung et al., 2009). Among these 11 TPS genes inspected, TPS8 (GRMZM2G038153) was induced (1.93- and 1.11-fold increased at 6 and 24 h, respectively) by R. padi infestation while TPS10 was repressed (84 and 43% decreased after 6 and 24 h R. padi feeding, respectively) in B73, but not in Mo17 (Figure S2A, Table S8). TPS2 (GRMZM2G046615) and TPS3 (GRMZM2G064406) were 4.98- and 3.91-fold up-regulated by R. padi herbivory after 6 h R. padi feeding in Mo17, but not in B73 (Figure S2B, Table S8). We further compared the expression level of TPS genes between B73 and Mo17. Without R. padi feeding, there were five genes (TPS2, TPS3, TPS7, TPS8, and TPS10) differentially expressed between these two lines: TPS2, TPS3, TPS8, and TPS10 showed higher levels in B73 compared to those in Mo17, while TPS7 showed the opposite pattern (Figure S3A, Table S8). After 6 h of R. padi feeding, only TPS7 and TPS8 showed different expression levels between B73 and Mo17 (Figure S3B, Table S8), and after 24 h of R. padi feeding TPS7 still expressed higher in B73 while the expression level of TPS10 was lower in B73 than in Mo17 (Figure S3C, Table S8).

The differentially spliced genes (DSGs) were investigated to determine whether AS might be involved in the response of B73 and Mo17 to R. padi feeding. 6 h after R. padi feeding on B73, 36 DSGs were identified, and there was no overlap between DEGs and DSGs at this time point (Figure 9A, Tables S3, S9); only one gene (among 59 DSGs) was both a DEG and a DSG after 24 h of R. padi feeding on B73 (Figure 9B, Tables S3, S9); among these 36 and 59 DSGs, only two genes showed changes in AS in B73 at both 6 and 24 h after R. padi feeding (Figure 9C, Table S9). Similarly, 59 and 39 DSGs were found in Mo17 after 6 and 24 h of R. padi feeding, respectively (Figure 9F, Table S9), and only a few genes were both regulated on the expression and AS level (Figures 9D,E, Tables S3, S9). Furthermore, six genes were found to be common between the DSGs identified in the samples of Mo17 6 and Mo17 24 h (Figure 9F, Table S9). To gain further insight into R. padi-induced dynamic AS responses in B73 and Mo17, pathway enrichment analysis was performed on the DSGs using MetGenMAP (Joung et al., 2009). No grouping of the DSGs in B73 6 and B73 24 h (Figure S4, Table S9) into enriched pathways showed no overlapping function. Notably, three pathways were in common among the enriched pathways in the DSGs of Mo17 6 and Mo17 24 h: synthesis of suberin, salicylate biosynthesis, and phenylpropanoid biosynthesis, and these are all defense-related (Figure S4, Table S9). These results indicate that B73 and Mo17 respond specifically to R. padi feeding on the AS level.

In order to gain further insight into the differences between these two maize lines, the DEGs and DSGs between B73 and Mo17 under control condition and at 6 and 24 h after R. padi treatment were determined. Consistent with the large differences in gene expression levels between these lines, 688 DSGs were identified between B73 0 and Mo17 0 h, and 159 genes were regulated both on expression and AS level (Figure 10A, Tables S6, S10). Similarly, after 6 and 24 h of R. padi infestation, B73 and Mo17 had 595 and 658 DSGs, respectively, and 99 and 122 genes were both DEGs and DSGs (Figures 10B,C, Tables S6, S10). The enriched pathways in the DSGs between B73 and Mo17 under control and after 6 h of R. padi feeding included flavin, spermine, and spermidine biosynthesis; terpenoid biosynthesis is the most enriched pathway in the DSGs identified between B73 and Mo17 after 24 h of R. padi herbivory (Figure S5, Table S10). Comparison of all these three groups of DSGs indicated that there were 181, 126, and 166 condition-specific DSGs between B73 and Mo17, while 340 genes were found consistently differentially spliced under all conditions (Figure 10D, Table S10). Furthermore, 10 TFs were found among these 340 DSGs, which belonged to the families of WRKY, bZIP, b HLH, GRAS, MYB, C2H2, MYB, and GATA (Table S10), and notably, none of these TFs showed differences in their expression levels between these two maize lines.

These results suggest that R. padi feeding induces specific AS patterns in different maize lines. Furthermore, AS difference also existed between B73 and Mo17 and is largely independent of the regulation of gene expression levels.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.