Two wheat genotypes of Yangmai 18 (Y18) and Yangmai 20 (Y20) were used for this study. The dry seeds of Y18 and Y20 were exposed to ⁷Li-ion beam (Y18LR and Y20LR: 83.6 keV μm^–1) and the mixed high-energy particle field (Y18CR and Y20CR: 1.5 GeV μm^–1) with the same dose of 100 GY. The irradiation treatments were performed by China Institute of Atomic Energy. The irradiated and untreated (wide type, WT) seeds were placed in plastic germination boxes (13 cm × 19 cm × 9 cm) containing four layers of moistened filter paper with sterile water (100 seeds/each box) and grown in a growth chamber with a 14/10 h and 24/20°C day/night light and temperature cycle. The germination rates were recorded after 3 days. After 7 days of growth, the seedling lengths were measured for 20 plants per replicate. Besides, the plant leaves of each treatment were randomly collected and mixed with more than 20 individuals as one biological replicate, frozen in liquid nitrogen immediately and stored at −80°C. Three biological duplicates of each treatment were set.

Total RNA of each sample was isolated using the RNA plant Plus Reagent Kit (TIANGEN, China). The quantity and concentration of each RNA sample were determined by NanoDrop 2000 (Thermo Fisher Scientific, United States). The integrity and purification of RNA samples were qualified using the 2100 Bioanalyzer instrument (Agilent Technologies, United States). Library preparation for RNA-seq was conducted using a MGIEasy mRNA Kit (MGI, China) according to manufacturer’s protocol. Finally, 18 cDNA libraries were sequenced on MGISEQ 2000 platform (GOOALGENE, China) and 150 bp pair-end reads were generated.

Raw reads produced by the sequencer were filtered to remove reads with low quality using Fastp (version 0.19.7) (Chen et al., 2018) and were further assessed for quality of using FASTQC toolkit (v0.11.9) (Brown et al., 2017). The clean reads of each sample were mapped to the wheat reference genome¹ using HISAT2 (v2.1.0) (Kim et al., 2015).

The uniquely aligned reads were used to identify SBSs and InDels variations between the assembled reads and the reference genome sequence using GATK2 software.² To obtain reliable mutations of each mutagenesis progeny, the heterozygous sites of wild type lines were removed in the following analysis. The total numbers of heterozygous and homozygous InDels, SBSs, and different types of SBSs, including transitions (Ti: purine > purine or pyrimidine > pyrimidine) and transversions (Tv: purine > pyrimidine or pyrimidine > purine), were counted for each replicate. Mutation rate was calculated as the average number of mutations per mutant divided by the average length (number of bases) of all genomic regions (Li et al., 2019). The frequency of mutations in every 10-Mb region of each chromosome was calculated and visualized by using Circos software.

The expression level of each gene was estimated by fragments per kilobase of transcript per million fragments mapped (FPKM). The differentially expressed genes (DEGs) between control and irradiated samples were identified using DESeq2 R package (Anders and Huber, 2010). The cutoff of DEGs were defined by using the standard as | log₂^{(fold change)} | ≥ 2 and P-value ≤ 0.05. Gene Ontology (GO) enrichment analysis of the DEGs was performed using singular enrichment analysis tool with FDR < 0.05 by agriGO v2.0, which assigned all DEGs into three principal categories, namely cellular component, molecular function, and biological process (Tian et al., 2017). Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed on the OmicShare platform³ and significant enrichment were selected at P < 0.05. Transcription factors (TFs) were predicted and classified into different families using the PlantTFDB.⁴ Venn Diagrams comparing the number of DEGs across different IR treatments were created using jvenn.⁵ Heatmaps of the gene expressions were illustrated by using TBtools (Chen C. et al., 2020).

