Multi-Omics Analysis of Vicia cracca Responses to Chronic Radiation Exposure in the Chernobyl Exclusion Zone

1. Introduction

Radioactive contamination of the environment can have detrimental effects on living organisms, including plants. As sessile organisms, plants must cope with radioactive contamination through complex response pathways and metabolic adjustments. The effects of ionising radiation on plants can vary depending on the dose and duration of exposure and the characteristics of a plant species [1,2,3,4]. At low doses, radiation can stimulate plant growth and enhance photosynthesis, while at higher doses, it can lead to DNA damage, chromosomal aberrations, and mutations, ultimately affecting plant growth and reproduction [1]. Plants can undergo epigenetic changes in response to ionising radiation, leading to alterations in gene expression and phenotype [2,3]. The study of the biological responses of plants to chronic radiation exposure in radioactively contaminated areas is important for understanding the fundamental mechanisms of plant radiobiological responses and developing effective strategies to remediate radioactively contaminated sites.

The relevance of plant research in radioactively contaminated areas is twofold. First, plants can be used as bioindicators of radioactive contamination in the environment [5]. The morphological and physiological changes in plants exposed to ionising radiation can be used to monitor the extent of contamination and assess the effectiveness of remediation strategies. Second, plants can be used in phytoremediation, a process in which plants remove, stabilise, or detoxify radioactive elements from the soil or water [6,7,8]. Phytoremediation has several advantages over traditional remediation methods, such as being cost-effective, sustainable, and environmentally friendly [9]. However, the success of phytoremediation depends on several factors, including the plant species used, the level and type of contamination, and the local environmental conditions [9].

The genus Vicia belongs to the Leguminosae family, commonly known as vetches, and includes several plant species of agricultural importance. The genus comprises nitrogen-fixing plants with great potential for human and animal nutrition [10]. As green manure, Vicia can improve soil structure and plant nutrition [10] and mitigate climate change consequences as a substitute for synthetic fertilisers [11]. It harbours a large community of rhizospheric microorganisms, which can help to increase phytoremediation efficiency and make a remarkable difference in the ecological restoration of polluted soils due to Vicia dual role as cover crops and phytoremediator plants [12]. In recent years, several studies have focused on the use of Vicia species for phytoremediation: V. sativa for phenols [13], Cd [14,15], and polychlorinated biphenyl [16]; V. faba for Cd [17,18], Pu [17] [18], oil [19], and Sr [20]; V. villosa for As [21]; and V. cracca as a heavy metal accumulator for Cu, Fe, Zn, and Pb [22].

V. cracca is a hardy, fast-growing plant common in the Chernobyl Exclusion Zone (CEZ). In areas contaminated by uranium mill tailings and radium production wastes, these plants showed a significant increase in chromosome aberrations and a decreased survival rate of sprouts [23], outlining this species as a sensitive indicator of radiation exposure. Therefore, understanding the biological responses of V. cracca to ionising radiation is essential for its role as a bioindicator and for the possible use of the Vicia genus for phytoremediation purposes.

In recent years, multi-omics approaches have emerged as powerful tools for studying the responses of plants to environmental stresses, including pollution [24]. These approaches integrate data from different "omics" levels, such as genomics, transcriptomics, proteomics, and metabolomics, in order to comprehensively understand the biological processes and molecular mechanisms of plant responses to pollution. In the context of natural areas, multi-omics approaches can be used to study the responses of plants to a variety of pollutants [24,25,26]. This approach helps identify key molecular pathways and biomarkers that can be used to monitor and mitigate the effects of pollution on plants and the environment.

The current article presents an integrated multi-omics analysis of V. cracca plants inhabiting radioactively contaminated areas in the CEZ. This approach can provide a complete picture of the molecular changes in plants experiencing chronic low-dose radiation exposure and help identify key molecular pathways and potential biomarkers. Studies on V. cracca using the multi-omics approach can therefore contribute to a better understanding of the biological responses of plants to ionising radiation and the development of effective strategies for phytoremediation of radioactively contaminated areas.

2. Results

2.1. Environmental conditions in experimental plots

Vicia cracca leaves were sampled in the Polesye State Radiation-Ecological Reserve (Republic of Belarus) at four plots with different levels of radioactive contamination in June 2021. V. cracca leaves were sampled at Lomysh (Lom, L) and Babchin (Bab, B) as reference plots outside the CEZ and at two radioactively contaminated plots, Kulazhin (Kul, K) and Masany (Mas, M) (Figure 1).

The experimental plots had significantly different levels of radionuclide contamination (Table 1, Table S1). The contaminated plot Masany was characterised by a high ambient dose rate of 3.2 μSv x h^-1. At the same time, the highest specific activity of ¹³⁷Cs (36 990 Bq × kg⁻¹) was observed in the plot Kulazhin. Specific activities of ⁹⁰Sr were similar at both contaminated plots, approximately 1 500 Bq × kg⁻¹ (Table 1, Table S1). Soil properties at the experimental plots were slightly different (Table 2): soils samples at all plots were slightly acidic, with low humus content (1–5%) and low variation in the content of important minerals, being typical for soils of this region. Contents of heavy metals in soils did not exceed permitted levels (Table S1).

Additionally, in three plots out of four, we received data on annual absorbed dose rates from thermoluminescent dosimeters. Annual absorbed dose rates were calculated as mGy x year^-1, reflected external γ-radiation exposure, and comprised 1.13 mGy for Babchin, 51.57 mGy for Kulazhin, and 27.51 mGy for Masany for 365 days.

2.2. Analysis of expression of selective genes in V. cracca populations

As a first approach, we performed targeted gene expression analysis of some genes of interest revealed in our previous studies on other plant species [25,26,27]. The products of these genes participate in redox processes, signal transduction, and photosynthesis. Several genes encoding chaperones and histones were included to test the earlier observations of the possible role of these molecules in responses to chronic radiation [25,26,27]. Interestingly, the studied genes were mostly upregulated at two contaminated plots compared to both controls (Figure 2). This observation includes NADPH-oxidase RBOH-F (Figure 2A), ascorbate peroxidase APX1 (Figure 2B), chlorophyll-binding protein CAB1 (Figure 2C), transcription factor HY5 (Figure 2D), histone H2B (Figure 2I), and chloroplast chaperonin CPN60А (Figure 2G). However, another chloroplast chaperonin CPN20 was downregulated at both contaminated plots (Figure 2G). The differential expression of V. cracca genes is shown in Figure 2. These results match patterns found earlier in other plant species in the Chernobyl exclusion zone [25,26,27]. While these gene functions seem to characterize a chronic ionizing radiation response profile, the specific role of each of these genes on the adaptive responses to chronic radiation remains to be investigated.

