The soybean (Glycine max L.) lines “Zhechun No. 3” (ZC3) and ZC-2 were used in this study. ZC-2 is a mutant that was identified among progeny derived from ZC3 seeds irradiated with 150 Gy gamma rays (Yuan et al., 2007). Seeds of each soybean line were sown in sterilized nutrient soil in plastic pots (35 cm diameter) and transferred to a pre-sterilized greenhouse. Three seedlings were retained in each pot. Twenty pots were prepared for each soybean line.

The C. truncatum isolate CT5 (Feng et al., 2021) was cultured on potato dextrose agar (PDA) at 25°C in the dark for 5 days. The mycelial suspension used as inoculum was prepared as follows. Ten mycelial disks (5 mm diameter) from an actively growing culture of isolate CT5 on PDA were added to potato dextrose broth (PDB) in 500 mL flasks and cultured for 4 days in an incubator at 25°C with shaking at 120 rpm in the dark. After incubation, the PDB was filtered through sterilized gauze. The mycelial pellets were rinsed with sterilized water at least six times to remove as much residual PDB as possible. The residual water was squeezed out and the pellets were weighed. The pellets were resuspended in sterilized water, fragmented in a blender at low speed for 20 s, and diluted with sterilized water to a concentration of 50 mg mL^–1. All steps were conducted on an ultra-clean workbench.

Healthy pods of uniform maturity (approximately 15 days after flowering) from each soybean line were sampled, quickly transported to the laboratory, and washed with sterile water. The pods were placed in a beaker containing 500 mL mycelial suspension, gently stirred for 10 s, and then removed from the suspension. After draining off the excess mycelial suspension, the inoculated pods were placed on filter paper in 15 cm petri dishes containing 3 mL sterilized water. As a control, pods were mock-inoculated using sterile water. The petri dishes were placed in an incubator at 25°C in the dark. Samples of inoculated and corresponding control pods for transcriptomic and metabolomic analysis were collected at 4, 8, 12, 24, and 48 h post-inoculation (HPI), immediately frozen in liquid nitrogen, and stored at −80°C. Three replicates for RNA-Seq and six replicates for UHPLC-MS were prepared. Each replicate comprised at least 15 pods. Symptoms of anthracnose on pods of soybean lines ZC-2 and ZC3 at 2 and 5 days after inoculation with Colletotrichum truncatum were recorded. Disease blotches were scanned, and their areas were calculated using the LA-S plant analysis system (Wseen, Shenzhen, China).

Total RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, United States) in accordance with the manufacturer’s protocol. The RNA quality was assessed using an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, United States) and RNase-free agarose gel electrophoresis. RNA library preparation and sequencing were performed by the Gene Denovo Biotechnology Co. (Guangzhou, China) on an Illumina HiSeq 2500 platform.

Raw reads were filtered with fastp v0.18.0 to obtain high-quality clean data. An index of the soybean reference genome (Wm82.a2) was constructed and paired-end clean reads were mapped using HISAT2.2.4 (Kim et al., 2015) with the default parameters. The fragments per kilobase of transcript per million mapped reads (FPKM) value was calculated to quantify the expression abundance and variation of genes using StringTie v1.3.1 (Pertea et al., 2015, 2016). Differential expression analysis was performed using DESeq2 and edgeR (Robinson et al., 2010; Love et al., 2014). Genes with false discovery rate < 0.05 and log₂ expression fold change ≥1 were considered to be differentially expressed genes (DEGs). Principal component analysis (PCA) was performed with the gmodels R package¹ to evaluate the relationships among samples and replicates. A statistical Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis for the DEGs was performed using KOBAS software (Kanehisa et al., 2008).

To determine the first sampling time point, the expression of nine disease-responsive genes was detected by qRT-PCR. In addition, 24 DEGs were selected for qRT-PCR analysis to validate the RNA-Seq data. The methods for qRT-PCR analysis and statistical tests were performed as described previously (Jin et al., 2021). The specific primers used are listed in Supplementary Table 1.