The total RNA of all 18 samples was used for the transcriptome analysis and was also used to make cDNas for quantitative real-time PCR (qRT-PCR) validation. The synthesis of first-strand cDNA was conducted using 1 μg of total RNA from each sample with an UEIris RT mix with Dnase (All-in-One) kit (US Everbright, China) according to the manufacturer’s instructions. qRT-PCR was performed using 2× SYBR Green qPCR Master mix (S2014, US Everbright, China) on a QuantStudio™ 7 Flex Real-Time PCR System (Applied Biosystems, United States). Gene-specific primers for qRT-PCR were designed via primer premier 5 and are listed in Supplementary Table 1. qRT-PCR was conducted in triplicate (technical repeats) with three biological replicates for each sample, and the relative gene expression levels were calculated using the 2^–ΔΔCt method.

The total genomic DNA of samples mentioned above were isolated using modified CTAB method (Doyle and Doyle, 1987). The MSAP procedure is performed to investigate global DNA methylation changes according to an established protocol (Tang et al., 2022). Two restriction enzyme combinations, EcoRI/HpaII and EcoRI/MspI (Thermo Fisher Scientific, United States), were used for digestion. After ligated to the HpaII/MspI adapter and EcoRI adapter, two consecutive PCRs were carried out to produce a specific DNA fragment fingerprint. Twenty primer combinations with three selective nucleotides for the EcoRI ends and two to four selective nucleotides for the HpaII/MspI ends were used (Supplementary Table 2). The PCR products were separated by fragment analyze automated CE system (AATI, United States) with the Quick Start Guide 96 Capillary DNF-900 dsDNA Reagent Kit, 35–500 bp (AATI, United States) (Li et al., 2020). The MSAP profiles showing reproducible results between replicates, ranging from 100 to 500 bp, were scored and transformed into a 1/0 binary matrices, where 1 indicates the presence and 0 the absence of a given fragment. Four types of MSAP bands were defined as non-methylation, hemi-methylation, and full-methylation (Table 2) according to Tang et al. (2022).

The germination and seedling height data were statistically analyzed by one-way ANOVA (SNK methods) with a significance level of P < 0.05 using SPSS 18.0 statistical software. The heatmap of gene expression was constructed basing on the Log2FC (RNA-seq) and 2^–ΔΔCt (qRT-PCR) value. The linear correlation coefficient between RNA-seq and qRT-PCR results was detected by using Excel 2010.

Under control condition, the germination rates for Y18 and Y20 were 97.67 and 96.33%, respectively. Both CR and LR irradiation treatments significantly reduced the seed vitality. The germination rates were significantly reduced to an average of 29.67% (23.67% in Y18CR and 35.67% in Y20CR) under CR (P < 0.05), and 51.00% (55.33% of Y18LR and 46.67% of Y20LR) under LR when compared with WT (Figure 1A), suggesting the damaging effect induced by CR is much higher than that induced by LR.

The seedlings of the two wheat genotypes displayed comparable height in the control condition at 7 days after germination (8.57 and 8.1 cm for Y18 and Y20, respectively). The seedling height of both genotypes were significantly reduced under both LR and CR treatments. For both genotypes, seedlings of CR treatment exhibited much more suppressed with an average seedling height reduction of 57.05% compared to the control (1.18 and 3.10 cm for Y18CR and Y20CR, respectively), whilst seedlings with LR treatment only have an average height reduction of 28.45% (6.95 and 5.13 cm for Y18LR and Y20LR, respectively) (Figure 1B). Interestingly, the two wheat genotypes showed clearly different sensitivity to the two types of IR irradiation. In particular, Y18 (18.87% height reduction) displayed to be less sensitive than Y20 (36.63%) to LR treatment. In contrast, Y18 growth was much more severely suppressed (86.19% height reduction) by CR treatment than Y20 (61.73% height reduction) (Figure 1B).

To identify the mutations caused by IR treatments, a total of 18 cDNA libraries (Supplementary Table 3) were constructed for sequencing. After a stringent quality filtering process, an average of 7.60 Gb high-quality clean data (Q30 > 89.49% and uniquely mapped rates > 83.92%) was retained for each sample, which represents ∼57 folds of the total length of the all predicted high-confidence genes in wheat (Supplementary Table 3). Correlation heat map analysis detected high correlations between the biological replicates (Supplementary Figure 1). These results indicated the overall reproducibility and quality of the assay.

