A Comparison of the Transcriptomes of Cowpeas in Response to Two Different Ionizing Radiations

In this study, gene expression changes in cowpea plants irradiated by two different types of radiation: proton-beams and gamma-rays were investigated. Seeds of the Okdang cultivar were exposed to 100, 200, and 300 Gy of gamma-rays and proton-beams. In transcriptome analysis, the 32, 75, and 69 differentially expressed genes (DEGs) at each dose of gamma-ray irradiation compared with that of the control were identified. A total of eight genes were commonly up-regulated for all gamma-ray doses. However, there were no down-regulated genes. In contrast, 168, 434, and 387 DEGs were identified for each dose of proton-beam irradiation compared with that of the control. A total of 61 DEGs were commonly up-regulated for all proton-beam doses. As a result of GO and KEGG analysis, the ranks of functional categories according to the number of DEGs were not the same in both treatments and were more diverse in terms of pathways in the proton-beam treatments than gamma-ray treatments. The number of genes related to defense, photosynthesis, reactive oxygen species (ROS), plant hormones, and transcription factors (TF) that were up-/down-regulated was higher in the proton beam treatment than that in gamma ray treatment. Proton-beam treatment had a distinct mutation spectrum and gene expression pattern compared to that of gamma-ray treatment. These results provide important information on the mechanism for gene regulation in response to two ionizing radiations in cowpeas.


Introduction
Ionizing radiation has been considered the most powerful source of mutagenesis for improving agricultural traits of various crops worldwide [1]. Since the 1960s, gamma-rays and X-rays have been commonly used to induce mutations during plant breeding [2]. In recent years, new mutant-derived cultivars have been developed using novel mutagens, such as cosmic rays and ion beams. These ionizing radiations have created various mutations that achieve high yield [3], early maturity [4], improvement of crop quality and nutritional traits [5], and resistance to biotic [6] and abiotic stress [7].
Gamma-rays are electromagnetic radiation with 0.2 keV µm -1 linear energy transfer (LET; the energy transferred per unit length), which can penetrate tissues [8]. The ion beams are composed of particles of various masses from protons to uranium atoms generated through particle accelerators and can have high-LET ranging from 22.5 keV µm -1 to 4000 keV µm -1 [8,9]. Compared with gamma-rays, ion beams show high relative biological