2.3. Whole-transcriptome analysis of V. cracca populations

In order to estimate the whole mRNA transcriptional profile of V. cracca populations in conditions of chronic radiation exposure, we performed RNA sequencing of leaf samples. V. cracca transcriptome was assembled de novo. After functional annotation, we revealed DEGs in samples from each radioactively contaminated plot compared to each control and performed the Gene Ontology analysis to functionally describe transcriptomic changes under chronic radiation exposure. Comparison of samples from K and M radioactively contaminated plots with each reference plot (B and L) revealed 4446 for M x B comparison, 4011 for K x B, 5426 for K x L, 5672 for K x B, and finally 4033 for L x B. To exclude natural heterogeneity and possible confounding effects, we used only those genes which met the following conditions for further analysis: (1) had unidirectional differential expression compared to both controls and (2) were not differentially expressed between the two controls themselves. Among them were 457 downregulated DEGs and 371 upregulated DEGs (Table S2, Figure 3).

GO terms for DEGs shared for both radioactively contaminated plots compared to both controls are provided in Table 3. Upregulated DEGs were associated with nucleosomes and protein dimerization, and downregulated DEGs were mainly photosynthesis-related.

For in-depth analysis, we identified the 25 most downregulated and 25 most upregulated DEGs (Table 4). Among the downregulated ones, there is an evident suppression of photosynthesis-related elements, specifically, the Calvin cycle, chlorophyll biogenesis, thiamine biosynthesis, and translation process. The induction of peroxidases, protein catabolism machinery, phytohormonal relays, and cell cycle-involved proteins is noted for upregulated genes. However, some elements of photosynthesis were found to be upregulated, which suggests an overall fine-tuning of the photosynthesis rate.

2.4. Non-targeted proteomic analysis of V. cracca populations

To reveal changes in the protein profile of irradiated samples, we performed a non-targeted proteomic analysis. One hundred and six proteins were successfully identified in the V. cracca samples (Table S3). Proteins whose abundance changed at least two-fold in at least one of the radiation x control comparisons are highlighted in Figure 4. Interestingly, among them are proteins related to photosynthesis, ribosome assembly, and several cytochromes. Some proteins were identified only in control samples (cytochrome f A0A0F6NLU, NAD(P)H-quinone oxidoreductase A0A0F6NLR2), while another cytochrome f A0A7T4XAP1 was revealed only in samples from radioactively contaminated plots.

In order to describe the possible functional interconnection between the identified proteins, we performed an analysis with the STRING tool (Figure 5), using as an input homologous proteins of Medicago truncatula. The hubs of translational responses reflect the interconnection of photosynthetic responses and ribosomal proteins, suggesting profound changes in energy metabolism and translation. These results overlap with the transcriptomic data, further supporting the notion that a strategy for coping with chronic ionizing radiation exposure involves the modulation of photosynthesis rate, possibly for adjusting toxic ROS levels which are produced during energetic processes involving electron transport chains, but also due to the water radiolysis [2, 25, 26, 28].

2.5. Non-targeted metabolomic analysis of V. cracca populations

To specifically describe plant molecular phenotype and to reveal the results of interactions of transcriptome and proteome, we performed non-targeted GC-MS analysis of metabolites in V. cracca leaves. The obtained metabolic profiles included around 750 compounds, of which 70% were identified. We filtered the metabolites with close QC values for all replicates for further analysis. Thus, 139 metabolites were selected for further statistical analysis, and their concentrations are presented in Table S4. The normalised data were used for the PCA analysis. The scores plot (Figure S1) shows five principal components (PC1-PC5) explaining 74.2% of variations from the sixteen samples of V. cracca. According to PC1, which explains 24.1% of the variance, the metabolite profiles of samples from control plots were grouped separately from samples of contaminated plots, which overlapped. Supervised PLS-DA analysis (Figure 6) showed that samples from radioactively contaminated plots K and M were separated from control samples and each other by Component 1 (Figure 6A). The reliability of the PLS-DA model was confirmed by R² (0.926) and Q2 (0.754) values obtained as a result of cross-validation. According to VIP scores, the greatest contribution to the separation was O-methyl-inositol, inositol-3-phosphate, and pinitol (Figure 6B).

To visualise the metabolomic changes in control and chronically irradiated V. cracca plants, the top 50 metabolites were analysed using Heatmap and hierarchical cluster analysis. An obvious separation was observed between control and chronically irradiated plants, where four subclasses corresponded to four studied plots (Figure 6C).

One-way ANOVA with Fisher’s HSD at the 5 % confidence level was used to identify metabolites with different concentrations on radioactively contaminated plots. Data from 2 control plots, L and B, were combined and used as one group. As a result, 20 metabolites were identified in V. cracca, which were affected by chronic radiation exposure (Table 5).

A total of 139 metabolites was used for volcano plot analysis to understand how metabolites change at each experimental plot relative to the control (Figure 7). Specifically, K had 10 significantly downregulated metabolites and 19 upregulated compared to control L+B (Figure 7A), while M had 11 downregulated metabolites and 15 upregulated (Figure 7B) with biological and statistical cut-offs for fold change |FC| > 2 and p-value of < 0.05. Metabolites related to the citrate cycle, galactose metabolism, and glyoxylate and dicarboxylate metabolisms such as galactinol, raffinose, isomaltose, aconitate, and citrate were downregulated at both contaminated plots K and M compared to the control.

Kyoto Encyclopedia of Genes and Genomes (KEGG) identification numbers were used for pathway analyses of significantly changed metabolites (Table 5) in MetaboAnalyst 5.0 against the Arabidopsis thaliana library as a reference. Metabolic pathways altered considerably by chronic radiation exposure were revealed (Figure 7C), which included galactose metabolism (impact value: 0.27, FDR = 0.002), tyrosine metabolism (impact value: 0.22, FDR = 0.007), and citrate cycle (TCA cycle) (impact value: 0.15, FDR = 0.01) (Figure 7C).

2.6. Multi-omics data integration