To compare the mutagenic effects of LR and CR on wheat, SBSs and InDels between M₁ and WT seedlings were detected. An average of 24,321 mutations were detected for each genotype under LR and CR treatments (Table 1). The number of SBSs are much higher than that of InDels in all samples, accounting 87.46–89.59% of the total mutations. In both wheat genotypes, CR treatment induced higher number of mutations (average 26,002 mutations) and mutation rate (1.79 × 10^–6/bp) than LR (average 22,641 mutations and mutation rate 1.56 × 10^–6/bp) (Table 1). LR induced an average of 11,432 heterozygous SBSs, which was higher than homozygous SBSs (8,396). Whilst the number of homozygous SBSs (average of 12,108) induced by CR is higher than that of heterozygous SBSs (10,985). Similarly, LR induced a higher number of heterozygous InDels (average of 1,680) than that of homozygous InDels (1,135), whereas the number of homozygous (average of 1,438) and heterozygous (1,472) InDels are comparable under CR treatment (Table 1).

Among the SBSs in all target samples, six possible types of substitution were detected. Ti mutations (purine to purine or pyrimidine to pyrimidine, including A/T to G/C and G/C to A/T) (average 65.78% of all SBSs) is much higher than Tv mutations (average 34.22% of all SBSs) (purine > pyrimidine or pyrimidine > purine, including A/T to C/G, A/T to T/A, G/C to T/A, and G/C to C/G) in the all samples. The ratio of Ti/Tv ranges from the lowest (∼1.89) in Y18CR to the highest (∼1.95) in Y20LR (Table 1).

The distribution of the identified mutations (SBSs and InDels) on wheat chromosomes were shown in Figure 2. As a result, CR induced mutations were ranged from 627 (chr4D) to 2,728 (chr6B) and 664 (chr4D) to 1,940 (chr3A) in Y18CR and Y20CR, respectively. On the other hand, LR-induced mutations were ranged from 533 (chr4D) to 2,269 (chr6B) and 626 (chr4D) to 1,705 (chr7A) in Y18LR and Y20LR, respectively. In Y18 and Y20 M₁ seedlings, CR induced an average of 1,226 (Y18CR) and 1,138 mutations (Y20CR) in each chromosome, while LR induced an average of 1,026 (Y18LR) and 1,033 (Y20LR) mutations in each chromosome. Most of the mutations in LR and CR irradiated wheat seedlings were found on the terminal regions of chromosomes in our study (Figure 2). Furthermore, the highest mutation number in Y18 M₁ (Y18LR and Y18CR) seedlings were detected at chr6B [accounting for 10.06% (2,269) and 10.11% (2,728) of total mutations, respectively] and chr6D [accounting for 9.83% (2,219) and 10.01% (2,701) of total mutations, respectively]. Instead, for mutations in Y20 M₁ seedlings, Chr. 7A [accounting for 7.50% (1,705) and 6.70% (1,676) in Y20LR and Y20CR, respectively] and chr3A [accounting for 6.30% (1,432) and 7.75% (1,940) in Y20LR and Y20CR, respectively] seem to be more prone to mutations (Figure 2).

We then investigated the frequencies of mutations (SBSs and InDels) in per 10 Mb regions of wheat genome. The result showed that a highest frequency mutation region, locating on 10–15 Mb of chr3A, was commonly detected in Y18CR and Y20LR, accounting for 0.55 and 0.70% of the total mutations, respectively. Instead, the highest frequency mutation regions of Y18LR and Y20CR were located on 470–473.59 Mb of chr6D and 5–10 Mb of chr3B, accounting for 0.44 and 0.62% of the total mutations, respectively (Figure 2).

To analyze the transcriptional changes related to IRs responses, we compared the transcriptomes of Y20 and Y18 M₁ seedlings exposed to LR and CR to the control without IRs treatment. The results showed that the number of DEGs varied across IR types and wheat genotypes (Figure 3A). In both wheat genotypes, LR resulted in relatively smaller numbers of DEGs (2,250 in Y18LR and 3,491 in Y20LR) than CR (Y18CR 10,618 and Y20CR 8,463). The number of upregulated DEGs were significantly higher than the downregulated DEGs in all samples, with the proportion of up-regulated DEGs ranged from 56.31% (Y18LR) to 83.93% (Y20LR) (Figure 3A).