Transcriptional Variations Induced by Two Different Ionizing Radiations in Cowpeas
To compare transcriptional variations induced by gamma-rays and proton-beams in the cowpeas, RNA-sequencing analysis was performed with a total of 21 samples of control and irradiation treatments with three different doses (100, 200, and 300 Gy) for the two radiations. Each treatment was repeated independently three times. As a result of RNA-seq (Table S1), 600,059,830 transcriptome short reads (average length 101 bp) were collected from all samples. In the transcriptome short reads, bases with a phred score (Q) of less than 20, representing base quality, were trimmed, and reads with a trimmed read length of less than 25 bp were removed. The final value of the total length of trimmed reads/total length of raw reads was 80.00% using this preprocessing method. As a result of calculating the expression value by mapping 537,045,136 cleaned reads to 42,287 reference transcripts of cowpeas (Vunguiculata_469_v1.1), the mapping rate was approximately 91.46%. Additionally, among the 42,287 standard genes used for the analysis, 40,873 genes were expressed, and among them, 38,385 (93.91%) were genes with functional descriptions. These data indicated that reliable transcriptome data were available for subsequent differential analysis.
DEGs were screened between the treatments using log2FC (fold change) > 1 and p-adj (adjusted p-value) < 0.05 (Tables S2-S12). The value of log2FC greater than 1 was defined as upregulation, and less than 1 was defined as downregulation. The number of DEGs selected between the comparisons is shown in Figure 1. The 32, 75, and 69 DEGs at each dosage (100 Gy, 200 Gy, and 300 Gy) of gamma-ray irradiation compared with the control were identified. In contrast, 168, 434, and 387 DEGs were selected for each dose of proton-beam irradiation compared with the control. In the case of both sources, the number of DEGs was highest in the 200 Gy treatment, and slightly decreased in the 300 Gy treatment. Compared to that of gamma-rays, at all doses, the number of up-regulated and down-regulated DEGs were significantly higher in proton-beam treatments. In the case of gamma-ray treatments (Figure 2a,b, Tables S8 and S9), 3, 26, and 28 DEGs were dose-specifically up-regulated in a dose-dependent manner for the 100 Gy, 200 Gy, and 300 Gy doses, and 1, 14, and 26 DEGs were down-regulated in a dose-dependent manner. Only eight genes were up-regulated in common at all doses, and these genes mainly encoded the Subtilisin-like serine endopeptidase family protein. In proton-beam treatments (Figure 2c,d, Tables S10 and S11), 26, 245, and 196 DEGs were up-regulated in a dose-dependent manner, and 6, 65, and 53 DEGs were down-regulated in a dosedependent manner with 100 Gy, 200 Gy, and 300 Gy doses, respectively. A total of 61 DEGs were commonly up-regulated at all doses, and these genes encoded proteins, heat shock transcription factor A6B, the chaperone DnaJ-domain superfamily, NAC transcription factor-like 9, and beta glucosidase 15. One commonly down-regulated DEG encoded the nodulin MtN21/EamA-like transporter family protein. Overall, 133 and 759 DEGs were detected in gamma-ray treatments and proton-beam treatments, respectively. Among these, 95 overlapping DEGs were commonly detected in at least one or more doses for both radiations (Table S12). These DEGs included the GDSL-like Lipase/Acylhydrolase superfamily protein, NAD(P)-binding Rossmann-fold superfamily protein, lipoxygenase 1, and Subtilisin-like serine endopeptidase family protein.
There were no common DEGs in all doses of both sources. However, different DEGs encoding the same protein have been identified, and the definition of this protein was a subtilisin-like serine endopeptidase family protein. Interestingly, the DEGs encoding the heat shock protein 90.1, the heat shock transcription factor A6B, and the Chaperone DnaJ-domain superfamily protein, which increased in a dose-dependent manner in the proton-beam treatment, were not differentially expressed at all doses, but were differentially expressed only at 300 Gy in the gamma-ray treatments.
The qRT-PCR was performed to validate RNA-seq. A total of 15 genes were selected from the DEGs that commonly detected in gamma-ray treatments and proton-beam treatments ( Figure S1). The relative expression levels of 15 genes obtained by qRT-PCR were similar to results of RNA-sequencing analysis.

Functional Categorization of DEGs Induced by Two Ionizing Radiations in Cowpeas
To further understand the function of DEGs, Gene Ontology (GO) term enrichment analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment were performed (Tables S2-S7). The GO term is presented in three independent categories Biological Process (BP), Molecular Function (MF), and Cellular Component (CC). Figures 3 and 4 showed the top 20 enriched GO terms of the up/down-regulated genes by gamma-rays and proton-beams, respectively.
In a result of each treatment excluding overlapping up-regulated and down-regulated terms, 100 Gy, 200 Gy, and 300 Gy of the gamma-ray treatments showed 21, 81, and 60 terms, and the proton beam treatment showed 97, 157, and 79 terms. Overall, protonbeam treatments exhibited more diverse GO terms compared to gamma-ray treatments.
Comparing the GO results of up-regulated genes in the 200 Gy gamma-ray and protonbeam treatments, the most DEGs were detected in metabolic process and biological process terms in BP from both sources. However, differences existed in the next level. For gammaray treatments, oxidation-reduction process and proteolysis were the most enriched terms and for proton-beam treatments, organic substance metabolic process and cellular metabolic process were the most enriched terms. In CC, the apoplast the extracellular region were the most enriched terms in gamma-ray treatments, and membrane part and membrane protein complex were the most enriched terms in proton-beam treatments. In MF, catalytic activity, hydrolase activity, cation binding, metal ion binding, and oxidoreductase activity were similarly over-enriched in both sources, except for the molecular function term in gamma-ray treatments. Tables 1 and 2 show the enriched KEGG pathways of the up/down-regulated genes by gamma-ray and proton-beam. The DEGs were most involved in biosynthesis of other secondary metabolite pathways in all treatments. In proton-beam treatments, the next major pathways involved were environmental adaptation and signal transduction in 100 Gy, carbohydrate metabolism and energy metabolism in 200 Gy, environmental adaptation, folding, and sorting and degradation in 300 Gy. However, in the case of gamma-ray treatments, metabolism of terpenoids and polyketides in 100 Gy, lipid metabolism in 200 Gy, and environmental adaptation in 300 Gy were the prominent pathways. In particular, genes involved in environmental adaptation (plant-pathogen interaction, circadian rhythmplant) pathway were generally up-regulated in all treatments of the proton beam treatment. However, in the gamma-ray treatment group, the number of up/down-regulated genes was the same or there were more down-regulated genes.   Table 1. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the up-and down-regulated genes of cowpeas exposed to different dose rates of gamma-rays.  Table 2. Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis of the up-and down-regulated genes in cowpeas exposed to different dose rates of proton-beams.