The PaintOmics4 tool was used to identify pathways enriched or underrepresented in contaminated plots on different functional levels for a holistic view of the data. Among them were general pathways of carbon (Figure S2) and fatty acid catabolism (Figure S3), as well as N-glycan biosynthesis (Figure S4). The multi-omics analysis confirmed changes in the citrate cycle, where fumarate, 2-oxoglutarate, and malonate accumulation was observed, and citrate concentrations decreased (Figure 8), and in pyruvate metabolism (Figure 9), where accumulation of pyruvate was noted accompanied by downregulation of glycolysis in general (Figure 9, Figure S2). Profound changes in photosynthesis accompanied the changes in energy metabolism (Figure 10). Among more specific pathways, we observed upregulation of the glyoxylate cycle (Figure S5), downregulation of biotin metabolism (Figure S6), and upregulation of certain branches of propanoate metabolism, which led to the β-alanine synthesis (Figure S7). Structural modifications of ribosomes were also evident (Figure S8), accompanied by significant changes in RNA polymerase subunits (Figure S9). Altogether, the results suggest that the modulation of energy production is indeed a central aspect of plants growing under chronic radiation exposure.

3. Discussion

3.1. Carbon metabolism and photosynthesis.

Cells of V. cracca exposed to chronic ionising radiation are changing their carbon cycle metabolism, including the upregulation of pyruvate (Figure 7, 9) and fumarate (Figure 7, 8) biosynthetic processes, which points to energetic and plastic rearrangements. It can be connected with asparagine synthase decreased protein abundance (Figure 4), as asparagine synthesis-deficient A. thaliana mutants showed an accumulation of alanine, GABA, pyruvate, and fumarate, indicating an alanine formation from pyruvate through the GABA shunt to consume excess ammonium in the absence of asparagine synthesis [30]. Indeed, alanine, aspartate, and glutamate metabolism are also enriched in our data (Figure 7). On the other hand, the decreased citrate concentrations (Figure 7) coupled with the upregulation of transcripts involved in its biosynthesis (Figure 8) points to several bypass possibilities of TCA due to redox changes in the cell. Organic acids represent the stored pools of fixed carbon. When the redox level in the cell increases, the TCA cycle in mitochondria is transformed into a partial cycle supplying citrate for the synthesis of 2-oxoglutarate and glutamate (citrate valve). At the same time, malate is accumulated and participates in the redox balance in different cell compartments (via the malate valve) [31,32]. Indeed, multi-omics analysis revealed higher 2-oxoglutarate concentrations and malate accumulation (Figure 8), which may support the activation of citrate and malate valves. The upregulation of the glyoxylate cycle intermediates, such as D-glycerate and glycolate (Figure S5), also suggests that the cells are trying to bypass some steps of the TCA cycle to generate key metabolic intermediates. It could be a response to the decreased availability of other metabolic intermediates, such as citrate [31].

Fatty acids are major components of cell membranes, and changes in their biosynthesis could alter the fluidity and permeability of the membranes. The downregulation of pathways leading to fatty acid biosynthesis (Figure S3) suggests that the cells may be reducing their production of fatty acids in response to chronic radiation exposure. The downregulation of fatty acid biosynthesis pathways may have implications for the overall lipid metabolism and membrane composition of plant cells [51]. The decrease in the production of propanoyl-CoA from 2-oxobutanoate (Figure S7) also confirms lipid metabolism downregulation in irradiated plants. The increase in the production of malonate semialdehyde from 3-hydroxy-propionyl-CoA (Figure S7) suggests that the cells may be relying more heavily on the degradation of odd-chain fatty acids or amino acids for energy production [33].

Reallocation of carbon can be seen as profound changes in galactose metabolism (Figure 7), indicated by the decrease of raffinose family oligosaccharides (RFOs) concentrations. RFOs are involved in antioxidant, membrane-stabilising, and signalling pathways [34]. However, the reduction of galactinol and raffinose concentrations may rather be associated with their alleged role in providing carbon skeletons to lipid synthesis or thylakoid membrane biogenesis [35]. Further depletion of some carbohydrates was also confirmed by the downregulation of various N-glycans (Figure S4), which are involved in N-glycosylation and participate in multiple processes, including development, protein folding, photosynthesis, and phytohormone homeostasis [36]. N-glycans and raffinose family oligosaccharides interact with lectins in response to pathogens to initiate defence reactions [37]. It is plausible to suggest that the downregulation of raffinose and N-glycan metabolism may be partially attributed to the attempt to avoid overreaction to continuous damage of cell wall structure by ionising radiation, which can be perceived as a pathogen attack and promote defence responses.

Changes in carbon availability can be related to decreased capacities of carbon assimilation. The chloroplasts of plants are particularly sensitive to ionising radiation due to the susceptibility of their light-absorbing pigments and genetic material to the negative effects of reactive oxygen species (ROS) [3839,40]. Several photosynthesis-related genes were downregulated (Table 4). The changes in the antenna complex LHCII and photosystems I and II (Figure 10) were observed, suggesting that the photosynthetic machinery may undergo alterations in response to chronic radiation exposure. Downregulated DEGs included ISPE_SOLLC, PUR5_VIGUN, PGL1B_ARATH, FTSI1_ARATH, CHLM_ARATH, and PGKH_WHEAT (Table 4), which are associated with chloroplast development and chlorophyll biosynthesis [41,42,43,44,45]. A decreased expression can be either associated with redox damage of the electron transport chain or developmental changes under chronic irradiation conditions. The reduced synthesis of cytochrome b6/f (Figure 4, 10B) and changes in electron transport chain components suggest a decrease in the efficiency of photosynthesis. The changes in Fd, FNP, and PC (Figure 10B) also imply electron transport chain dynamics alterations, which may impact the overall ATP production. The upregulation of ATPase components may indicate an increase in the cellular demand for ATP (Figure 10B), possibly due to the decreased efficiency of other energy-generating pathways like glycolysis, which is evidently downregulated (Figure 9). However, in the targeted gene expression analysis, we observed the upregulation of CAB1 homologue (Figure 2), which is associated with chlorophyll accumulation and stability of chloroplast membranes [46]. Expression of transcription factor HY5, which regulates, particularly, CAB1 expression [47] and is a key factor of chlorophyll accumulation, was increased in all contaminated plots (Figure 2). Induction of these genes is aligned with our previous results [25, 27] and points to the important role of HY5 in responses to chronic radiation and to the possible disruption of circadian rhytmes of irradiated plants [47,48,49].