Determination of the relative content of terpenoids was performed using UHPLC-MS. Extraction and sample preparation were performed in accordance with a protocol described previously (Chen et al., 2013; Doppler et al., 2016). The UHPLC separation was performed on an ExionLC System (Sciex Technologies, Redwood City, CA, United States). The mobile phase A was 0.1% formic acid in water and the mobile phase B was acetonitrile. The column temperature was set at 40°C. The auto-sampler temperature was set at 4°C and the injection volume was 2 μL. A Sciex QTrap 6500 + LC-MS/MS system (Sciex Technologies, Redwood City, CA, United States) was used for assay development. Typical ion-source parameters were as follows: IonSpray voltage, +5500/−4500 V; curtain gas, 35 psi; temperature, 400°C; ion source gas 1, 60 psi; ion source gas 2, 60 psi; and declustering potential, ±100 V.

SCIEX Analyst Work Station v1.6.3 software (Sciex Technologies, Redwood City, CA, United States) was employed for multiple-reaction monitoring for data acquisition and processing. The MS raw data (.wiff) files were converted to TXT format using MSconventer. An in-house R program and database were used for peak detection and annotation. The SIMCA 16.0.2 software package (Sartorius Stedim Data Analytics AB, Umeå, Sweden) was used for orthogonal projections to latent structures–discriminate analysis (OPLS-DA). The variable importance in the projection (VIP) value of the first principal component in the OPLS-DA analysis was determined. The metabolites with VIP > 1 and p < 0.05 (Student’s t-test) were considered to be significantly differentially accumulated metabolites (DAMs). The integrated analysis of transcriptomic and metabolomic data was conducted on the basis of the enriched KEGG pathways.

Eight terpenoids or phytohormones (simvastatin, 11-keto-beta-boswellic acid, cantharidin, picrocrocin, 3-indolebutyric acid, abscisic acid, methyl jasmonate, and trans-zeatin riboside) were randomly selected and tested for their antifungal activity and ability to induce anthracnose resistance. PDA medium was supplemented individually with these standard biochemical reagents (final concentration 10 μM); standard PDA medium was used as the control. Each plate was inoculated with a C. truncatum mycelium plug (5 mm in diameter). Each treatment comprised 20 plates. After culture at 25°C in the dark for 6 days, the plates were scanned and the colony area was calculated with the LA-S analysis system (Wseen, Shenzhen, China).

Healthy pods of ZC3 of uniform maturity were sampled. The pods were sprayed with an aqueous solution (10 μM) of one of the standard reagents, sterile water was used as the control. After moisturizing at 25°C for 1 h, the pods were rinsed with sterile water and inoculated with C. truncatum, as described above. Fifty pods were treated for each treatment. After incubation for 5 days, the disease blotches on the pods were scanned and the disease area was calculated with the LA-S analysis system (Wseen, Shenzhen, China).

Symptoms of anthracnose on pods of soybean lines ZC-2 and ZC3 were recorded at 2 and 5 days after inoculation with C. truncatum. At 2 days after inoculation, tiny blotches had appeared on the pods of both lines (Figure 1A). However, fewer ZC-2 pods had developed blotches than ZC3 pods (Figure 1B). The proportion of the ZC-2 pod surface covered by blotches was significantly lower than that for ZC3 (Figure 1C). At 5 days after inoculation, distinct brown or black blotches had developed on the pods of ZC-2 and ZC3. The majority of lesions on the ZC-2 pods were tiny spots, whereas large lesions were common on the ZC3 pods (Figure 1A). The proportion of the ZC-2 pod surface covered by blotches was significantly lower than that for ZC3 (Figure 1D). These results indicated that the anthracnose resistance of ZC-2 was stronger than that of ZC3.

The expression analysis of disease-responsive genes by qRT-PCR showed that PR1, PR9, and WRKY genes were significantly up-regulated at 4 HPI in ZC-2 and/or in ZC3. However, only two PR9 genes were significantly up-regulated at 2 HPI in ZC-2 (Figure 1E). Therefore, we selected 4 HPI as the first sampling time point. At 48 HPI, tiny blotches had appeared on the pods (Figure 1A). To avoid contamination of soybean RNA-Seq data with RNAs from C. truncatum, due to extensive C. truncatum colonization of the pods, 48 HPI was selected as the last sampling time point. We included three intermediate sampling time points (8, 12, and 24 HPI).

