Effects of cadmium stress on DNA methylation in soybean

Abstract

Epigenetics plays an important role in plant resistance to biotic and abiotic stresses, although few studies have focused on soybean epigenetics associated with cadmium stress. In this study, hydroponically grown soybean seedlings were treated with hydroponic nutrient solution containing various concentrations of CdCl₂ (0.5, 5, 15 mg L⁻¹) to investigate the effects of cadmium exposure on soybean genomic DNA methylation levels and patterns using methylation-sensitive amplified polymorphism (MSAP) analysis. Overall, the results revealed alterations of genomic DNA methylation levels and patterns in soybean seedlings under cadmium stress, with methylation level and methylation polymorphism rate both increasing in a dose-dependent manner with increasing cadmium concentration. The total methylation rates for genomic DNA of control (CK), S1 (0.5 mg L^−T CdCl₂), S2 (5 mg L⁻¹ CdCl₂) and S3 (15 mg L⁻¹ CdCl₂) groups were 27.92%, 29.78%, 33.2% and 35.79%, respectively, and methylation polymorphism rates were 14.46%, 16.57% and 18.15%, respectively. Thus, soybean cadmium stress-induced methylation changes mainly reflected increased genomic DNA methylation. After polymorphic methylated genomic DNA fragments were recovered from polyacrylamide gels, cloned and sequenced, sequence homology comparisons revealed that among 30 detected differentially methylated fragments, 15 sequences were highly homologous to genes encoding functional proteins already known to be associated with plant stress responses. These results indicate that soybean resists cadmium stress by increasing both the level of genomic DNA methylation and the genomic DNA remethylation rate. This work provides a methodological foundation for investigating molecular mechanisms underlying the soybean response to cadmium stress.

Keywords:

• DNA methylation
• methylation-sensitive amplified polymorphism
• homology analysis
• soybean
• Glycine max
• cadmium stresses

Introduction

Soybeans (Glycine max) are an important source of high-quality protein and vegetable oil, which both greatly benefit global human health and survival. However, abiotic stresses such as drought, high salinity, low temperature and exposure to heavy metals often lead to reduced soybean yields and crop losses [1]. Cadmium, a particularly harmful heavy metal and the most common heavy metal pollutant in farmland soil, has known detrimental effects on soybean growth, crop yield and quality. After cadmium is absorbed by plants, it accumulates in roots, stems, leaves and other tissues. Plants grown in soil with high cadmium levels exhibit root growth inhibition, plant dwarfism and chlorosis that directly compromise crop yield and quality [2].

DNA methylation, an important mechanism used by eukaryotic cells to regulate gene expression, maintain genome stability and sustain normal growth and development, has received considerable attention in recent years. In plants and many other eukaryotes, expression of genes within heterochromatin domains is determined in part by methylation of cytosines at position 5 of the pyrimidine ring (5-Me-C), with methylation of 20%–30% of cytosine residues in plant genomes typically observed under normal growth conditions [3]. Once DNA methylation of a given genomic cytosine residue occurs, it can be sustained for many generations as an epigenetic change that can trigger stress-induced gene evolution [4, 5]. Thus, stress induced by external environmental changes can lead to genomic DNA methylation site alterations that can alter the expression patterns of functionally adaptive genes. Therefore, the study of genomic DNA methylation changes that occur during the soybean response to cadmium exposure would enhance our understanding of key genes and molecular pathways involved in soybean adaptation to cadmium stress. This knowledge would be useful for guiding development of new soybean varieties with greater resistance to cadmium stress.

Methylation-sensitive amplified polymorphism (MSAP) analysis, based on amplified fragment length polymorphism (AFLP) analysis, is a common method used for generating and analyzing whole-genome methylation profiles. It has advantages of simple operation, high polymorphism-detection rate and effectiveness in detecting large numbers of cytosine methylation changes occurring at CCGG sites throughout an entire genome. In plants, this method has been effectively used to investigate methylation changes occurring in genomes of corn, rice, sorghum, pepper and other crops [6–9]. Specifically, MSAP technology has been used for detecting increases in genomic methylation level, upregulation of methyltransferase gene expression and nuclear chromatin structural changes in aquatic plants under cadmium stress. Such results have demonstrated that some plants can adapt to cadmium stress by regulating chromatin structure. In turn, chromatin structure can be controlled by altering genomic DNA methylation patterns within specific chromatin regions to control the balance between gene expression and gene silencing for groups of genes [10]. In this work, these methods were used to evaluate changes in genomic methylation profiles of soybean seedlings after exposure to cadmium stress, since few investigations with this focus have been reported to date.