Overall, these changes in energy metabolism suggest that the cells undergo significant alterations in response to chronic radiation exposure, possibly redirecting metabolic intermediates of the TCA cycle for the synthesis of protective compounds and relying more on fatty and amino acid catabolism for energy production than plants from the reference plots.

3.2. Chronic radiation stress response: DNA damage and changes in the redox status.

Chronic radiation exposure is a mutagenic factor capable of inducing DNA single- and double-strand breaks. Our previous research, however, revealed a very subtle DNA repair machinery response in plants growing under chronic radiation, such as Pinus sylvestris [25] and Capsella bursa-pastoris [26]. Deep sequencing of V. cracca transcriptome revealed upregulated DEGs involved in DNA repair and cell cycle control (Table 4). Among them are topoisomerase 3-alpha TOP3A_ARATH, involved in homologous recombination, and cyclin CCNB1_MEDSA. However, the canonical DNA damage response elements (such as ATR, ATM, and SOG1 [2]) were not identified, which may point to the involvement of identified DEGs in replication maintenance rather than in radiation-caused strand break repair. Nevertheless, if chronic radiation exposure affects DNA replication efficiency remains to be investigated. Upregulation of histones (Figure 2, Table 4, Table S2) and nucleosomal processes (Table 3) may confirm our previous hypothesis of the role of DNA packaging in coping with chronic DNA damage rather than keeping energy-demanding systems of DNA damage recognition constantly upregulated [25, 26]. Homologous recombination is a part of the response to DNA double-strand breaks [2], but other components of DNA damage response signalling [50] were not deregulated in our data (Table S2). Furthermore, a potential inducer of apoptosis CEP1_ARATH [51] was strongly upregulated (Table 4). Therefore, the data point to at least moderate DNA damage in plants from radioactively contaminated plots.

Chronic damage of macromolecules implies the need to protect DNA and protein integrity and conformation. Based on our previous findings of the involvement of histones and chaperones in responses to chronic radiation exposure [25,26], we performed targeted expression analyses of several chaperones (Figure 2). The upregulation of chloroplast chaperonin CPN60А homologue may be connected to decreased RuBisCO abundance (Figure 4). Rubisco deactivation may be a protective strategy under stress conditions, while CPN60A, which is involved in carboxylation, participates in photosynthesis maintenance under stress conditions, being upregulated, particularly under heat stress [52]. However, the CPN20 homologue encoding another chloroplast chaperonin was significantly downregulated in both radioactively contaminated plots (Figure 2). CPN20 negatively regulates certain branches of ABA signalling [53,54]; therefore, its downregulation may be connected to CPN20 signalling function rather than protein folding.

Finally, the redox imbalance in V. cracca is confirmed not only by TCA bypassing and photosynthesis reduction. The upregulation of peroxidases (Table 4, Figure 2) suggests increased ROS concentrations. Downregulation of biotin metabolism (Figure S6) may lead to biotin deficiency, which, in turn, is associated with H₂O₂ accumulation [55]. In turn, biotin deficiency can contribute to observed changes in carbon metabolism since biotin is a cofactor for some carboxylases, decarboxylases and transcarboxylases participating in fatty acid and carbohydrate metabolism [56]. Under radiation-induced ROS imbalance, accumulating secondary metabolites with antioxidant properties [1] is a promising coping strategy for chronic oxidative stress.

3.3. Secondary metabolites.

Indeed, several DEGs involved in the biosynthesis of phenolic compounds, which are important antioxidants, were induced in response to ionising radiation (Table S2), including PALY_MEDSA (Table 4), phenylalanine ammonia-lyase. This enzyme catalyses the first synthesis step of various phenylpropanoid compounds [57] and is often induced by ionising radiation [58,59,60]. The upregulation of the HY5 transcription factor gene (Figure 2) can also induce anthocyanin biosynthesis [47]. Meanwhile, the downregulation of genes involved in terpenoid biosynthesis, such as ISPE_SOLLC and CHLP_ARATH (Table 4), potentially indicates a compensatory role in the biosynthesis of other metabolites. In V. cracca from contaminated plots, increased concentrations of metabolites involved in tyrosine metabolism were observed (Figure 7). Tyrosine is a precursor of metabolites that perform antioxidant functions (tocopherols, rosmarinic acid), electron transfer in the electron transport chain (plastoquinone and ubiquinone), and plant protection against pathogens (alkaloids and glycosides) [61].

Nitrogen reallocation is essential for nitrogen-containing secondary metabolites biosyntheses such as alkaloids, cyanogenic glycosides, glucosinolates, and nonprotein amino acids [62]. Besides general changes in amino acid metabolic pathways, we observed a significant increase in the allantoin concentration in the contaminated plot Masany (Table 5) with the highest ambient dose rate and activities of ¹³⁷Cs and ⁹⁰Sr in soil. Allantoin is an intermediate product of the purine catabolic pathway, which acts in nitrogen metabolism and transport in plants [63] and accumulates in plants during environmental challenges increasing stress tolerance [64]. In chronically irradiated V. cracca plants, an increase in the concentrations of nicotinate and its derivative trigonelline was revealed (Figure 7), which may be associated with oxidative damage to cells and perturbations in cellular energy metabolism [65]. Several other secondary metabolites had a higher concentration on radioactively contaminated plots (Figure 7), including compounds with antioxidant and protective properties, such as pinitol [66], trigonelline [67], β-amyrin [68], hydrangeifolin [69], and piceid [70]. Therefore, at least partially, certain general metabolic pathways can be suppressed due to carbon and nitrogen reallocation for protective compound synthesis under chronic irradiation.

3.4. Protein catabolism and RNA processing.

The changes in ribosomal protein synthesis (Figure S8) suggest a complex response to chronic ionising radiation, with some ribosomal proteins being upregulated and others being downregulated. Such differential regulation patterns can be related to tissue-specific responses [71] or be an example of multidirectional ribosomal protein responses, as was shown for other stressors [72,73]. For instance, ribosomal proteins of large ribosomal subunits are highly responsive to stress and signalling molecules, indicating that they may have roles in stress responses besides housekeeping function [72]. Extraribosomal functions of some ribosomal plant proteins include microRNA biogenesis, anti-virus defence, and plant immunity, which may be associated with the widespread duplication of ribosomal protein genes in plants [74].