The cross-comparison between two wheat genotypes Y18 and Y20 showed that most of the IRs induced DEGs were genotype specific. About 41.99% (4,755) DEGs induced by CR were overlapped between Y18 and Y20, while, only 22.49% (1,054) DEGs induced by LR were overlapped between Y18 and Y20 (Figures 3B,D). Of these overlapped DEGs, 4,148 (79.74%) and 447 (8.59%) DEGs were exclusively expressed in response to CR and LR, respectively (Figure 3D), whilst the rest 607 DEGs (11.67%) was responsive to both CR and LR treatments in both genotypes, implying specific and common regulatory changes in wheat responses to different types of IR (Figure 3D).

To further dissect the differential response to IR treatments, the 4,775 CR commonly regulated and 1,054 LR commonly regulated DEGs were used for GO enrichment analysis. Of the 4,755 CR co-expressed DEGs, 65 GO terms were significantly enriched, of the 1,054 LR co-expressed DEGs, 19 GO terms were significantly enriched. Interestingly, most of the enriched GO terms in LR (15) were overlapped with those under CR. These overlapped GOs include nine biological processes (“response to stimulus,” “carbohydrate metabolic process,” “response to abiotic stimulus,” “response to stress,” “response to endogenous stimulus,” “secondary metabolic process,” “response to biotic stimulus,” “response to external stimulus,” and “biosynthetic process”) and nine cellular components (“extracellular region,” “external encapsulating structure,” “cell wall,” “plasma membrane,” “membrane,” and “vacuole”) (Supplementary Table 4).

The KEGG pathway enrichment analysis showed that, a total of 19 and 25 pathways were significantly enriched under CR and LR irradiation, respectively (Figure 4). Among these pathways, 10 common pathways were observed in both IR responses, which include “biosynthesis of secondary metabolites,” “metabolic pathways,” “fatty acid elongation,” “phenylpropanoid biosynthesis,” “cutin, suberine, and wax biosynthesis,” “biosynthesis of unsaturated fatty acids,” and “starch and sucrose metabolism.”

In addition to the above common pathways, nine pathways including “glyoxylate and dicarboxylate metabolism,” “carbon fixation in photosynthetic organisms,” “DNA replication,” “carbon metabolism,” “pentose and glucuronate interconversions,” and “photosynthesis – antenna proteins,” were uniquely enriched under CR treatment (Figure 4A), and 15 pathways were uniquely enriched under LR treatment, including “phenylalanine metabolism,” “amino sugar and nucleotide sugar metabolism,” “anthocyanin biosynthesis,” and “flavone and flavonol biosynthesis” (Figure 4B). Taken together, these results highlighted the involvement of both common and unique pathways in LR and CR treatments responses.

Among the detected DEGs, a total of 101 genes related to phenylpropanoid biosynthesis were found to be differentially expressed under CR and LR treatments (Figure 5A and Supplementary Table 5). The 101 phenylpropanoid related genes included 73 peroxidases, 10 beta-glucosidases, 9 phenylalanine ammonia-lyases, 3 cinnamoyl-CoA reductases, 2 4-coumarate-CoA ligases, 2 scopoletin glucosyltransferases, 1 shikimate O-hydroxycinnamoyl transferase, and 1 cinnamyl alcohol dehydrogenase (Figure 5A and Supplementary Table 5). Of the phenylpropanoid related genes, 21 genes were co-regulated by LR and CR, including 5 phenylalanine ammonia-lyases and 16 peroxidases (Supplementary Table 5). Among the phenylpropanoid related genes, 3 phenylalanine ammonia-lyase and 2 peroxidase genes were genotype specific and mainly upregulated in Y20 under CR and LR treatments, while downregulated in Y18 under CR and LR treatments. Moreover, 34 phenylpropanoid related genes, including 1 phenylalanine ammonia-lyase, 1 cinnamoyl-CoA reductase, 1 cinnamyl alcohol dehydrogenase, 5 beta-glucosidases, 2 scopoletin glucosyltransferases, and 24 peroxidases, were exclusively expressed under CR irradiation. Noteworthy, most of the 73 peroxidase genes were found significantly upregulated at one or both IR treatments (Figure 5B and Supplementary Table 5). In addition, 4 DEGs encoding 1 4-coumarate–CoA ligase, 1 beta-glucosidase, and 2 peroxidases were only differentially expressed under LR irradiation (Supplementary Table 5).

