Differential Contributions of MYCs to Insect Defense Reveals Flavonoids Alleviating Growth Inhibition Caused by Wounding in Arabidopsis

In Arabidopsis, basic helix–loop–helix transcription factors (TFs) MYC2, MYC3, and MYC4 are involved in many biological processes, such as defense against insects. We found that despite functional redundancy, MYC-related mutants displayed different resistance to cotton bollworm (Helicoverpa armigera). To screen out the most likely genes involved in defense against insects, we analyzed the correlation of gene expression with cotton bollworm resistance in wild-type (WT) and MYC-related mutants. In total, the expression of 94 genes in untreated plants and 545 genes in wounded plants were strongly correlated with insect resistance, and these genes were defined as MGAIs (MYC-related genes against insects). MYC3 had the greatest impact on the total expression of MGAIs. Gene ontology (GO) analysis revealed that besides the biosynthesis pathway of glucosinolates (GLSs), MGAIs, which are well-known defense compounds, were also enriched in flavonoid biosynthesis. Moreover, MYC3 dominantly affected the gene expression of flavonoid biosynthesis. Weighted gene co-expression network analysis (WGCNA) revealed that AAE18, which is involved in activating auxin precursor 2,4-dichlorophenoxybutyric acid (2,4-DB) and two other auxin response genes, was highly co-expressed with flavonoid biosynthesis genes. With wounding treatment, the WT plants exhibited better growth performance than chalcone synthase (CHS), which was defective in flavonoid biosynthesis. The data demonstrated dominant contributions of MYC3 to cotton bollworm resistance and imply that flavonoids might alleviate the growth inhibition caused by wounding in Arabidopsis.


INTRODUCTION
Because of their sessile lifestyle, plants have developed numerous strategies to cope with biotic and abiotic stresses from the environment, such as insect herbivores and pathogenic infections (Howe and Jander, 2008). Jasmonates (JAs) are the main regulators involved in defense response and plant development (Erb and Glauser, 2010;Okada et al., 2015;Chen et al., 2019;Wasternack and Strnad, 2019;Zander et al., 2020). The biosynthesis and signal transduction of JAs have been extensively studied in the past decades (Browse, 2009;Erb et al., 2012;Wasternack and Hause, 2013;Howe et al., 2018).
MYC transcription factors (TFs) that have a conserved basic helix-loop-helix (bHLH) motif, are the main regulators of the JA signaling pathway. Although MYC2, MYC3, and MYC4 have redundant functions, evidence has shown their differential contributions in Arabidopsis (Carretero-Paulet et al., 2010;Goossens et al., 2017;Wang et al., 2019). MYC2 is considered as the master regulator of most aspects of JA signaling transduction (Kazan and Manners, 2013;Yang et al., 2019), whereas MYC3 and MYC4 have specificity similar to that of MYC2 to bind the cis-element G-box, which assists in the activation of a gene expression (Fernandez-Calvo et al., 2011).
MYC TFs are well-known regulators of plant defense against insects. The triple mutant myc2 myc3 myc4 (mycT) is hypersensitive to the generalist herbivore Spodoptera littoralis (Fernandez-Calvo et al., 2011;Schweizer et al., 2013a). Of the redundant target genes regulated by MYCs, the direct defense genes and the differential contributions of MYC2, MYC3, and MYC4 in plant defense against herbivorous insects have not been systemically studied. Here, we found that myc-related mutants have different resistance to generalist herbivores (Helicoverpa armigera). Using weighted gene coexpression network analysis (WGCNA), we screened out 626 genes, of which the expressions were highly correlated with insect resistance; and further, they were named as MGAIs (MYC-related genes against insects). MYC3 contributes more to MGAIs expression than MYC2 and MYC4. MGAIs were significantly enriched in GLS and flavonoid biosynthesis.
We further found that flavonoids alleviate the growth inhibition caused by wounding treatment in Arabidopsis. The data elaborated the differential contributions of MYCs to insect resistance and provided a possible role for flavonoids in plant defense.
The larvae of cotton bollworm (H. armigera) and diamondback moth (Plutella xylostella) were obtained from the Institute of Zoology, Chinese Academy of Sciences, and reared on an artificial diet at25 • C with 70% relative humidity and a 14-h light/10-h dark photoperiod. For insect feeding tests, second instar H. armigera larvae and P. xylostella larvae 2 days post hatching were reared on Arabidopsis leaves (replaced with fresh leaves every day) for 4 days and weighed.