The changes in RNA polymerase biosynthesis are also multidirectional (Figure S9). Eukaryotic RNA polymerase II components are mostly upregulated (Figure S9), while this enzyme complex is responsible for the transcription of protein-coding genes, long non-coding RNAs, and small RNAs [75]. The downregulation of RNA polymerases I and III components (Figure S9) may suggest a decrease in the transcription of rRNA and tRNA genes [75]. The processing of rRNA is deregulated on the transcriptional level (Table 4). Overall, the changes in ribosomal protein synthesis and RNA polymerases suggest that the cells undergo significant alterations in their gene expression and translational patterns in response to chronic radiation exposure.

Plant growth and development are directly linked to ribosome synthesis and protein translation. Still, the precise molecular mechanisms and signals that connect cellular and environmental conditions with rRNA synthesis and ribosome assembly remain to be uncovered [76]. Yet, the ribosomal responses of V. cracca to chronic irradiation may be of particular interest since, lately, ribosomal proteins have been suggested as a valuable resource for manipulating the stress tolerance of crops [72].

4. Materials and Methods

4.1. Experimental plots and soil sampling.

Sampling was conducted at the Polesye State Radiation-Ecological Reserve (Republic of Belarus) at four plots with different levels of radioactive contamination in June 2021. V. cracca leaves were sampled at Lomysh (Lom, L) and Babchin (Bab, B) as reference plots outside the CEZ and at two radioactively contaminated plots, Kulazhin (Kul, K) and Masany (Mas, M) (Figure 1).

At each experimental plot, the ambient dose rate was determined using a dosimeter-radiometer MKS-02SA1 (SNIIP, Russia) at a 1 m from the surface, and α- and β-particle flux densities (min⁻¹ × cm⁻²) were determined 0.2 cm above the ground. Soil samples at each plot were collected by the envelope method from a rectangular area (1 × 2 m) inside the plot at a depth of up to 15 cm, and then the soil was mixed into a pooled sample. The activity concentrations of ¹³⁷Cs in the soil samples were measured with a γ-spectrometer CANBERRA (USA). We used the radiochemical method to determine ⁹⁰Sr specific activity in soil sampling. To estimate the annual absorbed dose of external γ-radiation at each experimental plot, we installed a thermoluminescent dosimeter DTL-02 (NPP Doza, Russia). The dosimeters were covered from weather exposure by a section of plastic tube installed over them, and laid on the ground inside the perimeter of a sample plot.

At each plot, soil properties and levels of heavy metal contamination were assessed, including pH, hydrolytic activity, cation exchange capacity, K, P, Ca, Mg, Na, and humus content in accordance with the ISO standard for soil quality [77]. Total concentrations of heavy metals (Cd, Cu, Co, Ni, Cr, Mn, Pb, Zn, As, and Mo) in the soil samples were measured using a plasma optical emission spectrometer (ICP-OES, Australia) in accordance with the ISO 11047 standard following treatments with a mixture of HNO₃, HCl, and HF [78].

4.2. Plant sampling.

V. cracca was sampled at the four plots by collecting pools of leaves from 5 to 7 plants (approximately 0.5 g total). Sampled plants visually belonged to the same developmental stages. The samples were immediately frozen in liquid nitrogen until analyses. For quantitative real-time PCR (qRT-PCR) gene expression analysis, we used 4 samples per plot. For transcriptome, proteome, and metabolome analyses, leaves V. cracca were sampled from the same plants for all three analyses to ensure correlation among all omics levels. For the transcriptomic analysis, we used 3 samples per plot. For proteomic analysis, we used 4 samples for each contaminated plot and 3 for each reference plot. For metabolomic analysis, 4 samples for each reference and radioactively contaminated plots were used. Sampling was performed around noon for all experimental plots (Table S1). The samples were collected under similar environmental conditions, and no rain occurred during the sampling period.

4.3. Gene expression by qRT-PCR.

Based on our previous research using other plant species growing in the CEZ [25,26,27] and data from other authors on plant responses to chronic radiation exposure [79,80], we selected several genes to study specific transcriptional responses of V. cracca to chronic irradiation (Table 6). The products of these genes are associated with the antioxidant system, signal transduction, and photosynthesis processes. The expression of genes encoding aquaporins was also studied, as these proteins maintain the redox balance of cytosol and hydrogen peroxide transport [81]. To test the hypothesis of the essential role of chaperones and histones in adaptive response to chronic irradiation, firstly raised in Pinus sylvestris L. [25] and Capsella bursa-pastoris L. [26], we included in the analysis several genes encoding chaperones and histones for V. cracca. CYP2 and ELF1A [82] were used as reference genes.

For primer design, cDNA sequences of A. thaliana were aligned against Fabaceae sequences using the NCBI BLAST tool. 5-10 sequences of homologous genes with the highest score were selected for each target gene. The sequences were processed in Сlustal Omega 1.2.2, revealing the conserved regions. Next, primers to the selected regions were constructed using Primer BLAST. The specificity of primer pairs was tested using qRT-PCR and 2% agarose gel separation. Primer sequences are provided in Table S5.

Total RNA from V. cracca leaves was isolated using GeneJET Plant RNA Purification Mini Kit (Thermo Fisher Scientific, USA). RNA concentration and quality were measured using NanoDrop OneC (Thermo Fisher Scientific, USA). 1 μg of total RNA was treated with DNase I (Thermo Fisher Scientific, USA). MMLV RT kit (Evrogen, Russia) was used to synthesise the first strand of cDNA. The resulting cDNA was diluted 1:10 in nuclease-free water and used as a template for the qRT-PCR.

qRT-PCR was performed using a DT-96 PCR machine (DNA-Technology, Russia) under the following thermal cycling conditions: initial cycle at 50 °C for 2 min, initial denaturation for 2 min at 95 °C, followed by 50 annealing/elongation cycles (95 °C for 15 s, 60 °C for 60 s). The volume of the reaction mixture for one reaction was 20 µl, which included 4 µl of the cDNA template, 2 µl of a 10 µM mixture of forward and reverse primers, 4 µl of the qRT-PCR mix-HS SYBR reaction mixture (Evrogen, Russia), 10 µl of nuclease-free water. To calculate changes in gene expression (RQ), we used the ∆∆Сt model [83]. A significant change in gene expression was considered a twofold increase or decrease in expression compared to both controls.

4.4. RNA sequencing and transcriptomic analysis.