In this study, in order to reveal epigenetics-based molecular mechanisms underlying the soybean response to cadmium stress, MSAP technology was used to analyze overall genomic DNA methylation levels and modes of methylation. Moreover, functional annotation analysis of DNA sequences corresponding to differentially methylated genomic DNA fragments obtained from plants under cadmium stress revealed the identities and putative functions of numerous genes with altered expression associated with cadmium stress. These results provide a methodological foundation to guide future studies of soybean cadmium stress response mechanisms. This information will be useful for guiding the development cadmium-resistant soybean varieties.

Materials and methods

Plant materials, enzymes and reagents

The soybean variety ‘Heinong 48’ was provided by Heilongjiang Academy of Agricultural Sciences. For 115 days from emergence to maturity, the active accumulated temperature should be gricultural Scienc2380 °C. It is an early maturing and high protein soybean variety. It has long leaves, purple flowers, gray hairs and light brown pods. The plant is 80–95 cm high, with 17 nodes of main stem and branches. The seeds are round, yellow and glossy, and the hilum is yellow. The 100 seeds weight is 22–25 g. From late April to early May, the yield per hectare is about 2480 480. From late April to early May, the yield per Table 1.

Table 1. Main reagents used in this study.

Reagent name	Manufacturer (Product model)
Tris base	CAISSON LABS, Beijing, China (CAS77-86-1)
EDTA Na2	BIOSHARP, Shanghai, China (CAS139-33-3)
Acrylamide	Beijing Hifi Technology Co., Ltd, Beijing, China (CAS79-06-1)
N’.N’-Methylene-bis-acrylamide	Aoran Fine Chemical Research Institute, Tianjin, China (CAS110-26-9)
Ammonium persulfate	Longxi Chemical Co., Ltd, Longxi, China (CAS7727-54-0)
Acetic acid glacial	Tianjin Fuyu Fine Chemical Co., Ltd, Tianjin, China (CAS64-19-7)
Silver nitrate	Shanghai No. 1 Reagent Factory, Shanghai, China (CAS7761-88-8)
Hpa II DNA Restriction Enzymes	TAKARA, Dalian, China (1053A)
Msp I DNA Restriction Enzymes	TAKARA, Dalian, China (1053A)
EcoR I DNA Restriction Enzymes	TAKARA, Dalian, China (1040A)
dNTP	TAKARA, Dalian, China (639132)
T4 DNA ligase	TAKARA, Dalian, China (2011A)
Taq DNA polymerase	TAKARA, Dalian, China (RR001A)
MiniBEST Plant DNA Extraction Kit	TAKARA, Dalian, China (9763)

Treatments and experimental grouping

Fully formed soybean seeds of uniform size were selected, soaked in distilled water for 30 min then soaked in 0.1% NaCl solution for 30 min. Next, seeds were sown in sterilized vermiculite in plastic pots that were placed in a greenhouse under controlled light and temperature (25 °C) conditions. After cotyledons fully emerged from seeds (with flattened appearance), plants were irrigated with 0.5× Hoagland nutrient solution every 2 days then were allowed to grow until the third compound leaf pair was visible. Next, soybean seedlings were selected based on similarity of growth vigor then selected seedlings were transferred to glass jars containing 1× Hoagland nutrient solution for cultivation.

Next, seedlings in jars were randomly assigned to the various treatment groups, with 3 jars used for each treatment (4 plants per jar). The nutrient solution was changed daily during treatment. After 24 h of culture in 1× Hoagland nutrient solution, seedlings were treated with or without CdCl₂ at various concentrations then plants were harvested after 8 days of treatment. Each treatment was as follows: Control group (CK): 1× Hoagland nutrient solution + sprayed distilled water; S1: 1× Hoagland nutrient solution + sprayed 0.5 mg L⁻¹ CdCl₂; S2: 1× Hoagland nutrient solution + sprayed 5 mg L⁻¹ CdCl₂; S3: 1× Hoagland nutrient solution + sprayed 15 mg L⁻¹ CdCl₂.

Extraction of soybean genomic DNA