Intermittent Wounding Treatment
After the plants were grown for 12 days, the leaves were wounded for the first time using microtips. The newly grown leaves were subsequently wounded in the same way every 3 days. After 3-4 rounds of the treatment, the aerial portions were harvested and used for biomass detection.

RNA Extraction, Reverse Transcription, and Quantitative Real-Time PCR
Total RNA was extracted from Arabidopsis using a TRIzol reagent (Invitrogen, Waltham, MA, United States). One microliter of DNase I (40 U/µl, Promega, Madison, WI, United States) was added to total RNA to remove genomic DNA, and1 µg of total RNA was used to prepare the first strand cDNA (Transgene, Beijing, China). Real-time polymerase chain reaction (qRT-PCR) was performed with a Bio-Rad CFX Connect qRT-PCR system (Bio-Rad, United States) using SYBR Green PCR Mix (Takara, Japan), according to the instructions of the manufacturer for the standard step amplification program. Arabidopsis S18 (At4g09800) was used as an internal standard. Biological triplicates were performed with technical duplicates. All nucleotide primers used in this study are listed in Supplementary Table 1.

Arabidopsis Sample Preparation and RNA Sequencing
For transcriptome sequencing, the second pair of true leaves of the 18-day wild-type (WT) (Col-0) and myc-related mutants (myc2-2, myc3, myc4, myc2/3, myc2/4, myc3/4, and mycT) were wounded and harvested 4 h post wounding. The corresponding untreated plants were used as controls (CK). Each sample contained three biological replicates. Total RNA was extracted and approximately 1 µg of total RNA was used to enrich poly(A) mRNA using oligo-dT magnetic beads (Invitrogen, Waltham, MA, United States), followed by fragmentation into 100-400 nt sizes; and the fragments were subsequently used to synthesize cDNA with random hexamer primers (Invitrogen, Waltham, MA, United States). RNA sequencing was performed on an Illumina HiSeqXten platform (Illumina, San Diego, CA, United States) at Majorbio (Shanghai, China).

Differential Gene Expression Analysis
Differentially expressed genes (DEGs) were determined by comparing expression levels in WT and myc-related mutants using the R package DESeq2 (Love et al., 2014). Only genes with relatively high expression levels (mean count > 1) were subjected to further analysis. Genes with a log2 fold change greater than 1 and adjusted P-value less than 0.05 in Bonferroni correction were defined as DEGs.

Correlation Calculation and MGAI Identification
We used the mean weights of cotton bollworms fed on mycrelated mutants to reflect the resistance of the corresponding plants. We calculated the Pearson correlation between cotton bollworm resistance (opposite numbers to average of weight of insects) and gene expression in different plants and tested the results by Student's t-test. We chose the genes with a correlation greater than 0.6 and a P value less than 0.05 as MGAIs.

Co-expressed Analysis
Weighted gene co-expression network analysis was performed to generate an unsigned co-expression network on all sample data based on topology overlap measurement (TOM) using the R package WGCNA (Langfelder and Horvath, 2008). Genes with an average TPM (transcripts per million) greater than one among all samples were used to generate the network. The soft threshold β was chosen to be 18. All the genes were clustered into 55 different modules (and a gray module containing genes that could not be classified to any module) by hierarchical clustering using a distance matrix (1-TOM).

Enrichment Analysis
All enrichment analyses were based on Fisher's exact test with Bonferroni correction. Gene ontology (GO) enrichment analysis was performed using GO annotations downloaded from TAIR 1 . GO enrichment, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment, and module enrichment analyses were performed using the R package cluster Profiler (Yu et al., 2012) and KEGGREST (Tenenbaum and Maintainer, 2020).