Total RNA was isolated using the GeneJET RNA Purification Kit (Thermo Fisher Scientific, USA) with polyvinylpyrrolidone addition. The quality and purity of isolated RNA were assessed using NanoDrop OneC (Thermo Fisher Scientific, USA) and horizontal gel electrophoresis. cDNA synthesis, library preparation, and sequencing were provided by "Evrogen" (Moscow, Russia). Using TruSeq mRNA Stranded reagent kit (Illumina, USA), poly(A+)-fraction enrichment and random primer cDNA synthesis were performed in 12 RNA samples. The resulting cDNA was used to prepare libraries compatible with Illumina sequencing technology. The quality of the obtained libraries was checked using Fragment Analyzer (Agilent, USA). Quantitative analysis was performed using the qRT-PCR method. After quality control, the pool of cDNA libraries was sequenced using Illumina NovaSeq 6000 (150 bp, pair-ended). FASTQ files were obtained using bcl2fastq v2.20 Conversion Software (Illumina, USA). The format for recording the quality data string was Phred 33. As a result, 1 090 239 068 raw reads were obtained.

The transcriptome data analysis was performed following suggestions previously described in [84]. Quality assessment was performed using FastQC v.0.11.9 and MultiQC v.1.10 software. Trimmomatic v.0.40 software was used for trimming using the following parameters: ILLUMINACLIP:TruSeq3-PE.fa:2:30:10; HEADCROP:11; LEADING:3; TRAILING:3; MAXINFO:130:1; MINLEN:50. The trimmed reads were deposited in Sequence Read Archive under accession number PRJNA958217, and used for de novo transcriptome assembly using Trinity v.2.9.1 at the Galaxy public server https://usegalaxy.eu [85]. In silico normalisation of paired and unpaired reads was performed, and unpaired readings were attached to the "left" paired read. Contigs smaller than 200 bp were filtered from further analysis. Assembly completeness was tested in BUSCO v2 [86]. The fabales_odb10 data set corresponding to representatives of the Fabaceae family was used for evaluation. The result showed 89.4% similarity of the assembly with the group of orthologues and a high percentage of duplications (84.0%), which may indicate the possible tetraploidy of the samples since V. cracca is known to possess either diploid or autotetraploid cytotypes [87,88]. The alignment of paired reads to the assembled transcriptome showed a satisfactory result (above 60%). Open reading frames were predicted for the assembled contigs in TransDecoder v.5.5.0. The search for protein (blastp) and translation (blastx) motifs were performed on sequences obtained from the Swiss-Prot database. Protein domains were identified in HMMER v2.2 [89] and Pfam [90]. Predictions of signal peptides, rRNA transcripts, and transmembrane regions were performed in SignalP v4.1 [91], RNAMMER v1.2 [92], and TmHMM v2.0, respectively. Functional annotation was carried out in Trinotate v.3.2.2 [93]. and InterProScan 5.59_91.0 [94]. As a result, gene annotation of contigs and terms of gene ontology (GO) were obtained.

The quantification algorithm kallisto v.0.46 [95] and DESeq2 statistics package [96] were used to quantify gene expression. Comparisons were made between contaminated and control plots and two control plots (Babchin-Lomysh) to assess heterogeneity (Figure S10). Genes with log₂FC ≥ |2| and adjusted p-value ≥ 0.05 (Benjamini-Hochberg method) were considered as differentially expressed. DAVID Bioinformatics tool [97] was used to identify significant GO terms. Venn diagrams were created using the web tool http://bioinformatics.psb.ugent.be/webtools/Venn.

4.5. Protein isolation and proteomic analysis.

For protein isolation, at each control plot, we used 3 replicates and 4 replicates for each radioactively contaminated plot. Two buffers were used to isolate and purify protein from V. cracca leaves: (1) extraction buffer (EB): 0.2 M TRIS-HCl, pH 9.0; 0.2 M KCl; 0.025 M EGTA, pH 8.0; 0.035 M MgCl2; 1% (V/V) PTE; 1% (V/V) detergent blend (Brij-35, Tween 20, Triton X-100, Igepal CA 630); 0.005M DTT; 1x cocktail of protease inhibitors; (2) lysis buffer (LB): 4% SDS; 100 mM Tris/HCl, pH 7.4; 1x cocktail of protease inhibitors. 100 mg of tissue per sample was ground in liquid nitrogen and dissolved in 750 μL of buffer EB on ice for 20 minutes, pipetting until complete dissolution of plant tissue aggregates. The suspension was centrifuged at maximum speed for 10 minutes at 4 °C. The supernatant was collected and stored on ice. The precipitate was redissolved in 500 μL of LB, incubated for 3 minutes at 95 °C, and then sonicated for 5 minutes to destroy the membranes. After centrifugation, the supernatants were combined, and 10% TCA in acetone pre-cooled to ‒20 °C was added to precipitate the proteins. This solution was incubated overnight at ‒20 °C and then centrifuged for 10 minutes. The protein precipitate was then washed twice with cold acetone and dried. Before the analysis, the protein pellet was stored at ‒80 °C.

Skolkovo Institute of Science and Technology provided proteomic analysis services. The precipitated protein was dissolved in 200 µL of sample dilution buffer, placed in an ultrasonic bath for 1 min, and centrifuged for 10 min at 14 000 g and 10 °C. The supernatant was transferred to clean, labelled vials. To estimate the amount of protein by the BSA method, the following mix was prepared: 5 µL of the supernatant obtained during lysis, 25 µL of water, 1 mL of the BSA reagent, and 20 µL of a 4% copper sulphate solution. To prepare a blank sample, 30 µL of water, 1 mL of the BSA reagent, and 20 µL of a 4% copper sulphate solution were mixed. Similarly, calibrant samples were prepared (solutions of BSA with a concentration of 1 µg/µL to 10 µg/µL) with BSA reagent. The samples were shaken and incubated for 20 min at 56 °C with stirring on a thermomixer. The samples were cooled to room temperature and analysed on a CLARIOstar reader (BMG Labtech, Germany). For further analysis, samples were diluted in 50 mM TEAB to a final protein concentration of 2 µg/µL.