Moreover, 24 genes encoding enzymes involved in ROS metabolism were also differentially expressed, mainly including glutathione S-transferase (GST), ascorbate peroxidases (APX), and glutathione peroxidase (GPX) (Figure 5C). Twenty-one GSTs were differentially expressed, and over half of them were up-regulated under LR or CR treatments. In which, 4 GSTs were co-upregulated at both IR treatments in two genotypes, while 3 and 3 GSTs were specifically up- and down-regulated by CR treatment, respectively. Additionally, 2 APXs and 1 GPX were only induced by CR treatment.

In addition to the phenylpropanoid biosynthesis pathway, KEGG enrichment results also showed that genes related to starch and sucrose metabolism were enriched significantly in both LR and CR treatments (Figure 4). A total of 31 DEGs were identified involved in four significant metabolic processes in “starch and sucrose metabolism,” which include “cellulose degradation” (22), “trehalose degradation” (2), “sucrose degradation” (3), and “starch degradation” (4) (Figure 6 and Supplementary Table 6). In the “sucrose degradation” process, 3 DEGs (2 sucrose synthases and 1 beta-fructofuranosidase) were significantly upregulated in the both IR treatments, while the DEGs for “cellulose degradation” and “trehalose degradation” processes were mainly upregulated under CR treatment. Notably, one beta-amulase gene involved in “starch degradation” process was co-upregulated in both IR treatments, while the expression abundance was higher in Y20 than that in Y18 under both CR and LR treatments (Figure 6 and Supplementary Table 6).

Transcription factors play crucial functions in stress adaptive signaling cascades and control the expression of numerous genes. A total of 285 TF-encoding DEGs (35 TF families) were found for CR irradiation, whilst only 63 TF-encoding DEGs (13 TF families) were found for LR irradiation (Figure 7A). The most abundant TF family is MYB (42 and 21 for CR and LR, respectively) family, followed by bHLH (38 and 7), ERF (20 and 5), and bZIP (16 and 2) (Figure 7A), sequentially. A total of 32 TF-encoding DEGs (11 TF families) were commonly detected under both IR treatments, including MYB (13), bHLH (4), GRAS (3), bZIP (2), MYB_related (2), HD-ZIP (1), GATA (1), B3 (1), LBD (2), ERF (2), and WRKY (1) (Figure 7B). Among these common TFs, majority were co-upregulated in both genotypes under LR and CR treatments (Figure 7B), with the exception of 2 MYB_related, 2 LBD, and 2 ERF TFs, which were down-regulated in CR and LR (Figure 7B). In addition, 190 TFs belonging to 35 families were induced only by CR treatment, such as TCP (10), AP2 (7), G2-like (7), Trihelix (5), HSF (3), and ARF (3). While 4 TFs, including 1 bHLH, 1 WRKY, and 2 MYB genes were specifically induced by LR (Supplementary Table 7).

To verify the accuracy of our RNA-seq, 15 DEGs in phenylpropanoid biosynthesis, starch and sucrose metabolism, and TFs were randomly selected (5 of each) for qRT-PCR validation. The results showed that the expression patterns were generally consistent with the RNA-seq data, suggesting that the high accuracy and reproducibility of our RNA-seq data (R² = 0.9011) (Supplementary Figure 2).