Clustering Analysis of DEGs Induced by Two Different Ionizing Radiations in Cowpeas
A clustering analysis was performed to confirm the expression pattern of genes using differential expression gene information in 100 Gy vs. the control, 200 Gy vs. the control, and 300 Gy vs. the control comparisons for the two sources. The DEGs of the two sources were classified into six clusters each. In gamma-ray treatments (Figure 5a), the 82 of 133 DEGs belong to cluster 2 (C2). They showed a higher expression level than that of the control at all doses, and the level increased at 200 Gy and decreased at 300 Gy. In the GO analysis of these genes, most were distributed in biological process in BP, apoplast and extracellular region in CC, catalytic activity in MF. Additionally, KEGG analysis revealed that most of the genes were involved in biosynthesis of other secondary metabolites, lipid metabolism, and signal transduction. In proton-beam treatments (Figure 5b), 328 of 759 DEGs belonged to cluster 1(C1) and had a similar pattern to that of C2 of gammaray treatments. The GO analysis showed that most of the genes were distributed in the biological process in BP, membrane in CC, and catalytic activity in MF. Additionally, KEGG analysis revealed that most of the genes were involved in biosynthesis of other secondary metabolites, similar to that of gamma rays, followed by energy metabolism and carbohydrate metabolism. However, the proportion of genes involved in signal transduction and lipid metabolism was relatively small.
In contrast to the abovementioned clusters, C4 (16 DEGs) of the gamma-ray treatments showed lower expression levels than did the control at all doses, and these DEGs distributed in protein phosphorylation and phosphorylation in BP, and ion binding and catalytic activity in MF. Additionally, they were mainly involved in environmental adaptation, biosynthesis of other secondary metabolites and lipid metabolism pathway. C2 of protonbeam treatments showed a similar trend, including 40 DEGs. It also showed distribution in peptidase regulator activity, endopeptidase regulator activity, peptidase inhibitor activity, and endopeptidase inhibitor activity in MF. Furthermore, they were mainly involved in biosynthesis of other secondary metabolites, and metabolism of terpenoids and polyketides. The environmental adaptation pathway, which was most often observed in the C4 of gamma-ray treatments, was not seen in C2 of proton-beam treatments, but most in C3.
Additionally, the C5 (118 DEGs) of proton-beam treatments and C3 (six DEGs) of gamma-ray treatments had a pattern of increasing expression levels with dose, and the term oxidoreductase activity in MF was most commonly observed in both. However, the C4 (94 DEGs) of proton-beam treatments and the C1 (26 DEGs) of gamma-ray treatments, which showed a pattern of decreasing expression levels with dose, represented the term major hydrolase activity and protein dimerization activity in MF, respectively. Detailed information on the various clusters of each source is shown in Tables S13 and S14.

Target Gene Analysis of DEGs Induced by Two Different Ionizing Radiations in Cowpeas
To determine which signaling pathways were activated in response to the two ionizing radiations, we investigated the significant differences in the expression of genes related to ROS, plant hormones, defense signals, photosynthesis, and plant transcription factors (Tables 3 and 4). For both sources, the largest number of DEGs were detected in plant hormones and transcription factors. For plant hormones, the genes related to SA and ABA were the most differentially expressed under both radiations. Interestingly, the AUX, CK, and BR genes were regulated only in proton-beam treatments, but not in gamma-ray treatments. Among them, CK-related DEGs were up-regulated at all doses, and none were down-regulated. For the TFs, 23 were regulated for both ionizing radiations. TALE was controlled only by gamma rays. bHLH, C3H, Dof, ERF, LBD, MIKC_MADS, and M-type_MADS were regulated by both sources. 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. Table 3. The differently expressed genes (DEGs) list as related to defense, photosynthesis, and ROS of cowpeas exposed to different dose rates of gamma-rays and proton-beams.