For protein hydrolysis, 10 µL of samples were transferred into 200 µL tubes, corresponding to 20 µg of protein. 10 µL of protein recovery mixture was added to the samples, shaken, rapidly centrifuged, then incubated for 30 min at 80 °C and cooled to room temperature. 100 µL 50 мМ TEAB and 2 μL of trypsin solution (trypsin:protein ratio 1:100) were added into the samples and incubated for 4 hours at 42 °C. Next, another 2 μL of trypsin solution was added (trypsin:protein ratio 1:50), incubated for 18 hours (overnight) at 37 °C, and then cooled to room temperature. Formic acid (2 μL) was added, shaken, and cooled for 10 minutes in a refrigerator at 4 °C, then centrifuged for 10 min at 14 000 g and 10 °C. The supernatant was transferred into glass inserts for HPLC-MS analysis and evaporated in a vacuum concentrator, then redissolved in 10 μl of 0.1% formic acid in water. Then samples were analysed by high-performance liquid chromatography with tandem mass spectrometry (HPLC-MS/MS) using UltiMate 3000 RSLCnano system coupled with Q Exactive HF hybrid quadrupole-Orbitrap mass spectrometer (Thermo Scientific, USA). Samples were analysed in three technical replicates. 2 μL per sample was separated on an Acclaim Pepmap C18 column (Thermo Scientific, USA) in a gradient elution mode. The gradient was formed by mobile phase A (0.1% formic acid in deionised water) and mobile phase B: (80% acetonitrile, 0.1% formic acid in deionised water). The column was washed with 2% mobile phase B for 1 min at a flow rate of 700 nL/min, then 1 min at 500 nL/min and 3 min at 400 nL/min. Further separation was carried out at a flow rate of 400 nL/min: the concentration of mobile phase B was linearly increased to 5% in 3 min, then the concentration of phase B was linearly increased to 30% in 70 min, then the concentration of phase B was linearly increased to 99% in 2 min, after a 5-minute wash at 99% buffer B, the concentration of buffer B was linearly reduced to the initial 2% over 3 minutes, and the column was washed with 2% mobile phase B for 8 minutes. The total duration of the analysis was 95 min.

Mass spectrometric analysis was performed on a Q-Exactive HF mass spectrometer in positive ionisation mode using a NESI source (Thermo Scientific, USA). The following parameters were set: emitter voltage 2.1 kV, capillary temperature 240 °C. Panoramic scanning was carried out in the mass range from 390 m/z to 1200 m/z at a resolution of 60 000. In tandem scanning, the measurement was carried out in the mass range from 200 m/z to the upper limit, which is determined automatically based on the mass of the precursor, but not more than 2 000 m/z, at a resolution of 15 000. Precursor ions were isolated in a window of ± 1 Da. The maximum number of ions allowed for isolation in the MS/MS mode was no more than 10, the minimum intensity of the precursor ion for tandem analysis was set as 50 000 units, and the normalised collision energy was 29. For tandem scanning, only ions from z = 2+ to z = 4+ in charging conditions. The maximum accumulation time for precursor ions was 50 ms, and for fragment ions, 150 ms. All measured precursors were dynamically excluded from the tandem MS/MS analysis for 10 s.

Proteins were identified using the MaxQuant v.1.6.17.0 software and the Andromeda search algorithm. Proteins were identified using the UniProt database with species indication. The following search parameters were set: the cleaving enzyme trypsin, the accuracy of mass determination of m/z peptides in MS ± 5 ppm, the accuracy of mass determination in MS/MS spectra ± 0.01 Da, and the possibility of skipping two trypsin cleavage sites. Proteins were considered reliably identified if at least two peptides were found. Oxidation of methionine, acetylation of the N-terminus of the protein, and modification of cysteine with chloroacetamide were considered possible and mandatory modifications of the peptides, respectively. The False Discovery Rate value of no more than 1.0% was used to validate comparisons of Peptide-Spectrum Matches spectra and peptides. Label-free protein quantification was based on iBAQ. The data were filtered to remove possible contaminants, and unreliably identified proteins and only proteins that occurred in two technical repetitions of the MS analysis were retained. We used only proteins with log2FC ≥ |1| for further functional evaluation compared to the control samples.

4.6. Metabolite extraction and metabolomic analysis.

Metabolite extraction, derivatisation, and GC-MS analysis were performed at Institute Jean-Pierre Bourgin (INRAE), as previously described [98]. The freeze-dried samples of V. cracca leaves were ground using metal beads and resuspended in 1 mL of frozen (-20 °C) water:acetonitrile:isopropanol (2:3:3) solution, containing ribitol at 4 mg.mL^-1 and extracted for 10 minutes at 4 °C by shaking at 1400 rpm in an Eppendorf thermomixer (Eppendorf, Germany). Insoluble material was removed by centrifugation at 20 000g for 5 minutes. 100 mL of supernatant were collected and dried for 4 hours in a Savant Speed-Vac (Thermo Fisher Scientific, USA) and stored at -80 °C. Three blank tubes underwent the same steps as the samples. Quality controls (QC) were prepared by pooling the same amount of each sample. They were injected at the beginning, the middle, and the end of the injection series to ensure the stability of the derivatisation. After drying for 1 hour after -80°C storage, 10 mL of 20 mg.mL^-1 methoxyamine in pyridine was added to the samples. The reaction was performed for 90 minutes at 28 °C under continuous shaking in the Eppendorf thermomixer. Then 50 mL of N-methyl-N-trimethylsilyl-trifluoroacetamide was added, and the reaction continued for 30 minutes at 38 °C. After cooling down, 45 mL of the derivatised sample was transferred to an Agilent vial for injection.

Metabolites were analysed by GC-MS 4 hours after derivatisation. One μL of the derivatised sample was injected in the split and splitless modes in Agilent 7890A gas chromatograph coupled to Agilent 5975C mass spectrometer, column Rxi-5SilMS from Restek (30 m with 10 m Integra-Guard column; Bellefonte, PA). The oven temperature ramp was 70 °C for 7 minutes, then 10 °C/ min to 330 °C for 4 minutes (run-length 36.5 minutes). Helium constant flow was 0.7 mL/min in splitless mode and 1 mL/min in split mode. Temperatures were as follows: injector ¼ 250 °C; transfer line ¼ 300°C; source ¼ 250 °C, and quadripole ¼ 150 °C. Samples were randomised. An alkane mix (C10, C12, C15, C19, C22, C28, C32, and C36) was injected in the middle of the queue for external retention index calibration. Five scans per second were acquired.