Cytosine methylation patterns in the leaves of two wheat genotypes under IR treatments were detected using 20 pairs of primers. An average of 383 clear and reproducible bands were detected (Table 2 and Supplementary Figure 3). Most of the CCGG sites were shown to be largely methylated with the values ranging between 74.15 and 83.29%. In the control, we observed slight variations in the numbers of methylated sites between the Y20 (78.07%) and Y18 (79.11%) genotypes. Compared with control, LR displayed increased DNA methylation levels in both genotypes (83.29 and 80.89% in Y18LR and Y20LR, respectively), while CR-treated samples decreased in the methylation levels (75.46 and 74.15% in Y20CR and Y18CR, respectively). Further analyses showed that fully methylated bands were more predominant than the hemi-methylated ones.

To further investigate the difference of wheat DNA methylation in response to IR treatments, 16 possible banding patterns between control and IR treatments were identified and classified into three groups: no change, hypomethylation, and hypermethylation. Results (Table 3) showed that ∼47.96% of the CCGG sites remained unchanged under IR treatments. The average percentage of hypomethylated bands was 28.85% under CR treatment, higher than that of LR treatment (23.45%). In contrast, the average percentage of hypermethylated bands was 21.36% under CR treatment, lower than that of LR treatment (30.43%) (Table 3). Additionally, different levels of DNA methylation patterns between Y18 and Y20 were also observed under both IR treatments. For example, the hypomethylation level in Y18CR was higher than that for Y20CR, while more hypomethylation events under LR occurred in Y20LR than Y18LR.

Ionizing radiation has been considered the most powerful source of mutagenesis for improving agricultural traits in various crops worldwide. The mutation spectrum and frequency are important factors during the selection of appropriate mutagens for mutation breeding and gene functional mutation (Kazama et al., 2017; Ichida et al., 2019; Jo and Kim, 2019; Li et al., 2019). The mutation frequencies in the genome of M₂ rice plants generated by gamma rays and carbon-ion beams irradiations were estimated to be 3.2 × 10^–8 and 2.4 × 10^–8/bp, respectively (Li et al., 2019). While the carbon-ion beams induced mutation frequencies reported by Yang et al. (2019) and Oono et al. (2020) were average 2.4 × 10^–7/bp (M₂) and 2.7 × 10^–7/bp (M₂ and M₃), respectively. More recently, proton beams and gamma rays induced mutation frequencies of M₂ mutants in rice were estimated ∼5.0 × 10^–7 and ∼7.0 × 10^–7/bp, respectively (Lee et al., 2021). In Arabidopsis, mutation rates induced by gamma rays (M₂–M₆) and carbon-ion irradiation (M₂) were 1.4–2.2 × 10^–7 and 0.9–1.2 × 10^–7/bp (Hase et al., 2020). However, the mutation rates obtained in the present study were higher than those estimated by these previous studies (Li et al., 2019; Yang et al., 2019; Hase et al., 2020; Oono et al., 2020; Lee et al., 2021). The present study demonstrates that LR and CR radiation are all effective for mutation induction in wheat, suggesting that CR and LR could be good mutagens for plant mutation breeding in the future. The higher mutation frequencies in this study could be partially explained by the higher LET applied in the present study, and partially by other factors such as the different species or mutant generation used in different studies.

It should be noted that RNA-seq is limited to the transcribed genetic regions, while those mutations located in the intron and non-coding regions and large size mutations are not covered. Therefore, it is understood that the calculated mutation rates in this study may be slightly lower than those based on whole genome sequencing. Whole genome sequencing may be necessary to compare the mutagenesis profiles induced by CR and LR irradiation in crops in the future.