After 8 days of growth, soybean seedlings were rinsed with deionized water then 0.1 g of soybean leaves was cut from the same node of each plant in each group (CK, S1, S2 and S3) then each leaf fragment was placed into a mg L^-1 centrifuge tube and frozen in liquid nitrogen. DNA was extracted using the SDS method [11, 12], then for each sample the DNA concentration was determined and the purity was assessed using a UV spectrophotometer and 1% agarose gel electrophoresis, respectively.

Methylation-sensitive amplified polymorphism (MSAP) analysis

The methylation-sensitive amplified polymorphism (MSAP) analysis procedure used here was based on the method reported by Pan Yajiao et al. [13], with modifications. Methylation-sensitive restriction enzymes HpaII and MspI were combined with restriction enzyme EcoRI (E + H, E + M) then genomic DNA for each experimental group was double-digested with each pair of enzymes. Sequences of adaptors and pre-amplification and selective amplification primers are shown in Table 2. The reaction system and reaction conditions are shown in Table 3. After a DNA denaturing step, PCR products were separated by size using 2% polyacrylamide gel electrophoresis, with detection of bands conducted using 0.1% silver nitrate staining and images recorded via digital imaging and electrophoretic bands recording.

Table 2. Adapter and primer sequences.

Adaptor and primer		Sequence (5’→3’)
Adaptors name	EcoRI adaptors (F)	CTCGTAGACTGCGTACC
	EcoRI adaptors (R)	AATTGGTACGCAGTCTAC
	H/M adaptors (F)	GATCATGAGTCCTGCT
	H/M adaptors (R)	CGAGCAGGACTCATGA
Preamplification primers	E0 (E + A)	GACTGCGTACCAATTCA
	HM0 (HM + 0)	ATCATGAGTCCTGCTCGG
Selecting amplification primers	E1 (E + AAC)	GACTGCGTACCAATTCAAC
	E2 (E + AAG)	GACTGCGTACCAATTCAAG
	E3 (E + ACA)	GACTGCGTACCAATTCACA
	E4 (E + ACT)	GACTGCGTACCAATTCACT
	E5 (E + ACG)	GACTGCGTACCAATTCACG
	E6 (E + AGC)	GACTGCGTACCAATTCAGC
	E7 (E + AGT)	GACTGCGTACCAATTCAGT
	E8 (E + AGA)	GACTGCGTACCAATTCAGA
	E9 (E + ATC)	GACTGCGTACCAATTCATC
	HM1 (HM + TAC)	ATCATGAGTCCTGCTCGGTAC
	HM2 (HM + TCT)	ATCATGAGTCCTGCTCGGTCT
	HM3 (HM + TCG)	ATCATGAGTCCTGCTCGGTCG
	HM4 (HM + TCC)	ATCATGAGTCCTGCTCGGTCC
	HM5 (HM + TGA)	ATCATGAGTCCTGCTCGGTGA
	HM6 (HM + TGT)	ATCATGAGTCCTGCTCGGTGT
	HM7 (HM + TGC)	ATCATGAGTCCTGCTCGGTGC
	HM8 (HM + TTA)	ATCATGAGTCCTGCTCGGTT
	HM9 (HM + TTC)	ATCATGAGTCCTGCTCGGTTC
	HM10 (HM + TTG)	ATCATGAGTCCTGCTCGGTTG

Notes: H/M: Hpa II and Msp I.

Table 3. Reaction system in the analysis of methylation sensitive amplified polymorphism (MSAP).

System	Reaction system	Reaction conditions
Restriction enzyme digestion	EcoR I 5U, Hpa II/Msp I 5U, 10 × Buffer 3 μL, DNA 400 ng. Added appropriate amount of ddH2O to a total capacity of 30 μl.	Digestion at 37 °C for 6 h, inactivation at 65 °C for 20 min and storage at 4 °C until use.
Ligation	Double-cut DNA product 20 µL, EcoR Iadaptors 5 pmol, H/M adaptors 50 pmol, T4 DNA ligase 1U, 10 × Buffer 4 µL. Added appropriate amount of ddH2O to a total capacity of 40 μL.	Ligation reactions proceeded at 16 °C for 13 h followed by inactivation at 65 °C for 20 min and storage at −20 °C for later use.
Pre-amplification	Ligated product 2 µL, E0 primer 5 pmol, H/M0 primer 5 pmol, dNTP (2.5 mmol L^−1) 1.6 µL, 10 × Buffer 2 µL, Taq polymerase 0.5U. Added appropriate amount of ddH2O to a total capacity of 20 μL.	Pre-denaturation at 94 °C for 5 min followed by 30 cycles of (denaturation at 94 °C for 30 s, annealing at 56 °C for 30 s, extension at 72 °C for 1 min) and a final extension step of 72 °C for 5 min. Completed reactions were stored at −20 °C for later use.
Selective amplification	Preamplification product diluted(10×) 3 µL, E primer 5 pmol, H/M primer 5 pmol, dNTP (2.5 mmol L^−1) 1.6 µL, 10 × Buffer 2 µL, Taq polymerase 0.5U. Added appropriate amount of ddH2O to a total capacity of 20 μL	Pre-denaturation 94 °C 5 min, denaturation 94 °C 30 s, annealing 65 °C 30 s, extension 72 °C 1 min, then a temperature increase of 0.7 °C per cycle for 13 cycles; followed by 23 cycles of (94 °C 30 s, 56 °C 30 s, 72 °C 1 min) and a final extension step at 72 °C for 10 min. Amplification reactions were stored at −20 °C.