Raw Agilent data files were converted to NetCDF format and analysed with AMDIS 32 [99]. A home AMDIS retention index/mass spectra library, built from the NIST, Golm, and Fiehn databases and standard compounds, was used for metabolite identification. Peak areas were determined more precisely using the Target-Lynx software (Waters, France) after converting the NetCDF file to MassLynx format. Metabolites were normalised to the ribitol internal standard and to dry weight. A subset of metabolites was quantified absolutely using a response coefficient in both injection modes to ensure the absence of saturation.

Statistical analysis was done with Metaboanalyst 5.0, a web-based metabolomic processing software, using multivariate and hierarchical clustering statistical approaches [100]. If a metabolite was not identified in one of the plots but was presented in others, the missing values were replaced by 1·10^-15. In other cases, the remaining missing values were replaced by a limit of detection (1/5 of the minimum positive value of each variable). Before statistical analysis, all data were subjected to sample normalisation by median and auto-scaling. Principal component analysis (PCA) and partial least squares–discriminant analysis (PLS-DA) was done to determine the metabolomic variability within experimental plots. Metabolite change measurements were subjected to one-way ANOVA. Significant differences between individual means were determined using Fisher’s HSD (honestly significant difference) pairwise comparison test at the 5 % confidence level. To visualise the relationship among samples and among the top 50 metabolites, hierarchical clustering heatmaps with Euclidean distances as the similarity measure and Ward’s linkage as the clustering algorithm were conducted. Volcano plots were also generated for all detected metabolites, and the significant metabolites with p ≤ 0.05 and FC ≥ 2 were highlighted based on Student’s t-test. Significant metabolites were used for the pathway analysis, combining powerful pathway enrichment analysis results with pathway topology analysis to identify the most relevant pathways involved in the plant response to chronic radiation exposure. The Arabidopsis thaliana pathway library was used as a reference set.

4.7. Multi-omics data integration.

Multi-omics analysis was conducted with the web tool PaintOmics 4 [101] [101]. The closest to V. cracca model species in KEGG databases, Medicago truncatula L., was used as a proxy for multi-omics analysis. Galaxy public server (https://usegalaxy.eu/) [85] was used to perform blastp search of the differentially expressed V. cracca transcripts and proteins against the M. truncatula proteome (retrieved from Uniprot database, entry UP000002051). RStudio v.4.2.2 was used to prepare data for the multi-omics analysis. Background data on transcriptomics, proteomics, and metabolomics were provided for PaintOmics 4 software as log₂FC in irradiated samples compared to controls (pairs K x B, K x L, M x B, and M x L). From metabolomic data, log₂FC was provided only for those metabolites with KEGG ID. As relevant features for PaintOmics 4 tool, we used differentially expressed genes (DEGs) and differentially expressed proteins in both contaminated plots. As relevant metabolites, we chose those significantly different from the control for at least one experimental plot. After selecting and preparing datasets, the following parameters were used to run PaintOmics4: species: Medicago truncatula L.; databases: KEGG, MapMan, Reactome; clusters presented in data: generate automatically; metabolite class activity threshold: generate automatically.

5. Conclusions

Chronic exposure to ionising radiation induces complex and multidirectional changes in the Vicia cracca molecular phenotype. Plants exposed to chronic ionising radiation significantly altered their metabolism and gene expression patterns. At the level of carbon metabolism, the upregulation of glycolate, pyruvate, and fumarate biosynthesis pathways and the downregulation of citrate synthesis suggests that the cells are trying to bypass some steps of the TCA cycle to generate key metabolic intermediates. At the level of energy metabolism, the downregulation of glycolysis pathways, alterations in electron transport chain components, and changes in the photosynthetic machinery suggests a decrease in the cellular ability to generate ATP. The downregulation of fatty acid biosynthesis pathways and alterations in the glyoxylate cycle indicates a shift in the cellular metabolism away from fatty acid biosynthesis and towards using odd-chain fatty acids or amino acids as alternative energy sources. Profound changes in carbon metabolism may indicate the reallocation of carbon skeletons for secondary metabolite production to increase the antioxidant capacity of plants. Those effects were accompanied by suppression of ribosomal biogenesis and translation. Overall, the changes observed in carbon metabolism, energy metabolism, nitrogen metabolism, and gene expression suggest a complex biochemical response of Vicia cracca plants in the Chernobyl exclusion zone to chronic radiation exposure.

Moreover, chronically irradiated plants showed signs of redox imbalance and the activation of protective metabolic pathways. The upregulation of histones and few DNA repair machinery genes suggests that the plant relies on DNA packaging and homologous recombination to cope with chronic DNA damage. Yet, it is still unknown if the observed changes would be linked to strand break repair or DNA replication errors. The upregulation of peroxidases and genes involved in phenolic compound biosynthesis suggests increased ROS concentrations and a compensatory increase in antioxidant production. The nitrogen redirection towards the biosynthesis of secondary metabolites and profound changes in carbon metabolism indicate resource reallocation to cope with environmental conditions. The changes in ribosomal protein synthesis and RNA polymerases further suggest that the plant undergoes significant gene expression and translational pattern alterations.

Several molecules can be promising as radiation exposure biomarkers and for testing the resilience of Vicia plants to other stressors. Among them are chaperonin CPN60А, histones, including H2B and H4, transcription factor HY5, and chlorophyll-binding protein CAB1, as their genes show consistent upregulation for different plant species growing under chronic radiation exposure. Expression of the phenylalanine ammonia-lyase (PAL) gene and the abundance of its protein can be used to approximately assess the activation of secondary metabolites biosynthesis under chronic irradiation. Metabolites that seem responsive to chronic irradiation include the TCA components such as citrate, pyruvate, malonate, fumarate, and various inositol metabolites and RFOs. Especially interesting for future research is the involvement of sugar signalling in photosynthesis suppression under radiation exposure and the role of certain ribosomal proteins as a possible source for radiosensitivity manipulations. Future research could also explore the molecular mechanisms (probably epigenetic) underlying the observed changes in DNA repair machinery and histone expression.

In conclusion, the study of the transcriptomic responses of V. cracca to chronic radiation exposure highlights its ability to adapt to stressful environmental conditions through complex and multidirectional changes in gene expression and metabolism. Comparing responses of different plant species to chronic radiation exposure could reveal universal stress-tolerant mechanisms that could be manipulated to improve the resilience of crops. The findings of this study have implications for agricultural practices and ecosystem management strategies, as well as for the development of stress-tolerant crops.