Category
Gene     The photosynthesis target genes were obtained only in 200 Gy of the proton-beam treatments. Most of these genes were matched to carbon fixation in photosynthetic organisms. However, both the ROS and defense genes were significantly up-regulated by 300 Gy in the proton beam treatments compared to that of gamma rays. The ROS-related genes included glutathione S-transferase TAU 8, 19, and peroxidase superfamily proteins. Furthermore, the defense-related genes included cation exchanger 1, cation exchanger 3, and β,-1,3-glucanase 1.

Discussion
Gamma rays have been commonly used as a mutagen to improve agricultural traits of crops [2]. Recently, it has been reported that proton-beam irradiation induced a different mutation spectrum compared to that of gamma-rays [14]. This means that the development of mutation breeding technology using proton-beams can lead to expansion of the genetic resource pool. Gamma-rays and ion beams caused different gene expression in plants that lead to changes in type, structure, and activity of proteins [23,24]. Therefore, transcriptome analysis can be a very useful approach to understanding physiological characteristics in response to stimuli. Several mutation breeding studies have been attempted to further understand the response to ionizing radiation through transcriptome analysis [24]. However, existing gene expression information is insufficient to fully understand the molecular mechanisms of gamma-rays and proton-beams. Therefore, this study was conducted to investigate the transcriptional variations and the functions of the DEGs induced by proton-beams and gamma-rays in cowpeas.
As a result of analyzing transcriptional variation in cowpeas irradiated with two ionizing radiations, there were clear differences in the numbers and types of DEGs. The numbers of DEGs were significantly higher with proton-beam treatments compared to that of gamma-rays. In both treatments, the largest number of DEGs were detected at 200 Gy, and the number of DEGs decreased at the next higher dose (300 Gy). The number of up-regulated DEGs was higher at all doses and sources than that of down-regulated DEGs. However, the opposite tendency was observed in rice irradiated with gamma rays and ion beams [25]. These DEGs were involved in a wide variety of biological functions such as plant development, abiotic stress responses, signaling, and plant secondary metabolism. Among them, genes encoding heat shock protein 90.1, the heat shock transcription factor A6B, and the Chaperone DnaJ-domain superfamily protein were dose-dependently expressed at all doses in the proton-beam treatments and differentially expressed only at 300 Gy of gamma-ray treatments. In rice, cosmic-ray and ion-beam treatments induced up-regulation of heat shock protein (HSP) genes, but gamma-ray treatments induced down-regulation of HSP genes [25]. The HSP and heat shock transcription factor were known to protect macromolecules to counteract stress, break down and recycle components of irreversibly damaged cells, and promote proper adaptive responses [26,27]. It was suggested that the misfolded proteins were generated or accumulated more in the protonbeam treatments than the gamma-ray treatments, and the degree increased in proportion to the dose. There were no common DEGs in all treatments of both sources. However, the different DEGs encoding the same protein have been identified, and the definition of this protein was a subtilisin-like serine endopeptidase family protein. This enzyme in Arabidopsis thaliana has been reported to be involved in the regulation of stomatal density and distribution [28] and auxin-induced lateral root formation [29]. Additionally, this gene was associated with stress-induced senescence in wheat [30].
The LET of an ion beam, such as a proton beam, is higher than that of a gamma ray, and this ion beam results in a higher mutation frequency and different mutation spectrum compared to that of gamma rays [11,31]. Furthermore, radiation with higher LET is known to cause more serious and more complex damage to DNA [32]. In the GO and KEGG analysis, the two treatments showed a different spectrum of terms and pathways, and a different number of related genes. In general, more terms and various pathways were revealed in the proton-beam treatments than the gamma-ray treatments. In other words, gene expression in cowpeas was more diversely and complexly regulated by a proton-beam with a higher LET than by gamma-rays, and the main mechanisms of each source were different. Additionally, it was complexly regulated at high doses (200 and 300 Gy) compared to low doses (100 Gy) in each source. Specifically, the DEGs involved in carbohydrate metabolism and energy metabolism pathways were excessively accumulated at 200 Gy and decreased at 300 Gy in the proton-beam treatments. However, the DEGs involved in environmental adaptation and folding, sorting, and degradation pathways were most prominent in all treatments (especially 300 Gy) of proton-beam compared to gamma-ray treatments. These results were similar to those of Lemna minor exposed to gamma rays, which changed the gene expression mechanism from acclimatization to a survival response by controlling energy distribution as the dose increased [33]. In the gamma-ray treatments, no corresponding trend was observed, which was considered a difference in crops.
As a result of confirming the expression pattern of DEGs through cluster analysis, most DEGs belonged to cluster 2 in the gamma-ray treatments. In proton-beam treatments, most of the DEGs belonged to cluster 1, which showed a similar pattern as did cluster 2 of gamma-rays. In both sources, the expression level of genes increased at 200 Gy and decreased at 300 Gy. This suggested that both sources may have a turning point in which gene expression tends to change with the dose of ionizing radiation. However, the function of the gene was not the same. The lipid metabolism (linoleic acid metabolism) pathway mainly observed in C2 of gamma-rays was observed in C5 of proton-beams. Additionally, many genes involved in energy metabolism (carbon fixation in photosynthetic organisms, photosynthesis) were observed in C1 of the proton beams but were not prominent in any clusters of gamma rays. This means that the response mechanism varies by controlling the expression of major genes according to the source and dose.
According to target gene analysis, the cowpeas stimulated by two ionizing radiations were most affected by plant hormones and TFs related genes. Hormones are involved in a wide range of biological processes and act as a defense in response to stress. ABA exhibits a prominent defense response against abiotic stresses, such as salinity and drought [34,35], and SA, JA, and ET are effective against biotic stresses, such as various pathogens and pests [36]. Additionally, AUX, GA, and CK are important hormones that promote plant growth [37]. In this study, ABA, SA, JA, ETH, and GA were regulated in gamma-ray treatments, and AUX, CK, and BR, including these were additionally regulated in the proton-beam treatments. In particular, the CK-related genes were only up-regulated in the proton-beam treatment. This hormone is known to be involved in cell division, leaf senescence, and nitrogen metabolism in plants [38]. The up-regulated cytokinin oxidase 5 is an enzyme that breaks down the hormone. The plant with increased cytokinin oxidase activity promoted root development and delayed shoot development [39]. When corn was treated to cold stress, the activity of this enzyme was increased at the root and shoot tips, and this enzyme was suggested to play a role in controlling the growth and development of organs [40]. Hormonal defense responses to adverse environmental conditions are known to rely on the complex crosstalk of signaling pathways rather than individual roles [35,37]. There is very little information on the role of hormones in responding to ionizing radiation. However, it is clear that compared to gamma-ray treatment, cowpea plants appear to require more complex hormonal mechanisms, including AUX, CK, and BR, to optimize development and growth by proton beam treatments.
TFs play an important role in the complex signaling processes, from the recognition of a stimulus to the expression of genes that respond to it [41]. These are classified according to the type of DNA-binding domains, and regulate metabolism by promoting or inhibiting specific gene expression in plants [42]. In this study, more diverse types of TFs were controlled by proton-beam treatments than by gamma-ray treatments. The Dof, ERF, and MADS were down-regulated in both sources, and most TFs tended to be up-regulated. In particular, bZIP, NAC, HSF, WRKY, and MYB were specifically up-regulated by the proton beams. However, these tended to be down-regulated overall by both ionizing radiation in rice irradiated with gamma-rays and carbon ion beams [25]. The AP2/EREBP was up-regulated by 800 Gy of gamma-rays, and C2H2 zinc finger and WRKY tended to be upregulated at 100 Gy and down-regulated at 800 Gy in Arabidopsis plants [17]. In other words, plant metabolism is complexly regulated by varying the combination and concentration of TFs according to the ionizing radiation and the type of crop. WRKY controls the seed size, pathogen defense, senescence, and trichome development [43,44], and bZIP is involved in pathogen defense, light and stress signaling, seed maturation, and flower development [42]. DOF regulates metabolism in response to the environment, as well as tissue differentiation and seed development [45]. Additionally, gene expression regulation of TFs can ultimately lead to increased resistance to abiotic stress. The overexpression of the NAC gene significantly increased the resistance to drought and salt in rice without a growth delay [46]. Furthermore, ERF proteins, a subfamily of the APETALA2 (AP2)/ethyleneresponsive-element-binding protein (EREBP), made tobacco plants resistant to drought and osmotic stress without having a significant effect on growth and development [47]. However, overexpressed DOF leads to growth defects [48]. Summarizing these results, it is proposed that TFs perform overlapping functions, such as plant development and stress tolerance, by controlling the combination and concentration in response to the two ionizing radiations and the specific role of individual TFs is not clear. However, it is clear that these TFs reacted more to proton-beams than gamma rays.
Notably, in the case of the photosynthesis target genes, they were all up-regulated only in 200 Gy proton-beam treatment. Most were involved in carbon fixation in photosynthetic organism's pathway. This result was similar in that the low-dose gamma-rays increased the photosynthetic efficiency of Arabidopsis thaliana, and high-dose gamma-rays did not differ from those of the control group [49]. However, N+ -beam implantation induced biological damage through downregulation of photosynthesis-related genes in plants [50]. Our results indicated that gamma-ray and proton-beam irradiation in the cowpeas did not negatively affect the photosynthesis system. Rather, it implied that certain doses had the ability to generate more energy to respond or adapt to the stimulation of ionizing radiation.
Stresses, such as ionizing radiation, soil salinity, and drought, cause damage to plants, including oxidative stress caused by ROS, and affect plant physiology [51,52]. Plants activate defense mechanisms to protect themselves from this damage [52]. In ROS-related genes, glutathione S-transferase TAU and peroxidase superfamily proteins are induced by biotic and abiotic stress and play an important role in resistance to oxidative stress [53][54][55]. These were significantly up-regulated in the proton-beam treatments compared to the gamma-ray treatments. In other words, it is predicted that cowpeas were subjected to more oxidative stress by the proton-beam than by the gamma-ray. Similarly, the MDA and antioxidant enzyme activity of cowpeas irradiated with proton beam was higher than that of gamma rays [56]. The defense-related genes tend to be the same as the ROS-related genes. These include cation exchangers and calreticulin, which play a role in resistance to biotic and abiotic stress and plant protection against oxidative stress [57,58], and beta-1,3glucanase, which retards fungal growth and reducing fruit spoilage [59].
When these results are summarized, it was concluded that the proton beam was a greater stressor than were gamma rays, which induces more disruption of normal protein production and accumulation of ROS. On the other hand, because many genes are involved in energy metabolism and environmental adaptation pathways were overexpressed in the proton-beam treatments, it is suggested that damage may have been minimized in terms of plant growth and development, and rather positive effects may have occurred. Additionally, a potential ability to induce positive changes in agricultural traits, such as flowering time, seed size, and resistance to biological and non-biological stress, can be proposed by the various genetic variations that respond to proton-beams.