Data collection

After conducting polyacrylamide gel electrophoresis and staining to reveal different reaction product bands, numbers of electrophoretic bands within each lane for each primer combination were counted. The presence of a band was recorded as a ‘1’ and absence of a band was recorded as a ‘0’ and results of triplicate experiments were analyzed to determine the statistical significance.

Analysis of DNA methylation level and pattern change

Sensitivities of homolyase HpaII and MspI to methylation status of DNA restriction sites differ such that MSAP analysis of double-digested restriction fragments can reveal the methylation status of an individual CCGG site, as well as the overall degree of restriction site methylation throughout an entire genome. DNA samples were digested with EcoRI/HpaII(H) and EcoRI/MspI (M), resulting in generation of four types of methylation band patterns (Types I-IV). For a Type I band pattern, bands were found in both H and M gel lanes, indicating that there was no methylation at the CCGG-site. For a type II band pattern, a band was found in the H-lane but not in the M-lane, indicating that single strand lateral methylation (hemimethylation/CHG methylation) occurred at the CCGG-site. For a type III band pattern, no band was found in the H-lane, but a band was found in the M-lane, indicating that the CCGG-site had undergone double-stranded medial methylation (full methylation/CG methylation). For a type IV band pattern, no band was found in the H-lane or M-lane, indicating that the CCGG-site had undergone double-stranded lateral methylation. Ultimately, only 3 types band patterns (type I, II and III) were detected on polyacrylamide gels in this study.

Results and discussion

Effects of cadmium stress on DNA methylation levels

For whole-genome MSAP analysis, 42 primer pair combinations were used for analyzing the soybean DNA methylation patterns in CK (control) S1 (0.5 mg L⁻¹), S2 (5 mg L^−,) and S3 (15 mg L⁻¹) samples, with relative frequencies of amplified bands analyzed statistically between groups. It can be seen from Table 4 that under Cd stress, MSAP absent target ratios (that indicate highly methylated sites) of soybean DNA of seedlings treated with S1, S2 and S3 (29.78%, 33.2%, 35.79%, respectively) were higher than the corresponding CK group ratio (27.92%). Total methylation (double-chain methylation) rates were 22.37%, 23.9%, 25.2% and 25.73% and half methylation rates were 5.55%, 5.89%, 8% and 10.06%, respectively. Thus, Cd stress exposure of groups S1, S2 and S3 appeared to lead to increased genomic DNA cytosine methylation, causing the rate of full methylation to exceed the half methylation rate. Thus, from these results it can be inferred that the predominant methylation status of CCGG sites in the genomes of the soybean plants under Cd stress was that of full methylation (double-stranded). Moreover, as compared with the control group, the total methylation level of soybean increased after Cd exposure such that as Cd concentration increased, the total DNA methylation level in soybean leaves gradually increased, indicating the existence of a significant dose-dependent relationship between the increase in total methylation level and Cd concentration.

Table 4. Effects of cadmium stress on soybean leaf genomic DNA methylation level.

MSAP amplified bands type	CK	S1	S2	S3
I type	1740	1613	1527	1455
II type	134	135	183	228
III type	540	549	576	583
Total bands	2414	2297	2286	2266
Total bands with methylation	674	684	759	811
Ratio of hemi-methylation (%)	5.55	5.89	8	10.06
Ratio of global methylation (%)	22.37	23.9	25.2	25.73
Ratio of total methylation (%)	27.92	29.78	33.2	35.79