Sixty RNA libraries were sequenced. More than 2.38 billion raw reads were obtained, comprising approximately 34.19–51.85 million raw reads for each library. After performing quality control on the data, 99.59–99.72% of the reads were retained in each library. We mapped 32.83–49.98 million reads to the reference genome with a mapping rate of 95.83–96.95%, of which 93.74–95.10% were uniquely mapped and 90.35–92.36% were mapped in exon regions (Supplementary Table 2). The expression levels for 56,725 genes were calculated using the FPKM method, which resulted in identification of 33,431 expressed genes (FPKM ≥ 1), including 684 novel genes (Supplementary Table 3).

The PCA showed that the three biological replicates for each sample were grouped together, which was indicative of high reproducibility, and thus that the experimental protocol and analysis methods were stable and reliable (Figure 2). All inoculated samples (enclosed by the green dashed ellipse in Figure 2) and all control samples (enclosed by the orange dashed ellipse in Figure 2) were grouped separately, and the inoculated samples were grouped to the upper right of the control samples. In addition, each inoculated sample was located to the upper right of the corresponding control sample (the corresponding samples are connected by a black line in Figure 2). These results indicated that the soybean pod transcriptome was reprogrammed in response to C. truncatum infection at each sampling time point, and this reprogramming showed a certain degree of consistency.

A total of 13,462 DEGs in the ten comparison groups were screened, of which 10,254 were detected in ZC-2 and 9,813 in ZC3 (Supplementary Table 4). The DEG statistics for each comparison group are shown in Figure 3A. From 4 to 48 HPI, the number of DEGs showed a U-shaped distribution in ZC-2 and ZC (Figure 3A). That is, many more DEGs were detected at 4 and 48 HPI than at the intervening time points, which indicated that ZC-2 and ZC3 showed two large-scale response periods to C. truncatum infection. At 4 HPI, the pod cells had detected C. truncatum infection and initiated the resistance response. At 48 HPI, a fresh large-scale resistance response may have been stimulated to respond to the progression of infection with the appearance of tiny lesions on the pods at approximately 48 HPI. In addition, a higher number of DEGs were detected at these two time points for ZC-2 than for ZC3. Except for the comparison group at 8 HPI in ZC3, each comparison group contained a higher frequency of up-regulated DEGs than down-regulated DEGs, especially in the comparisons for ZC-2 (Figure 3A).

Representation of these comparative results as Venn diagrams revealed that both unique and shared DEGs occurred between and among pairs in ZC-2 and ZC3, and the DEGs of ZC-2 and ZC3 showed general correspondence (Figure 3B). Overall, 6605 genes, which accounted for 66.4 and 67.3% of the total number of DEGs in ZC-2 and ZC3, respectively, were differentially expressed in both ZC-2 and ZC3 in response to C. truncatum infection. For both ZC-2 and ZC3, a majority of DEGs were shared DEGs, especially those from the comparison groups at 8, 12, and 24 HPI. The similar DEGs distribution and the high proportion of shared DEGs between ZC-2 and ZC3 implied that the response patterns of ZC-2 and ZC3 to C. truncatum infection were highly similar on account of their highly similar genetic background.

Twenty-four DEGs were selected to validate the RNA-Seq expression levels by qRT-PCR. The results showed a degree of reproducibility between transcript abundances assayed using RNA-Seq and the expression profiles revealed by qRT-PCR (Figure 4). However, inconsistency in abundances was suggestive of differences in sensitivity between the two methods.