The myc-Related Mutants Display Differential Resistance to Cotton Bollworm
To clarify the differential contributions of MYC2, MYC3, and MYC4 in plant defense against herbivorous insects, we tested the resistance of the myc-related mutants to that of the generalist insect, H. armigera, and the specialist insect, P. xylostella. For H. armigera larvae, the single mutant myc3 and the myc3 related double mutants (myc2 myc3, myc3, and myc4) displayed attenuated resistance, and the attenuation was more serious in the triple mutant (myc2 myc3 myc4, and mycT); however, the resistance of the other single mutants (myc2 and myc4) and double mutant (myc2 and myc4) were similar to that of the WT ( Figure 1A). None of the single mutants (myc2, myc3, and myc4) showed any obvious change in resistance to P. xylostella larvae. All the double mutants displayed reduced resistance; however, no obvious differences among the double mutants were observed, and mycT displayed the highest sensitivity to P. xylostella larvae ( Figure 1B). These results indicated that MYC3 played a dominant role in resistance to cotton bollworm, while the three MYCs had highly redundant functions in resistance to P. xylostella.
Many plant genes are induced by wounding damage caused by chewing insects and are thought to be important in insect defense. We used wounding pre-treated plants (24 h post wounding) for insect feeding to test whether wounding pre-treatment could improve plant resistance to insects. When plants were pretreated with wounding, the resistance of WT but not mycT to H. armigera was enhanced but not mycT ( Figure 1C). However, the resistance to P. xylostella was barely affected by pretreatment with wounding in either WT or mycT ( Figure 1D). These results indicate that the enhanced resistance to H. armigera caused by wounding is regulated by MYCs. In panels (C,D), the plants were wounded 1 day before being used for feeding. Data are means ± SEM (n = 25-35) analyzed by one-way ANOVA and Tukey's honest significant difference (HSD) test. Different letters indicate significant differences, P-value < 0.05.

MYC2, MYC3, and MYC4 Regulate Different Biological Processes
To gain a deeper insight into the plant defense regulated by MYCs, RNA-seq data of the WT and myc-related mutants both with and without wounding were analyzed. Owing to the redundant function of MYCs, we compared double the mutants (myc2, myc3, myc2, myc4, myc3, and myc4) with mycT to analyze the effects of single MYC on gene regulation. Genes with significantly higher expressions in myc3 and myc4 than in mycT were defined as the MYC2-affected genes, and so on, to the MYC3-and MYC4-affected genes. Genes that were significantly higher in the WT than in mycT were defined as MYC-affected genes (Figure 2A).
Without wounding, the expression of 478 genes was significantly higher in WT than in mycT, and with the wounding treatment, 538 genes were found to have higher expression in WT than in mycT. Venn diagram analysis revealed 162 overlapping genes (Figure 2B and Supplementary Table 2). The 316 MYC-affected genes from the untreated samples were enriched in flavonoid and anthocyanin biosynthesis processes. The overlapping genes were enriched in "GLS biosynthetic process" and those only from the wounded sample analysis were mainly enriched in "response to jasmonic acid" (Figure 2C  and Supplementary Table 3). These results revealed a general profile of MYC-regulated genes. We selected six MYC-affected genes and analyzed their expression in different mutants by RT-qPCR, and the results are consistent to those of RNA-seq assay (Supplementary Figure 1).

Comparisons of gene expression in double mutants and
mycT both with and without wounding revealed 508 MYC2affected, 372 MYC3-affected, and 463 MYC4-affected genes (Supplementary Figure 2A and Supplementary Table 4). A total of 172 genes, whose expression was affected by all the three MYCs (Figure 2D), were significantly enriched in "GLS biological process" (Supplementary Figure 2B and Supplementary Table 5), suggesting functional redundancy of MYCs in this process. Genes only affected by MYC2 (132) were enriched in "response to jasmonic acid" and "response to auxin, " whereas genes only affected by MYC3 (143) were enriched in "flavonoids biosynthetic process" and "response to red light"; and genes only affected by MYC4 (132) were enriched in "trichome patterning" (Figure 2E and Supplementary Table 5). These results suggest that despite functional redundancy, MYCs differentially regulate gene expression in multiple biological processes.