Notes: Type I:Bands were found in both H and M lanes (denoted as ‘1,1’), indicating that CCGG site was not methylated; Type II:A band in H lane and no band in M lane(denoted as ‘1,0’), which indicated that single strand lateral methylation (hemimethylation/CHG methylation) occurred in CCGG site; Type III:No band in H lane and band in M lane(denoted as ‘0,1’), which indicated that CCGG site had double strand medial methylation (full methylation/CG methylation). Total bands = I+II + III; Total bands with methylation = II + III; Ratio of hemi-methylation = II/total bands; Ratio of global methylation = III/total bands; Ratio of methylation polymorphism = Total methylated bands/total bands.

Plants are known to adapt to external environmental stresses by employing mechanisms that involve epigenetic modifications (such as DNA methylation). DNA methylation is thought to be a relatively conserved form of epigenetic modification that plays important roles in maintaining normal plant growth and development and in plant adaptation to adverse external environmental conditions, although a change in methylation level usually has a negative impact on plants [14]. It has been suggested that the transcriptional initiation interaction between promoter and RNA polymerase is interrupted if the cytosine in the promoter region of the coding gene is methylated, resulting in gene expression silencing, as evidenced by observations that demethylation of gene promoters can activate genes and promote gene expression [15]. Thus, investigating changes in plant genome DNA methylation levels and patterns would greatly enhance our understanding of molecular and metabolic mechanisms underlying plant stress resistance, growth and development. In turn, such knowledge could be used to guide the development of soybean varieties with greater resistance to cadmium stress.

When plants are under heavy metal stress, they can regulate gene expression by altering the methylation of genomic DNA sites so as to silence or activate specific genes in order to cope with the stress [16]. In fact, studies have demonstrated that heavy metal stress is directly related to altered plant genomic DNA methylation levels and methylation patterns [17]. For example, after exposures to cadmium and lead stresses, the total methylation levels of rice and wheat genomes increased [17, 18], while the genomic DNA methylation levels in rapeseed seedlings after cadmium and chromium stress exposures exceeded the corresponding control levels and were positively correlated with the metal concentrations [19]. Conversely, during heavy metal stress due to exposures to heavy metals Ni, Cd and Cr, the changes in the overall cytosine methylation levels of other plant species, such as clover and hemp, decreased [20]. Here our results showed that the total methylated genomic DNA content of soybean seedlings after exposure to different Cd concentrations increased on average as compared with the control, with the total methylation level increasing in a dose-dependent manner with increasing Cd concentration. More specifically, although the total methylation level of soybean genomic DNA treated with 0.5 mg L⁻ Cd (29.78%) was barely greater than the control level (27.92%), the total genomic DNA methylation levels of soybean plants treated with 5 mg L^−C and 15 mg L^−a Cd were significantly greater than the control levels, as observed in rapeseed [19]. It has been suggested that plants may maintain genome stability through methylation to avoid oxidative damage caused by heavy metal stress [21]. In this study, the rates of total methylation and semi methylation of soybean genomic DNA in the control group were 22.37% and 5.55%, respectively, while the corresponding rates for plants exposed to Cd at a concentration of 0.5 mg L^−. were 23.9% and 27.73%, respectively, and to Cd at a concentration of 15 mg L⁻² were 5.89% and 10.06%, respectively. These results thus indicated that altered methylation of soybean genomic DNA in response to plant exposure to Cd was mainly due to full methylation of CCGG sites. Nevertheless, the change in the semi methylation level after Cd exposure exceeded that of full methylation, indicating that under Cd stress the soybean genome tended to undergo unstable hemimethylation as an important protective mechanism to cope with stress.

Effects of cadmium stress on DNA methylation patterns