To explore the biological functions of the DEGs and the resistance mechanisms of ZC-2 and ZC3 in response to C. truncatum infection, KEGG pathway enrichment analysis was performed. The enriched pathways in ZC-2 and ZC3 were roughly identical, indicating that the defense response mechanisms of the two lines to C. truncatum infection were highly similar (Figure 5). The enriched pathways were predominantly involved in amino acid metabolism, carbohydrate metabolism, terpenoid and polyketide metabolism, phenylpropanoid biosynthesis, flavonoid biosynthesis, isoflavonoid biosynthesis, plant–pathogen interaction, and signal transduction, all of which are important components of the resistance response mechanism. In particular, the enriched pathways that showed temporal continuity, such as “cyanoamino acid metabolism,” “phenylpropanoid biosynthesis,” “flavonoid biosynthesis,” “isoflavonoid biosynthesis,” “alpha-linolenic acid metabolism,” and “starch and sucrose metabolism,” which were highly enriched in ZC-2 and ZC3, may be vital for soybean resistance to C. truncatum infection. The pathways specifically enriched at 4 and 8 HPI in both ZC-2 and ZC3 were largely involved in amino acid and carbohydrate metabolism. The pathways “mitogen-activated protein kinase (MAPK) signaling pathway–plant,” “plant hormone signal transduction,” “plant–pathogen interaction,” and pathways involved in metabolism of terpenoids and polyketides were responsive to C. truncatum infection more promptly and/or more drastically in ZC-2 than in ZC3. Therefore, these pathways may be crucial factors in the enhanced resistance of ZC-2 to anthracnose compared with that of ZC3 (Figure 5).

In addition to the transcription factors involved in signal transduction, we analyzed other important differentially expressed transcription factors. Among these genes, WRKY, bHLH, and MYB genes were the predominant transcription factors responsive to C. truncatum infection. However, the response patterns of these three families of transcription factors to C. truncatum infection were not identical. For example, almost all WRKY genes were up-regulated, whereas almost one-third of the bHLH genes and half of the MYB genes were up-regulated (Figure 8A and Supplementary Figure 2).

In addition, certain transcription factors responded differentially to infection by C. truncatum between ZC-2 and ZC3. For instance, many WRKY and MYB transcription factors responded to C. truncatum infection specifically in ZC-2, especially at 4 and 8 HPI, including WRKY6, WRKY15, WRKY33, WRKY40, WRKY42, MYB86, and MYB108. The majority of differentially expressed bHLH genes were identified in both ZC-2 and ZC3. Most were strongly up-regulated in ZC-2 but down-regulated in ZC3, such as bHLH25 and bHLH35, or relatively weakly up-regulated in ZC3, such as bHLH36 and bHLH130, especially at 4 HPI. Several other transcription factor families showed a positive response to C. truncatum infection, such as HSF (heat stress transcription factor), bZIP (homeobox-leucine zipper protein), and GATA genes (Figure 8A). In conclusion, the differential responses of the aforementioned transcription factors to C. truncatum infection between ZC-2 and ZC-3 may be a contributing factor in the differential resistance to anthracnose observed between the two soybean lines.

A large number of putative R and PR genes were identified in response to C. truncatum infection, of which most, especially putative R genes, were strongly up-regulated. Furthermore, most of these putative R and PR genes showed similar response patterns to C. truncatum infection (Figures 8B,C and Supplementary Figure 3). However, certain putative R genes, such as PTI1 (pto-interacting protein 1), RPP13 (disease resistance protein RPS2), RGA2 (disease resistance protein RGA2), RPS6 (disease resistance protein RPS6), ULP2B (disease resistance protein ULP2B), and PR genes, such as CHI (chitinase) and PR14 (lipid transfer protein) genes, responded more promptly and strongly to C. truncatum infection in ZC-2 than in ZC3; this was especially the case for ULP2B genes, which were differentially expressed specifically in ZC-2 (Figures 8B,C). This response may be a crucial factor in the stronger resistance of ZC-2 to anthracnose than that of ZC3.

The KEGG pathway enrichment analysis revealed that pathways involved in the metabolism of terpenoids and polyketides, including “terpenoid backbone biosynthesis,” “carotenoid biosynthesis,” “zeatin biosynthesis,” and “sesquiterpenoid and triterpenoid biosynthesis,” responded to C. truncatum infection more promptly in ZC-2 than in ZC3. Thus, we integrated these pathways and determined the associated DEGs. Most of the DEGs involved in terpenoid backbone biosynthesis and zeatin biosynthesis were up-regulated, whereas the majority of DEGs involved in carotenoid biosynthesis were down-regulated (Figures 9, 10). Most DEGs involved in the metabolism of terpenoids and polyketides showed similar response patterns to C. truncatum infection in ZC-2 and ZC3, whereas genes such as GGPS (geranylgeranyl pyrophosphate synthase), FDPS (farnesyl diphosphate synthase), and TPS14 (terpene synthase 14) (all involved in terpenoid backbone biosynthesis), IPT (adenylate isopentenyltransferase) (involved in zeatin biosynthesis) (Figure 9), and CYP707A (abscisic acid 8′-hydroxylase) and AOG (abscisate beta-glucosyltransferase) (involved in carotenoid biosynthesis) (Figure 10), exhibited a more rapid and/or more intense and/or more prolonged response to C. truncatum infection. Therefore, it was speculated that timeous and intense adjustment of terpenoid backbone and plant hormone metabolism may contribute to improved soybean resistance to anthracnose.