The distribution of IR-induced mutations on chromosomes have been discussed extensively (Li et al., 2016; Tan et al., 2019; Xiong et al., 2019; Yang et al., 2019). Fast-neutron (Li et al., 2016), gamma rays, and carbon-ion beams induced mutations in rice (Yang et al., 2019) were evenly distributed across the genome. In contrast, CR- and LR-induced mutations are biased toward one or both ends of most chromosomes, which is consistent with the previous studies in wheat that LR-induced mutations in Jing411 and Heyou1 were found on the terminal regions of chromosomes (Xiong et al., 2019). The detection method may be the potential reason of this inconsistent result. Whole genome sequencing and resequencing were used in the studies of Li et al. (2016) and Yang et al. (2019), respectively. However, RNA-seq was used in our and Xiong’s studies, which is limited to the transcribed genetic regions and associated with the distribution of genes (Figure 2). Monroe et al. (2022) revealed the natural genetic variation in mutation bias, which associates with GC content, methylated cytosines and gene structure (Monroe et al., 2022). Similarly, Weng et al. (2019) found that spontaneous mutation profile in Arabidopsis do not occur evenly across the genome: they are biased toward G: C to A: T transitions (Weng et al., 2019). Similar results were observed in the gamma rays, carbon-ion beams, LR, and proton beams induced mutations (Guo et al., 2019; Li et al., 2019; Xiong et al., 2019; Yang et al., 2019; Oono et al., 2020; Choi et al., 2021; Lee et al., 2021), which is also consistent with the present study. Therefore, the uneven distribution of IR induced mutations may be related to mutational bias. On the other hand, the differences between different materials suggest that genotypes may be involved in this uneven distribution, but more evidence and experiments are needed to confirm this.

When plants are irradiated by IR, the intracellular water is decomposed into ROS, such as superoxide anion radicals, hydroxyl radicals, hydrogen peroxide (Wang et al., 2018). A significant part (about 70–80%) of the IR-related DNA damage is caused by ROS formed during radiolysis of water and only 20–30% of the damage is due to the direct absorption of high-energy IR quanta by the target DNA molecules (Gudkov et al., 2019). To alleviate oxidative damage, plant significantly activated the antioxidant defense system using ROS detoxification of antioxidants (e.g., peroxidases and catalases) and osmotic adjustment substance (e.g., soluble sugar and proline) to maintain cellular ROS (Luo et al., 2019). For example, the exposure of barley seeds to gamma rays results in the accumulation of hydrogen peroxide in seedlings and activating the antioxidant system (Volkova et al., 2019). An RNA-seq-based study of genes differentially expressed during LR-irradiation revealed that LR-induced DEGs are associated with “phenylpropanoid biosynthesis” and “antioxidant process” (Xiong et al., 2019). Variations in phenylpropanoid and flavonoid biosynthesis pathways in response to UV-B radiation were found in date palm leaves (Maher et al., 2021).

Phenylalanine ammonia-lyase catalyzes the first step in the phenylpropanoid pathway, which plays an important role in the production of antioxidant phenolic compounds such as flavonoids and tannins. Acute gamma irradiation promoted the accumulation of the H₂O₂ and malondialdehyde content of rice plants, and increased the enzyme activities of phenylalanine ammonia-lyase, SOD, POD, CAT, and APX (Choi et al., 2021). Consistently, in this study, 5 phenylalanine ammonia-lyases and 16 PODs involving phenylpropanoid biosynthesis pathway were co-regulated by LR and CR.

Glucose is one of the soluble sugars in plants, which plays an important role in the plant osmotic-regulation under stresses (Khaleghi et al., 2019). The increase of soluble sugar content and the enrichment of “starch and sucrose metabolism” pathway were also reported in rice under gamma rays, cosmic rays, and carbon-ion beams treatments (Hwang et al., 2014), which is consistent with our result. The contents of soluble sugar, such as sucrose and D-glucose-6-phosphate, were increased in Porphyra haitanensis under UV-B exposure (Fu et al., 2021). Combined with these studies, we suggested that the accumulation of soluble sugar and the activation of ROS scavenging processes might be the important processes for plants in response to IR irradiation.

Under LR and CR exposure, we observed clear differences in the numbers and types of DEGs in two wheat genotypes irradiated with CR and LR (Figure 3). The number of DEG increased with the LET level, and showed a higher number of DEGs in CR treated samples than that of LR (Figure 3), implicated the type-specific responses in plants to IR irradiation exists.