By comparing the MSAP results for the soybean plants under Cd stress and the results for the control group, 15 band pattern differences were detected (Tables 5 and 6, Figure 1). It can be seen from Tables 5 and 6 that the numbers of type A (methylation type) loci in the genomic DNA of S1, S2 and S3 Cd-exposed plant groups were 182, 232 and 273, respectively, accounting for 7.74%, 9.9% and 11.69% of the total methylated polymorphic loci, respectively. Type B (demethylation type) loci numbers of 150, 149 and 143 were detected, respectively, accounting for 6.38%, 6.35% and 6.12% of total methylated polymorphic loci, respectively. Type C (indefinite type) loci numbers accounted for 0.34% of the total methylated polymorphic loci. Type D (with unchanged methylation status) polymorphic rate accounted for 85.54%, 83.43% and 81.85% of the total methylation polymorphic loci, respectively. As compared with CK (control), the methylation polymorphism rates for S1, S2 and S3 Cd-treated groups were 14.46%, 16.57% and 18.15%, respectively. It can be seen that with increasing Cd concentration, the methylation polymorphism rates of soybean genomic DNA increased, while the proportion of demethylated sites remained almost unchanged. Meanwhile, the number of methylated polymorphic sites gradually increased and ultimately exceeded the numbers of demethylated polymorphic sites. Therefore, the ratio of methylation to demethylation (net methylation) of soybean genomic DNA increased with increasing Cd concentration, indicating that high-level Cd exposure led to an increase in total methylation level of soybean genomic DNA.

Figure 1. Soybean leaf DNA methylation banding patterns. H1–H4: EcoR I+Hpa II combination; M1–M4: EcoR I+Msp I combination; 1: Control group leaf; 2: S1 group leaf; 3: S3 group leaf; 4: S4 group leaf; D: Monomorphism, D1: no methylation occurred, D2: Hemimethylation, D3: total methylation; A: Methylation, A1, A2 and A3: Remethylation, A4 and A5:Hypermethylation; B: Demethylation, B1, B2, B3, B4 and B5: Demethylation ; C: Uncertain type, C1 and C2: Uncertain type.

Table 5. Changes of methylation pattern of soybean genome under cadmium stress.

Type					Methylation status		Sites
Status type	Control		Stress		Control	Stress	S1	S2	S3
A1	1	1	0	1	CCGG GGCC	CCGG GGCC	71	93	94
A2	1	1	1	0	CCGG GGCC	CCGG CCGG GGCC GGCC	54	78	109
A3	1	1	0	0	CCGG GGCC	CCGG GGCC	23	24	32
A4	1	0	0	0	CCGG CCGG GGCC GGCC	CCGG GGCC	12	15	21
A5	0	1	0	0	CCGG GGCC	CCGG GGCC	22	22	17
B1	0	1	1	1	CCGG GGCC	CCGG GGCC	55	46	50
B2	0	0	1	0	CCGG GGCC	CCGG CCGG GGCC GGCC	18	25	34
B3	1	0	1	1	CCGG CCGG GGCC GGCC	CCGG GGCC	28	30	20
B4	0	0	0	1	CCGG GGCC	CCGG GGCC	28	38	29
B5	0	0	1	1	CCGG GGCC	CCGG GGCC	21	10	10
C1	0	1	1	0	CCGG GGCC	CCGG CCGG GGCC GGCC	6	6	5
C2	1	0	0	1	CCGG CCGG GGCC GGCC	CCGG GGCC	2	2	3
D1	1	1	1	1	CCGG GGCC	CCGG GGCC	1507	1441	1375
D2	1	0	1	0	CCGG CCGG GGCC GGCC	CCGG CCGG GGCC GGCC	57	74	80
D3	0	1	0	1	CCGG GGCC	CCGG GGCC	448	443	457

Notes: The underline ‘_’ indicated methylated bases. H: the lane of EcoRI + Hpa II enzyme combinations; M: the lane of EcoR I + Msp I enzyme combinations; C and CC represented methylated cytosine; 1: band present; 0: band absent.

Table 6. Effects of cadmium stress on soybean leaf genomic DNA methylation.

	S1	S2	S3
Total amplified bands (A + B + C + D)	2352	2347	2336
Total polymorphic bands (A + B + C)	340	389	424
Polymorphic bands ratio of polymorphic bands (%)	14.46	16.57	18.15
Type A band ratio of polymorphic bands (%)	7.74	9.9	11.69
Type B band ratio of polymorphic bands (%)	6.38	6.35	6.12
Type C band ratio of polymorphic bands (%)	0.34	0.34	0.34
Type D band ratio of polymorphic bands (%)	85.54	83.43	81.85

Plant regulation of gene expression mainly employs two methylation modes, remethylation and demethylation, with a change in the methylation mode in vivo involved in the protection of plants from stress. In this study, MSAP was conducted using 42 pairs of primers to amplify soybean genomic DNA, with results showing increasing soybean genomic DNA methylation polymorphism rates with increasing Cd concentration (0.5, 5 and 15 mg L⁻ⁱ), with methylation ratios of 7.74%, 9.9% and 11.69% and demethylation ratios of 6.38%, 6.35% and 6.12% observed, respectively. Thus, these results demonstrated that the main DNA methylation mode change in the soybean genome after Cd exposure mainly involved DNA remethylation.