Plant Materials and Radiation Treatments
A cowpea (Vigna unguiculata L. Walp) cultivar (Okdang) obtained from the Jeollanamdo Agricultural Research and Extension Services (JARES, Naju, Korea) was used in this study. The Okdang cultivar has an erect plant type with an intermediate plant habit and a high lodging resistance [60]. Okdang was irradiated with two different types of radiations, proton-beams and gamma-rays.

Proton-Beam Irradiation
The seeds were irradiated with a 57 MeV proton-beam at 100 MeV using the proton linear accelerator at the Korea Multi-purpose Accelerator Complex (KOMAC) in Gyeongju, Korea. Seeds of the Okdang cultivar were exposed to three different doses of proton beams (100, 200, and 300 Gy). A total of 200 seeds were used for each treatment.

Gamma-Ray Irradiation
Gamma-ray irradiation was conducted using the low-level irradiation facility containing 60 CO as a source at the Korea Atomic Energy Research Institute. The irradiation doses were the same as the proton-beam irradiation. A total 200 seeds were used for each treatment.

RNA Extraction
At 4 weeks after planting, young leaves from 21 plants were pooled for RNA extraction. Three independent pooled samples were prepared for each treatment. Isolation of total RNA was performed according to the manufacturer's instructions for the TRIzol reagent (Invitrogen, Carlsbad, CA, USA). RNA was quantified based on absorbance at 260 nm measured with a Nanodrop 2000 spectrophotometer (Thermo Scientific, Waltham, MA, USA). After the extracted RNA was stained with Dyne LoadingSTAR (DYNEBIO Inc., Seongnam, Korea), the integrity of the sample was confirmed using 1.5% agarose gel electrophoresis.