In cowpea, the proton-beam treatment induced more DEGs than that of gamma-rays, and more diverse in terms of pathways were observed in the proton-beam treatment than gamma treatment. The “oxidation-reduction process” and “proteolysis” were the most enriched terms in GR treatment, and “substance metabolic process” and “cellular metabolic process” were the most enriched terms in proton-beam treatments (Kang et al., 2021). In duckweed plants, genes related to “anthocyanin accumulation” and “ATP synthesis” were specifically regulated in response to uranium treatment. While genes involved in “DNA damage and repair” and “calcium signaling” were specifically regulated in response to gamma radiation (Fu et al., 2019). In addition, the divergent response pathways in rice plants response to gamma rays, cosmic-ray, and carbon-ion beams treatments were also reported (Hwang et al., 2014), in which, genes related to “lipid metabolic process” and “gibberellin metabolic process” are greatly enriched in gamma rays treatment, while the stimulus responsive genes were significantly enriched in the carbon-ion beams irradiated plants. The molecular function pathways, such as “transcription regulator activity” and “oxidoreductase activity” were mainly enriched in the cosmic rays treatment (Hwang et al., 2014). These results suggest that plants have specific response pathways to different types of IR and that these response pathways are species-dependent (Hwang et al., 2014; Fu et al., 2019; Xiong et al., 2019; Kang et al., 2021).

Transcription factors are crucial components in signal transduction, and directly control the expression of specific sets of downstream stress-responsive genes (Luo et al., 2019). They are triggered by various signal transduction pathways and can bind to cis-acting elements directly or indirectly to modulate the transcription efficiency of target genes (Maher et al., 2021). In the previous study, TFs, such as MYB, WRKY, NAC, and bHLHs, were identified as regulatory proteins that are involved in regulating the expression of other genes that participate in the UV-B stress response (Maher et al., 2021). Consistently, these TFs were also detected in this study, and several IR type specific TFs were also identified between CR and LR, such as ARR-B, GRF, and HSF. Similar results were also detected in cowpeas, in which, proton-beam treatment induced more types of TFs than that of gamma-ray, and TFs, including ARR-B, B3, bZIP, C2H2, CO-like, DBB, G2-like, GRF, HSF, MYB, MYB_related, NAC, Trihelix, WOX, and WRKY, were regulated only by the proton beam when compared with gamma rays treatment (Kang et al., 2021). In rice, ARR-B and PHOR1 TFs were specifically induced by carbon-ion beams and gamma rays, respectively (Hwang et al., 2014). The results indicated that plants complexly regulated by varying the combination and concentration of TFs according to the IR and the crop species.

Plant response to DNA damage through epigenetic modifications has been well documented (Kim, 2019). Epigenetic mechanisms, such as histone modification and DNA methylation, involving response to radiation, have also been demonstrated (Shi et al., 2014; Georgieva et al., 2017; Pan et al., 2017; Zhao et al., 2018; Horemans et al., 2019; Marfil et al., 2019). Gamma rays have been shown to cause local or global changes in the chromatin structure, including an active (H4K12ac, H3K36me3, H3K4me3, H3K4ac, and H3K27ac) or repressive (H3K27me3) chromatin state (Pan et al., 2017). In addition, chernobyl radio-contaminated Pine trees (Volkova et al., 2018) and soybean seedlings (Georgieva et al., 2017) showed an increase in global DNA methylation, while a significant decrease in Arabidopsis (Horemans et al., 2018). UV-B induced hypermethylation in the grapevine (Marfil et al., 2019), while hypomethylation in sweet wormwood (Pandey and Pandey-Rai, 2015). In this study, DNA methylation levels increased in both wheat genotypes under LR treatment, while decreased under CR treatment. These contrasting reports imply that different methylation mechanisms for radiation response may exist in plant species.

Further analysis of DNA methylation patterns showed that DNA hypomethylation and hypermethylation are concurrent in the both LR and CR treatments, with LR inducing a higher rate of DNA hypermethylation than hypomethylation. Notably, this pattern was reversed under CR irradiation. Similar results were reported in rice under carbon-ion radiation (Zhao et al., 2018), which showed that low-dose carbon-ion radiation (0.01, 0.2, or 1 Gy) induced higher proportion of hypermethylation than hypomethylation, whereas high-dose carbon-ion radiation (2, 5, or 20 Gy) induced more hypomethylation than hypermethylation. These results suggest that the epigenetic regulation patterns were quite complex and may vary across species depending on the different types of IRs.