MYC3 Has the Most Significant Impact on the Expression of Defense Genes Against Cotton Bollworm
In addition to defense genes, MYCs also regulate multiple genes involved in growth, development, etc. To screen for MYC-regulated genes in defense against cotton bollworm, we analyzed the correlation between gene expression and cotton bollworm resistance in WT and myc-related mutants. The expression of 94 genes in untreated plants and 545 genes in wounded plants were found to be strongly correlated with insect resistance (Pearson correlation coefficient, r > 0.6) (Supplementary Table 6), and these genes were termed as MGAIs (Supplementary Figure 3). Interestingly, the Venn diagram showed that only 13 MGAIs in untreated and wounded plants overlapped (Figure 3A and Supplementary Table 7). MGAIs from the untreated samples were significantly enriched in the "flavonoid biosynthesis process, " whereas those from the wounding treated samples were enriched in "GLS biosynthetic process" (Figure 3B and Supplementary Table 8), consistent with a previous study showing that GLSs are well-known antiinsect compounds (Wittstock and Gershenzon, 2002;Mewis et al., 2006;Schlaeppi et al., 2008). GLS synthetic gene expression was largely downregulated in mycT but was not significantly reduced in single or double mutants (Supplementary Figure 4A and Supplementary Table 9), which further confirmed the redundant function of the three MYCs in the regulation of GLS biosynthesis. RT-qPCR analysis of the selected GLS synthesis genes was consistent with the RNA-seq results (Supplementary Figure 4B). Furthermore, the GLS content in the plants was highly correlated with their resistance to the cotton bollworm (Supplementary Figures 4C,D).
Of the 626 total MGAIs, the expressions of 273 genes were remarkably affected by MYCs, for their expression in either untreated or wounding-treated mycT was significantly downregulated (Supplementary   Supplementary Table 11. Data are analyzed by one-way ANOVA and Tukey's HSD test. Different letters indicate significant differences, P-value < 0.05. were more than MYC2(114)-and MYC4(92)-affected genes ( Figure 3C, Supplementary Figure 5, and Supplementary  Table 10). For the 94 MGAIs in the untreated plants, the total expression levels were significantly reduced in mycT; and notably, of the three single mutants, the lowest expressions were observed in myc3 (Figure 3D and Supplementary Table 11). Similar trends in gene expression were observed for the 545 MGAIs in wounded plants ( Figure 3E and Supplementary Table 11). These results suggest the dominant function of MYC3 on the expression of MGAIs in either untreated or wounded plants.

Classification of Gene Modules Highly Correlated With Plant Resistance to Cotton Bollworm
Genes with highly correlated expression are helpful in determining their association with biological processes. A total of 16,660 genes with an average TPM value greater than one across all 48 tested groups were subjected to WGCNA. In total, 10,991 genes were divided into 55 modules (Supplementary Figure 6). The correlation between the module eigengenes and cotton bollworm resistance is presented in Figure 4A. The eigengenes of two modules (dark magenta and maroon) in untreated plants and of six modules (plum1, orange, floral white, dark green, dark orange 2, and green yellow) in the wounded plants were highly correlated with plant resistance (r > 0.6) ( Figure 4A). The eigengenes of only two modules (sky blue and dark slate blue) in the untreated and wounded plants showed a similar positive correlation with resistance.
In the untreated plants, genes whose expression was affected by the three MYCs were mainly distributed across eight modules, such as dark magenta, of which eigengenes are highly correlated with resistance to cotton bollworm (Figures 4A,B). The distributions of MYC2-, MYC3-, and MYC4-affected genes across the eight modules were quite different, reflecting their distinct role in regulating gene expression in untreated plants. For example, the genes classified as dark magenta were primarily affected by MYC3 (Figure 4B and Supplementary Table 12) and were enriched in the flavonoid biosynthetic process (Figure 4D and Supplementary Table 13).
In the wounded plants, seven modules, such as dark slate blue, were significantly affected by MYCs. Interestingly, the distributions of MYC2-, MYC3-, and MYC4-affected genes across these seven modules were similar, indicating that the functions of MYC2, MYC3, and MYC4 on gene expression in wounding response are highly redundant (Figure 4C and Supplementary  Table 12). The MYC-affected genes of dark orange 2, plum1, and floral white modules were all enriched in "the response to jasmonic acid" (Figure 4E and Supplementary Table 13), reflecting the dominant role of JA in the wounding response regulated by MYCs.
Furthermore, the total MYC-affected genes of dark slate blue in both untreated and wounded plants were enriched in the GLS biosynthetic process. The expression of dark slate blue in untreated and wounded plants was positively correlated with cotton bollworm resistance (Figure 4A), suggesting the important role of GLS in resistance against cotton bollworm in both untreated and wounded plants (Figures 4D,E).

Flavonoid Synthesis and Growth-Related Genes Are Highly Co-expressed and Mainly Regulated by MYC3
In Arabidopsis, 19 genes were involved in flavonoid synthesis according to the KEGG annotation, and the expression of most genes was positively associated with resistance to cotton bollworm in untreated plants ( Figure 5A). The expressions of most of these genes were found to be suppressed in all the  Supplementary Table 14. Data are analyzed by one-way ANOVA and Tukey's HSD test. Different letters indicate significant differences, P-value < 0.05. (C) Co-expression network of flavonoid synthesis and growth-related genes. Blue and green stand for flavonoids synthesis genes and growth-related genes, respectively. (D) Scatter plot displays the co-expression of AAE18 with the eight flavonoid synthesis genes. X-axis stands for the AAE18 expressions (TPM), and Y -axis stands for the expressions (TPM) of the indicated genes in flavonoids synthesis. Dots represent the gene expressions in untreated samples, and crosses represent the gene expressions in wounding treated samples. (E,F) Images (E) and Biomass (F) of WT and CHS with intermittent wounding treatments. When plant has grown for 12 days, the first wounding treatment is performed, and the wounding treatment is performed every 3 days for a total of three times. (E) Scale bar, 1 cm. (F) Data are means ± SEM (n = 18) and analyzed by two-way ANOVA and Tukey's HSD test. *P < 0.05, **P < 0.01, ****P < 0.0001.
MYC3-related mutants (mycT, myc3 myc4, and myc3) but not in myc2, myc4, and myc2 myc4 ( Figure 5A). Moreover, the total expression of flavonoid synthesis genes was significantly reduced in all the MYC3-related mutants. These results indicate that flavonoid synthesis genes were mainly regulated by MYC3 ( Figure 5B and Supplementary Table 14). In addition to flavonoid synthesis genes, we noticed that some growth-related genes (AAE18, PUR7, and CGA1) were also classified in the dark magenta module (Supplementary Figure 7), of which AAE18 (Acyl-activating enzyme 18) is involved in the metabolism of auxin precursors to active auxins (Wiszniewski et al., 2009). CGA1 (cytokinin-responsive GATA factor 1) and PUR7 (purin 7) are auxin response genes that are involved in plant growth (Senecoff et al., 1996;Zubo et al., 2018). These three genes were highly co-expressed with flavonoid synthesis genes ( Figure 5C). We selected eight flavonoid synthesis genes to individually analyze their co-expression with AAE18 using a scatter plot. The results, along with Pearson correlations, showed that the AAE18 expression was highly correlated with all the selected genes ( Figure 5D). This implies that flavonoids might act as linkers between growth and defense.

Flavonoids Might Alleviate the Growth Inhibition Caused by Intermittent Wounding in Arabidopsis
To decipher the role of flavonoids in plant defense, we tested a flavonoid synthesis mutant, chs, for resistance to cotton bollworm. Although the high concentration of flavonoids (kaempferol-3,7-dirhamnoside) in the artificial diet was toxic to Pieris brassicae (Onkokesung et al., 2014), chs displayed similar resistance to H. armigera, P. xylostella, and Spodoptera frugiperda, as the WT (Supplementary Figure 8A). In addition, the inhibition of root growth and the induction of gene expression by the JA treatment were similar in chs and WT ( Supplementary  Figures 8B-D), indicating that the JA response was not affected in chs.
A tradeoff between plant defense and growth has been reported (Major et al., 2017;Wasternack, 2017). Active defense responses usually inhibit plant growth. Some plants exhibit a strong compensatory capacity for growth by increasing the photosynthetic rate and compensatory regrowth, and improving the utilization rate of stored nutrients to compensate for the damage caused by stress (Stowe et al., 2000;Turley et al., 2013). The expression of flavonoid synthesis genes was highly correlated with that of growth-related genes ( Figure 5C), which inspired the authors to determine whether flavonoids are involved in growth regulation when plants are under biotic stress.
Herbivorous insect actions are intermittent rather than persistent. To mimic insect feeding, we treated WT and chs with intermittent wounding. The plant sizes of the WT and chs were similar under normal conditions. After 22 days of intermittent wounding treatment, both types of plants were smaller than the untreated control; and, notably, the chs plants were much smaller than those of WT ( Figure 5E and Supplementary  Figure 9). Consistently, the plant biomass was reduced by intermittent wounding, and the reduction was more significant in chs ( Figure 5F). These results showed that intermittent wounding caused more significant growth inhibition in chs than in WT. A plausible explanation is that flavonoids alleviated the growth inhibition caused by wounding in Arabidopsis.

DISCUSSION
Jasmonate signaling plays an essential role in coordinating plant defense and development. MYC2/MYC3/MYC4 are the main TFs and have functional redundancy in JA signaling, such as plant defense and floral transition (Fernandez-Calvo et al., 2011;Wang et al., 2017;Van Moerkercke et al., 2019). However, the GO enrichment analysis of genes affected only by MYC2, MYC3, or MYC4 revealed the differential functions of individual MYC in many biological processes ( Figure 2E). MYC2 might function predominantly in plant growth and development, and flavonoid synthesis is primarily regulated by MYC3, while MYC4 is involved in trichome formation.
We found 94 MGAIs in the untreated plants, suggesting their possible involvement in plant constitutive defense. In addition, the 545 MGAIs in the wounded plants indicated their potential role in wounding-induced defense. Interestingly, only 13 MGAIs in untreated and wounded plants overlapped (Figure 3A and Supplementary Table 7), suggesting that constitutive and wounding-induced defense against cotton bollworm followed different mechanisms in Arabidopsis. Among the three single mutants, the expression of MGAIs in myc3 was the lowest, and the cotton bollworm gained the highest weight when raised on myc3, indicating the dominant contribution of MYC3 to plant defense against cotton bollworm.
Phenylpropanoids are abundant phenolic secondary plant metabolites derived from the aromatic amino acid phenylalanine in plants. Inducible phenylpropanoids and their derivatives, such as flavonoids, are tightly coupled with many plant biotic stresses, namely, pathogen infections and insect herbivory (Odonbayar et al., 2016;Ranjan et al., 2019;Vanholme et al., 2019;Seybold et al., 2020). Accumulated flavonoids are believed to affect auxin transport and affect plant growth (Besseau et al., 2007;Dong et al., 2020). Furthermore, higher concentrations of flavonoids have been reported to inhibit the growth of some insects, such as whitefly (Onkokesung et al., 2014;Xia et al., 2021). Flavonoid accumulation could be easily induced by various abiotic or biotic stresses, such as ultraviolet B or temperature (Petrussa et al., 2013). This indicates that flavonoids might have important roles in plants for adaptation to stresses. Although in this study the lack of flavonoids had little effect on Arabidopsis resistance to insects, flavonoids might be beneficial in relieving the growth inhibition caused by wounding (Supplementary Figure 8A and Figure 5F). Coupled with the evidence that flavonoid synthesis genes are highly co-expressed with auxin response genes (Figures 5C,D), flavonoids may be implicated in coordinating auxin signaling to help plants recover from wounding or insect attack.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available. This data can be found here: All Illumina sequencing data have been deposited in National Center for Biotechnology Information (NCBI) database: Short Read Archives (SRA) database (BioProject ID PRJNA720476, Run ID 17877327-17877374).

AUTHOR CONTRIBUTIONS
D-DW, PL, and Y-BM designed the study, analyzed the data, and wrote the manuscript with the input of all the authors. D-DW performed most of the experiments with the assistance of Q-YC, X-YC, Z-WY, and M-YW. PL performed the sequence analysis. All authors contributed to the article and approved the submitted version.