cDNA Library Construction and Massively Parallel Sequencing
RNA-Seq paired end libraries were prepared using the Illumina TruSeq RNA Sample Preparation Kit v2 (catalog #RS-122-2001, Illumina, San Diego, CA, USA). With total RNA, mRNA was purified using poly (A) selection or rRNA depleted, then RNA was chemically fragmented and converted into single-stranded cDNA using random hexamer priming. Next, the second strand was generated to create double-stranded cDNA. A library was constructed with blunt-end cDNA fragments from ds-cDNA. Then, A-base was added to the blunt-end ready them for ligation of sequencing adapters. After the size selection of ligates, the ligated cDNA fragments that contained adapter sequences were enhanced via PCR using adapter specific primers. The library was quantified with a KAPA library quantification kit (KK4854, Kapa Biosystems, Cape Town, South Africa) following the manufacturer's instructions. Each library was sequenced on the Illumina Hiseq2000 platform.

Preprocessing and Short Read Mapping
Sequence data, for which the quality of bp was indicated by Q ≥ 20, were extracted by SolexaQA [61]. Trimming resulted in reads with a mean length of 90.28 bp across all samples and a minimum length of 25 bp. Trimmed reads were mapped using the RNAseq mapping algorithm implemented in bowtie2 (v2.1.0) software [62] to the reference transcripts of Vigna unguiculata (v1.1) downloaded from the Phytozome database (http: //phytozome.jgi.doe.gov/ (accessed on 27 February 2021)), allowing all aligning with a maximum of two mismatches. The number of mapped clean reads for each gene was counted and then normalized with DESeq package in R [63] to avoid bias because of different sequencing amounts.

Identification of Differentially Expressed Genes (DEGs)
DEGs were identified by a ≥two-fold change in the number of mapped reads, a binomial test with a false discovery rate (FDR) ≤0.01, and a read count ≥ 100 between samples. The FDR was applied to identify the threshold p-value for multiple tests and was calculated using DESeq. Correlation analysis and hierarchical clustering was performed to group the genes according to patterns of expression using the AMAP library in R [64].

Functional Enrichment Analysis
To functionally annotate each DEG gene list, GO (Gene Ontology) analysis was conducted based on the sequence similarity (e-value cut off ≤ 1 × 10 −10 , best hits) of proteins in the Gene Ontology database [65]. The significance level was set to 0.05 using Blast2GO [66] and classified into functional categories: Biological Process (BP), Cellular Component (CC), and Molecular Function (MF). Additionally, the KEGG pathway was studied using the sequence similarity (e-value cut off ≤ 1 × 10 −10 , best hits) with the Vigna radiata protein in the KEGG database [67].

Target Gene Analysis
To obtain candidate genes related to ROS and hormones in cowpeas, protein sequences of Arabidopsis were used as a query, and 42,287 cowpea transcripts as a subject were compared and analyzed using the BLASTX (e-value ≤ 1 × 10 −10 , identity ≥ 50%). Repair system and photosynthesis-related genes were summarized by using the mung bean (Vigna radiata) gene information in the KEGG database. The comparative analysis used a BLAST tool and determined a gene aligned to the gene through a filter (e-value ≤ 1 × 10 −10 , identity ≥ 50%). The information of putative candidates of transcription factor (TFs) in cowpeas was obtained from PlantTFDB v3.0 [68].

Quantitative Reverse-Transcription PCR Validation of DEGs
For the validation analysis of mRNA expression levels quantitative real-time reverse transcription PCR (qRT-PCR) was conducted. The first-strand cDNA was synthesized from DNase-treated total RNA following the instructions of the SuperScript™ III First-Strand Synthesis SuperMix (Invitrogen, Carlsbad, CA, USA). Gene-specific primers were designed based on the nucleotide sequences of the selected DEGs for qRT-PCR analysis using Primer3 (v2.3.5) software (Table S15). The qRT-PCR were performed with a Bio-Rad iQ™ SYBR Green Supermix Kit (Invitrogen) using a StepOne Real-Time PCR System (Applied Biosystems, Foster City, CA, USA). The reaction mixtures (total volume of 20 µL) contained 100 ng of cDNA, each primer at 300 µM, 6 µL of ddH 2 O, and 10 µL of Bio-Rad iQ™ SYBR Green Supermix. The PCR conditions were set as follows: 95 • C for 3 min; 95 • C for 15 s, 60 • C for 60 s for 40 cycles; and the melt curve analysis was performed to confirm the absence of multiple products and the formation of primer dimers. The relative gene expression values were calculated by the 2 − Ct method [69]. The elongation factor 1β (EF-1β) gene was chosen as a reference gene for normalization of target gene expression in cowpeas [70]. For each of the three biological replicates, two technical replicates were included. For statistical analysis, the results were subjected to analysis of variance (ANOVA) using the SPSS version 25 software package (IBM, Armonk, NY, USA). Separation of means was performed using the Duncan's multiple range test (DMRT) with significance at p < 0.05.

Conclusions
After irradiation of three different doses of gamma rays and proton beams to the cowpea seeds, transcriptional variations were investigated. The numbers of DEGs were significantly higher in the proton beams when compared to gamma ray treatments. These DEGs were involved in a wide variety of biological functions, such as plant development, abiotic stress responses, signaling, and plant secondary metabolism. However, there were no common DEGs in all treatments for both sources. As a result of the GO and KEGG analysis, the two treatments showed a different spectrum of terms and pathways, and a difference in the number related genes. In other words, the gene expressions of cowpeas in response to the two ionizing radiations were highly diverse and complex, and their specific mechanisms were different. These results provided insights into genes that were regulated by gamma rays and proton beams and their functional relationships.
Supplementary Materials: The following are available online at https://www.mdpi.com/2223-7 747/10/3/567/s1, Figure S1: Quantitative Reverse-Transcription PCR validation of differently expressed genes (DEGs). Blue and orange bar graph represent relative expression level from RNA-seq and qRT-PCR, respectively. Error bars represent the SE for three independent replicates, Table S1: Read statistics for RNA sequencing of cowpeas exposed to gamma-rays and proton beams, Table S2: The information and functional analysis results of differentially expressed genes (DEGs) in cowpeas exposed to 100 Gy gamma-ray, Table S3: The information and functional analysis results of differentially expressed genes (DEGs) in cowpeas exposed to 200 Gy gamma-ray, Table S4: The information and functional analysis results of differentially expressed genes (DEGs) in cowpeas exposed to 300 Gy gamma-ray, Table S5: The information and functional analysis results of differentially expressed genes (DEGs) in cowpeas exposed to 100 Gy proton-beam, Table S6: The information and functional analysis results of differentially expressed genes (DEGs) in cowpeas exposed to 200 Gy proton-beam, Table S7: The information and functional analysis results of differentially expressed genes (DEGs) in cowpeas exposed to 300 Gy proton-beam, Table S8: The differently expressed genes (DEGs) list for each area of Venn diagram created using up-regulated DEGs in gamma-ray treatments, Table S9: The differently expressed genes (DEGs) list for each area of Venn diagram created using down-regulated DEGs in gamma-ray treatments, Table S10: The differently expressed genes (DEGs) list for each area of Venn diagram created using up-regulated DEGs in proton-beam treatments, Table S11: The differently expressed genes (DEGs) list for each area of Venn diagram created using down-regulated DEGs in proton-beam treatments Table S12: The differently expressed genes (DEGs) list for each area of Venn diagram created using DEGs in gamma-ray and proton-beam treatments, Table S13: The comprehensive results of clustering analysis performed using differently expressed genes (DEGs) in gamma-ray treatments, Table S14: The comprehensive results of clustering analysis performed using differently expressed genes (DEGs) in proton-beam treatments, Table S15: Primers used for qRT-PCR analysis for validation of differently expressed genes (DEGs).