In higher plants, DNA methylation pattern changes vary according to the type of heavy metal stress; for example, specific genomic DNA sequences in root cells of clover and marijuana exposed to heavy metals Ni, CD or Cr tend to undergo demethylation [20]. Conversely, methylation of some plant gene promoters and coding region sequences can silence the expression of associated genes during normal plant growth and development [22]. Usually, demethylation of plant DNA is thought to activate the expression of genes that induce synthesis of protective enzymes to resist harmful environmental conditions, since methylation cannot activate gene expression but can inhibit the transcription and expression of some genes. Nevertheless, after plants are exposed to certain stressors, their genomes can undergo remethylation as a defensive mechanism, mainly because methylation can inhibit the transcription of transposons and exogenous DNA and thus reduce the risk of damage resulting from genomic non-allelic translocation and recombination [23]. For example, Luo Jun et al. [24] showed that when different varieties of maize were subjected to Pb stress, DNA methylation patterns of resistant maize materials tended to be biased toward remethylation, while maize materials with lower resistance tended to be biased toward demethylation. As another example, Yang Jinlan et al. [25] found that under Cd stress, radish plants mainly initiated a dynamic Cd stress response that was biased toward demethylation to counter toxic Cd effects. In this study, absent target ratios of DNA methylation polymorphism after soybean seedling exposures to 0.5, 5, and 15 mg L⁻¹ Cd were 14.46%, 16.57% and 18.15%, respectively. Thus, the DNA methylation mode changes were mainly due to DNA methylation that may have shut down the expression and inhibited the transcription of some related genes through methylation, thus enhancing the resistance of soybean plants to toxic Cd exposures.

Sequencing and homology analysis of methylated DNA fragments