In addition to plant hormones, the relative contents of 109 terpenoid compounds (4 iridoids, 23 monoterpenoids, 23 sesquiterpenoids, 24 diterpenoids, and 35 triterpenoids) were detected by UHPLC-MS (Supplementary Table 5). Forty-three terpenoid compounds (8 monoterpenoids, 10 sesquiterpenoids, 11 diterpenoids, and 14 triterpenoids) were significantly differentially accumulated (VIP > 1, p < 0.05) in at least one comparison group (Figure 11), of which 28 were detected in ZC-2 and 24 in ZC3. The differential accumulation of these 43 terpenoid compounds was mostly plant material-specific and time point-specific. Among these compounds, 19 and 15 compounds were differentially accumulated only in ZC-2 and ZC3, respectively, and 31 terpenoid compounds were differentially accumulated only in one comparison group. Nine terpenoid compounds were significantly differentially accumulated in both ZC-2 and ZC3 (red font in Figure 11). Except for kaurenoic acid and ginsenoside F3, the remaining seven terpenoid compounds responded to C. truncatum infection earlier in ZC-2 than in ZC3. Furthermore, more differentially accumulated terpenoids were identified in ZC-2 than in ZC3 at 4 HPI (Figure 11). These results demonstrated that terpenoid metabolism was highly sensitive to C. truncatum infection and showed timely adjustments. In addition, approximately two-thirds of the differentially accumulated terpenoid compounds in ZC-2 were down-regulated, whereas approximately two-thirds of those differentially accumulated in ZC3 were up-regulated. These results indicated that, in addition to the difference in temporal response of terpenoid metabolism between ZC-2 and ZC3, the trend for change in content also differed.

In addition to certain plant hormones (cis-zeatin, DL-dihydrozeatin, trans-zeatin-riboside, gibberellin A4, gibberellin A7, and abscisic acid), only a small number of terpenoid compounds, including myrcene, picrocrocin, alpha-caryophyllene, abietic acid, kaurenoic acid, terahymanol, taraxerol, and germanicol, were annotated in the pathways involved in terpenoid metabolism (Figures 9, 10). Genes encoding enzymes that directly catalyze the metabolism of these differentially accumulated plant hormone compounds and terpenoid compounds were not differentially expressed. The present results indicate that the transcriptomic data and UHPLC-MS data were not strongly correlated. The metabolism of these plant hormones and terpenoid compounds is closely associated with dimethylallyl-pyrophosphate (dimethylallyl-PP), geranyl-PP, (E,E)-farnesyl-PP, and geranylgeranyl-PP. The genes GGPS, FDPS, and TPS14, which exhibited a more rapid and/or more intense and/or more prolonged response to C. truncatum infection, are directly involved in the metabolism of dimethylallyl-PP, geranyl-PP, (E,E)-farnesyl-PP, and geranylgeranyl-PP. In this regard, the transcriptomic data and the UHPLC-MS data were somewhat similar.

The in vitro antifungal activity test revealed that the terpenoids simvastatin, picrocrocin, and trans-zeatin riboside, and the plant hormone 3-indolebutyric acid significantly inhibited growth of C. truncatum; of these compounds, simvastatin had the strongest inhibitory effect (Figure 12A and Supplementary Figure 4). In addition, exogenous application of the plant hormones 3-indolebutyric acid, abscisic acid, and methyl jasmonate, and the terpenoids simvastatin and picrocrocin significantly increased the resistance of soybean pods to anthracnose; in particular, 3-indolebutyric acid most strongly improved resistance (Figure 12B and Supplementary Figure 5). Thus, 3-indolebutyric acid, simvastatin, and picrocrocin inhibited C. truncatum growth and enhanced the anthracnose resistance of soybean pods.