To identify the functions of the proteins encoded by differentially methylated genomic fragments sequenced in this work, the NCBI Blast search tool (http://www.ncbi.nlm.nih.gov/) was used for conducting Blast homology searches of DNA sequences obtained from 30 differentially methylated bands. The results (Table 7) showed that 9 of 30 sequences encoded proteins that were homologous with proteins in the database with unknown functions (E3, E7, E11, E12, E16, E18, E22, E24, E28). Meanwhile, six fragments encoded proteins with no homologues in the database (E6, E13, E15, E21, E23, e25), while the remaining 15 sequences had high homology with genes encoding proteins with known functions as follows: CCR3 serine/threonine-protein kinase (E1), wnk8 serine/threonine-protein kinase (E26), ZF HD homologous protein (E4), envelope protein (E8), atscl13 protein (member of GRAS family) (E27), malic acid transporter (E29), glycosamine-rich RNA binding protein (E2), antisense mediated mRNA degradation regulator (Upf1) (E5), GDP-l-fucosynthase (E10), protein translation initiation factor Sui1 (E14, E30), yipf1 (E19, E20) homologous protein, galactosaldehyde (E9) and microtubulin (E17). Notably, analysis of the functions of these proteins indicated they were mainly associated with cell metabolic pathways and plant stress resistance-associated pathways.

Table 7. Blast result of DNA differential methylation bands.

Fragment number	Fragment length	Variation type	GenBank Login ID	Blast
E1	213	A4	gb|KHN33070.1|	Putative serine/threonine-protein kinase-like protein CCR3 [Glycine soja]
E2	183	A4	ref|XP_003555886.1|	PREDICTED: glycine-rich RNA-binding protein 3, mitochondrial-like isoform X1 [Glycine max]
E3	131	B4	gb|ACU22922.1|	unknown [Glycine max]
E4	267	B3	gb|KHN13637.1|	ZF-HD homeobox protein [Glycine soja]
E5	193	B3	ref|XP_003517385.1|	PREDICTED: regulator of nonsense transcripts 1 homolog [Glycine max]
E7	266	A3	gb|ACU20615.1|	unknown [Glycine max]
E8	285	A3	gb|AAO73522.1|	Envelope-like protein [Glycine max]
E9	368	B4	ref|XP_003551855.1|	PREDICTED:galactolipase DONGLE, chloroplastic-like [Glycine max]
E10	203	B3	ref|XP_003552287.1	PREDICTED: putative GDP-L-fucose synthase 2-like [Glycine max]
E11	166	B4	gb|ACU19162.1|	unknown [Glycine max]
E12	329	A3	emb|CEM21753.1|	unnamed protein product
E14	240	B4	ref|XP_006586085.1|	PREDICTED: protein translation factor SUI1 homolog 2 [Glycine max]
E16	130	B5	ref|XP_003555498.1|	PREDICTED: uncharacterized protein LOC100812185 [Glycine max]
E17	122	B5	gb|ACU22922.1|	unknown [Glycine max]
E18	337	A3	gb|KHN00766.1|	DNA repair helicase UVH6 [Glycine soja]
E19	285	B1	ref|XP_003519363.1	PREDICTED: protein YIPF1 homolog isoform X1 [Glycine max]
E20	285	B1	ref|XP_003519363.1|	PREDICTED: protein YIPF1 homolog isoform X1 [Glycine max]
E22	235	B1	ref|XP_006595596.1|	PREDICTED: uncharacterized protein LOC100804748 isoform X4 [Glycine max]
E24	131	B5	gb|ACU22922.1|	unknown [Glycine max]
E26	293	A1	ref|XP_003542159.1	PREDICTED: serine/threonine-protein kinase WNK8 [Glycine max]
E27	145	B1	ref|XP_003534887.1|	Scarecrow-like protein 13-like [Glycine max]
E28	375	B3	ref|XP_006577601.1|	uncharacterized protein LOC102670534 [Glycine max]
E29	469	B1	ref|XP_003555230.1|	Aluminum-activated malate transporter 4-like isoform X1 [Glycine max]
E30	191	B4	ref|XP_006586085.1|	PREDICTED: protein translation factor SUI1 homolog 2 [Glycine max]

Functional cluster analysis revealed that proteins encoded by 15 differentially methylated genomic DNA fragments had homology to proteins with known functions that were associated with plant responses to salt, drought, low temperature and heavy metal stresses. For example, the UPF1 protein is an important regulator of nonsense-mediated mRNA decay (NMD). NMD selectively degrades mRNA that contains premature translation termination codons (PTCs) to minimize cell damage caused by accumulation of C-terminally truncated proteins and is a key mechanism for controlling post-transcriptional mRNA quality in eukaryotic cells [26]. PTCs, which mainly result from DNA- or RNA-level base changes that introduce translational stop codons into mRNA, are detected by UPF1 factor, the substrate of NMD [27]. A combination of alternative splicing (AS) and nonsense-mediated RNA degradation (AS/NMD) is employed as a post-transcriptional monitoring mechanism for regulating gene expression. Studies have shown that AS/NMD can induce the expression and transcription of many genes, including transcription factors and stress response genes that participate in normal plant growth, development and adaptation to stresses, [28], with mRNA degradation playing an essential role in maintaining plant health and vigor. In this work we found that demethylation of a gene homologous to UPF1 (E5) may play a role in the transcriptional regulation of mRNA degradation factor induced by plant exposure to cadmium as an adaptive mechanism employed by soybean plants to resist cadmium stress [29].

Conclusions

The total methylation and methylation polymorphism ratios of soybean genomic DNA increased after cadmium exposure in a dose-dependent manner with increasing Cd concentration triggering a change in methylation mode that increased genomic DNA methylation overall. Of the 30 differentially methylated DNA fragments characterized here, 15 had sequences that were highly homologous to genes encoding proteins with known functions that were mainly associated with plant stress responses. Strategies targeting such genes may be useful for developing soybean varieties with greater resistance to cadmium stress.

Ethics statement

There are no researches conducted on animals or humans.

Acknowledgements

We hereby thank the Heilongjiang Academy of Agricultural Sciences for providing the soybean material for breeding.

Data availability statement

The data that support the findings of this study are available from the authors upon reasonable request.

Disclosure statement

There are no conflicts to declare.

Funding

This study was supported by the Harbin Science and Technology Innovation Talent Research Project (RC2013QN002103), the Heilongjiang Provincial Postdoctoral Funding Project (LBH-Z13185), the Heilongjiang Provincial Department of Education Science and Technology Research Project (12511154) and the Harbin Normal University Doctoral Research Fund (KGB200903), Harbin Normal University graduate innovative scientific research project “Cloning and functional study of soybean squalene coxidase GmSQE3 gene”, Postgraduate research project of Harbin Normal University (project number: HSDSSCX2021-03; grant ID 10,000RMB).

Funding

The author(s) reported there is no funding associated with the work featured in this article.