Metabolomic analysis of the effects of cadmium and copper treatment in Oryza sativa L. using untargeted liquid chromatography coupled to high resolution mass spectrometry and all-ion fragmentation^†

Abstract

While the knowledge of plant metabolomes has increased in the last few years, their response to the presence of toxicants is still poorly understood. Here, we analyse the metabolomic changes in Japanese rice (Oryza sativa var. Japonica) upon exposure to heavy metals (Cd(ii) and Cu(ii)) in concentrations from 10 to 1000 μM. After harvesting, rice metabolites were extracted from aerial parts of the plants and analysed by HPLC (HILIC TSK gel amide-80 column) coupled to a mass spectrometer quadrupole-Orbitrap (Q-Exactive). Full scan and all ion fragmentation (AIF) mass spectrometry modes were used during the analysis. The proposed untargeted metabolomics data analysis strategy is based on the application of the multivariate curve resolution alternating least squares (MCR-ALS) method for feature detection, allowing the simultaneous resolution of pure chromatographic profiles and mass spectra of all metabolites present in the analysed rice extracts. All-ion fragmentation data were used to confirm the identification of MCR-ALS resolved metabolites. A total of 112 metabolites were detected, and 97 of them were subsequently identified and confirmed. Pathway analysis of the observed metabolic changes suggested an underlying similarity of the responses of the plant to Cd(ii) and Cu(ii), although the former treatment appeared to be the more severe of the two. In both cases, secondary metabolism and amino acid-, purine-, carbon- and glycerolipid-metabolism pathways were affected, in a pattern consistent with reduction in plant growth and/or photosynthetic capacity and with induction of defence mechanisms to reduce cell damage.

Graphical Abstract

Evaluation of the effect of metal exposure (cadmium/copper) on metabolic pathways of rice.

Significance to metallomics

The assessment of the effects caused in crop plants by exposure to pollutants (i.e. metals) is of key importance. In this work, metabolic changes caused by the presence of different concentrations of cadmium and copper in irrigation water were evaluated taking Japanese rice (Oryza sativa L.) as a model organism. From the results, the affected pathways could be determined to allow a biological interpretation regarding the effects on plant growth and defense.

Introduction

Heavy metals are important constituents of the Earth’s crust and geological processes. Some natural processes, for instance erosion of the underground geological material, and emissions from volcanoes or forest fires, contribute to their presence in the environment. Nevertheless, during the last century, some anthropogenic activities such as mining, industry or agriculture, have altered heavy metal distribution on the earth surface, increasing significantly the levels of these pollutants in the environment. Due to the scarcity of arable land around the planet, especially in industrialized countries, contaminated soils must be used by farmers and their heavy metal pollution may cause several environmental problems and risks to human health, including contamination of edible plants.^1–3 Plants are affected by heavy metal pollutants because they are easily absorbed by roots and readily translocated to the aerial parts of the plants.¹ Plant responses to heavy metal stress include numerous defense mechanisms and stress-inducible reactions, such as synthesis or metal-binding peptides, antioxidative mechanisms produced by the generation of reactive oxygen species and formation of highly active signaling compounds.^4,5

Biomarker discovery appears as a powerful tool to study the effects of heavy metal pollution in plants. A biomarker can be defined as a substance used as an indicator of a biological state that can cause alterations in organisms.⁶ The identification of these biomarkers is useful to study pathogenic processes, risk and progression of diseases and risk assessment and toxicity of pollutants.^7,8 When biomarkers are specific to a given pollutant, they can be used as a monitoring index for long-term monitoring studies.^6–8 It is clear that the need for sensitive and specific biomarkers is increasing and much effort is being made in the research related to its discovery. Omics sciences (genomics, proteomics, transcriptomics and metabolomics) can provide an untargeted, global knowledge of organisms' physiology and of their responses to environmental inputs and, consequently, they are really useful in biomarker research. In particular, the use of metabolomics in biomarker discovery is based on the assumption that biotic and abiotic stress causes disruptions of biochemical pathways leading to a metabolic fingerprint characteristic of each stressor or group of stressors.^7,9 There are two principal metabolomics approaches: targeted and untargeted. The first one is only focused on analysing selected molecular classes, whereas untargeted metabolomics aims to screen the entire metabolite content of biological samples. In this work, the untargeted approach is used because it enables the simultaneous profiling of the largest number of metabolites present in the sample. Also, untargeted metabolomics provides the possibility of finding which metabolites show changes in their concentration under a treatment and elucidate previously unexplored biological pathways.⁷

Considering that the samples analysed in plant metabolomics can be diverse and complex, the analytical techniques used in this field should have high separation power. The most frequently utilized techniques have been nuclear magnetic resonance (NMR) and separation techniques (gas chromatography (GC), liquid chromatography (LC) and capillary electrophoresis (CE)) coupled to mass spectrometry (MS). NMR is less sensitive than MS-based techniques, and its ability to detect low abundance metabolites is limited. For this reason, MS-based techniques are nowadays generally the chosen option.^10–12 In untargeted metabolomics studies, the high complexity of the samples and the presence of a wide variety of metabolites, many of them at low concentrations and with very different physicochemical properties, make metabolite identification and determination really challenging.^6,13 MS-based techniques allow a putative identification of the metabolites by comparing the m/z of the molecular ion with the theoretical metabolite exact mass. To ensure complete identification of the metabolites, currently high-resolution mass spectrometry (HRMS) and tandem mass spectrometry (MS/MS), as well as comparison with standards are used. In this work, the selected methodology is LC-HRMS based on hydrophilic interaction liquid chromatography (HILIC). This strategy allows the analysis of a broad range of compounds without the need for chemical derivatization.¹³ For the analysis of untargeted compounds HRMS and all-ion fragmentation (AIF) without selecting precursor ions provide accurate mass and structural information of the unknown metabolites.

Plant metabolomics is growing as an essential System Biology Tool in plant science and, in particular, for crop-enhancing projects. Plants produce extremely complex metabolomes, with large numbers of metabolites presenting a huge variety of structures and relative abundances. These metabolites play important roles in plant growth, development, and response to environment changes. At the same time, metabolome composition constitutes the chemical base of crop yield and quality, determining the nutritional properties both for humans and livestock.¹⁴ Rice (Oryza sativa L.) is a plant of remarkable alimentary and economic importance, being one of the cereals most consumed by the world population.^15,16 The cultivar Japonica Nipponbare was selected for this work, as it is one of the best-known rice varieties, with a relatively short period of growth and easy genetic modification.^5,11

Untargeted metabolomics, particularly from plants, produces extremely complex datasets, the processing of which constitutes a crucial step of the whole analytical process. Multivariate data analysis tools allow for a reliable evaluation of metabolite concentration changes. In a previous work, the authors demonstrated the usefulness of two different untargeted metabolomics data analysis approaches for the study of Japanese rice under heavy metal stress.¹⁷ The main objective of this previous work was the preliminary evaluation and validation of the proposed data analysis tool. In contrast, the goal of the present work is focused on the analytical identification and confirmation of the metabolites of Japanese rice showing altered concentrations due to cadmium and copper treatments. For this purpose, an untargeted LC-HRMS metabolomic approach has been used. Multivariate curve resolution alternating least squares (MCR-ALS) has been employed to resolve and detect the most important metabolites. Then, all-ion fragmentation (AIF) has been applied to confirm metabolite identifications, allowing a biological interpretation of the main pathways affected.

Experimental

Reagents

Cadmium chloride hydrate (≥98.0%), copper(ii) sulphate pentahydrate (≥98.0%), ammonium acetate (≥98.0%) and LC-MS grade water, acetonitrile (≥99.9%), methanol (≥99.9%) and acetic acid were supplied by Sigma-Aldrich (Steinheim, Germany). Chloroform was obtained from Carlo Erba (Peypin, France). Piperazine-N,N′-bis(2-ethanesulfonic acid) (PIPES) (≥99.0%) was used as an internal standard (Sigma-Aldrich, Steinheim, Germany).

1000 μM stock solutions of cadmium (Cd(ii)) and copper (Cu(ii)) were prepared weekly by the dissolution of appropriate amounts of cadmium chloride hydrate and copper(ii) sulphate salts. Solutions containing 10, 50 and 100 μM metals were prepared weekly by diluting the 1000 μM stock solutions. All these solutions were stored at 6 °C until their use.

Water used for plant watering, for preparing cadmium and copper solutions, and during the extraction procedure was purified using an Elix 3 coupled to a Milli-Q system (Millipore, Belford, MA, USA), and filtered through a 0.22 μm nylon filter integrated into the Milli-Q system.

Plant growth, stress treatment and metabolite extraction

Plant growth, stress treatment and metabolite extraction were performed using a procedure previously described.¹⁷ First, Oryza sativa var. Japonica Nipponbare seeds, obtained from the Centre for Research in Agricultural Genomics (CRAG), were incubated for two days at 30 °C in a wet environment. After this period, plants were grown on an Environmental Test Chamber MLE-352H (Panasonic®) for 22 days under white fluorescent light. Temperature, relative humidity, and light long-day conditions at the chamber were set as described in Fig. S1 in the ESI†.

During the first ten days of growth, rice plants were watered with Milli-Q water three times per week. After that, the treated plant samples were subjected to irrigation water containing different concentrations of Cd(ii) or Cu(ii), whereas the control plant samples were watered only with Milli-Q water until harvest. Metal concentrations were 10, 50, 100 and 1000 μM. After harvest, the aerial parts of the rice samples were frozen at liquid nitrogen temperature in order to quench metabolism. Samples were stored at −80 °C until extraction.

Before extraction, the aerial parts of the rice samples were ground to a fine powder using a liquid nitrogen mortar and lyophilized for 24 h until dryness. Metabolite extraction was carried out by dispersing 40 mg of the dried tissue in 1 mL of MeOH in a 2.0 mL Eppendorf tube. Then, the mixture was vortexed for 1 min and sonicated for 10 min; this step was repeated twice. After centrifuging for 20 min at 14 100 × g, a 750 μL aliquot of the supernatant was transferred to a 1.5 mL Eppendorf tube. Then, 500 μL of chloroform and 400 μL of water were added. After that, the mixture was vortexed for 1 min, incubated for 15 min at −4 °C, and centrifuged for 20 min at 14 100 × g. Finally, the aqueous fraction was transferred to a 1.5 mL Eppendorf tube, evaporated to dryness under nitrogen gas, and reconstituted with 450 μL of acetonitrile/water (1 : 1 v/v). For internal standard quantification, 50 μL of 50 mg L⁻¹ solution of the internal standard (PIPES) were added to the extract. All of the extracts were stored at −80 °C until they were analyzed. Before injection, the samples were filtered through 0.2 μm nylon filters (Pall Life Sciences, Port Washington, NY, USA).

LC-MS analysis

Chromatographic separation was carried on an Accela UHPLC system (Thermo Scientific, Hemel Hempstead, UK) using the method described elsewhere.¹⁷ The LC column was an HILIC TSK gel amide-80 column (250 × 2.0 mm i.d., 5 μm) with a guard column (10 × 2.0 mm i.d., 5 μm) of the same material provided by Tosoh Bioscience (Tokyo, Japan). Since the main aim was to separate metabolites, which are highly polar molecules, the use of HILIC columns was recommended for their analysis. An elution gradient was produced using solvent A (acetonitrile) and solvent B (acetic acid:ammonium acetate buffer 3 mM at pH 5.5) as follows: 0–3 min, isocratic gradient at 5% B; 3–27 min, linear gradient from 5 to 70% B; 27–30 min, isocratic gradient at 70% B; 30–32 min back to the initial conditions at 5% B; and from 32 to 40 min, at 5% B. The mobile phase flow rate was 0.15 mL min⁻¹ and the injection volume was 5 μL.

A Q-Exactive (Thermo Fisher Scientific, Hemel Hempstead, UK) equipped with a quadrupole-Orbitrap mass analyser was used as a mass spectrometer. The ionization source employed was a heated electrospray (HESI) in negative ion mode. Mass spectra were acquired in profile mode at a resolution of 70 000 FWHM (full width half maximum) at m/z 400. Working parameters were as follows: electrospray voltage, 3.0 kV; sheath gas flow rate, 25 arbitrary units (a.u.); auxiliary gas flow rate, 10 a.u.; heated capillary temperature, 300 °C; S-lens level, 60%; automatic gain control (AGC), 3 × 10⁶; and the maximum injection time was set at 200 ms with two microscans/scan. The full scan mass range was from m/z 90 to 1000. All ion fragmentation (AIF) was also performed with a normalized collision energy (NCE) of 35 eV.

Analysis of metal content in rice samples

After harvesting, the aerial parts and roots of the treated samples were lyophilized for 24 h until dryness and ground to a fine powder. After that, the plants were subjected to Teflon digestion. In the case of the aerial parts, 3 mL of HNO₃ and 1 mL of H₂O₂ were added to 100 mg of sample, and then the mixture was digested in a Teflon reactor for three days. The digestion of root samples is explained in the ESI.† Finally, the samples were diluted with 30 mL of water. Metal concentrations were determined using inductively coupled plasma-mass spectrometry (ICP-MS). An Agilent 7500ce model ICP-MS system (Agilent, Santa Clara, CA, USA) from the scientific and technological centers of the University of Barcelona (CCiTUB) was used for the analysis. The RF power was set at 1550 W, isotopes used during the analysis were Cu⁶³ and Cd¹¹¹ and the internal standard was Rh.

Data analysis strategy

Thermo Fisher raw chromatographic data files (raw format) were converted to the standard CDF format by the FileConverter function of Xcalibur™ 2.2.44 software (Thermo Scientific, Hemel Hempstead, UK). These data files were then imported into the MATLAB environment (release 2014b, The Mathworks Inc, Natick, MA, USA) by using the appropriate functions of the MATLAB Bioinformatics Toolbox (4.3.1.version) and in-house built routines. Each sample was represented by a data matrix containing the acquired retention times on 1020 rows (from 0 to 40 min) and the detected m/z values on 38 000 columns from 90 to 1000 m/z. Since in the Q-Exactive mass spectrometer, the mass spectra were acquired at high resolution in the Orbitrap mass analyser, the obtained matrices contained information from 90 to 1000 m/z with an accuracy of ±0.0001. The high amount of information obtained makes it difficult to process the dataset. To reduce the computer storage requirements and facilitate calculations, the total number of columns (i.e. m/z values) was reduced by using a binning approach (grouping mass values into a number of bins within a particular m/z range, in this case 0.05). Therefore, the final data matrix for each sample had 1020 rows (retention time from 0 to 40 min) and 18 200 columns (from m/z 90 to 1000 with a resolving power of 0.05 m/z units).¹⁸

The data analysis strategy is shown in Fig. 1. The first step was the identification of chromatographic regions able to discriminate between control and treated samples by using partial least squares-discriminant analysis (PLS-DA) to evaluate total ion current (TIC) chromatograms. Only the PLS-DA regions showing a high discrimination power between samples were then analysed by means of MCR-ALS, which resolved the pure elution profiles and mass spectra of the metabolites present in the analysed sample. The accurate mass of the MCR-ALS resolved mass spectra was verified by checking the HPLC-HRMS raw data. Then, metabolites were identified using their exact mass and confirmed by their AIF mass spectra. Finally, the statistical significance of the metabolite concentration changes between sample groups was assessed for the identified metabolites.

Fig. 1

MCR-ALS based untargeted metabolomics data analysis workflow.

Selection of relevant chromatographic regions

The first step in metabolite identification was the selection of relevant chromatographic regions where metabolite concentration changes can be expected. This selection was performed applying PLS-DA to the TIC chromatograms of the analysed rice samples. Prior to PLS-DA analysis, TIC chromatograms were normalized by dividing them by the peak area of the internal standard (PIPES) in the considered sample. Then, baseline correction (using a weighted least squares method)¹⁹ and peak alignment (using the correlation optimized warping (COW) method)²⁰ were performed. Finally, TIC chromatograms were only mean-centered (autoscaling was not used to avoid giving baseline noisy regions too much influence in the model). PLS-DA²¹ is a multivariate regression method oriented to discriminate among different groups of samples. The TIC chromatograms of every sample (X, predictor variables) were related with a vector describing the class membership (y, predicted variable).^22,23 In this case, for every metal class, the following were selected: control samples (class 0) and samples treated with Cd(ii) or Cu(ii) (class I). Apart from discriminating among different groups of samples, PLS-DA also provides information about which are the most important variables (in this case, chromatographic retention times) for discrimination. For instance, variable importance on projection (VIP) scores can be used for this purpose.²³ VIP scores are a weighted sum of squared PLS variables which measure the importance of each predictor variable to the final PLS model.²³ The “greater than one” rule is often used as a variable selection criterion, because the average of squared VIP scores is equal to 1. This means that variables with a VIP score greater than 1 can be considered relevant for sample discrimination.²⁴ In the particular case described in this work, control samples were distinguished from treated samples (10, 50, 100 and 1000 μM). Metabolites of interest are those that have different concentrations in control samples compared to treated samples. Therefore, VIP plots should spotlight chromatographic regions (retention times) that allow discriminating between samples because of the differences in the concentrations of their metabolites. PLS-DA results were assessed by leave-one-out cross-validation.²⁵ PLS Toolbox 7.8 working under MATLAB was used in all calculations.

Multivariate curve resolution by alternating least squares

Elution profiles and mass spectra of the metabolites present in every selected chromatographic region were resolved by the MCR-ALS method. MCR-ALS is a powerful chemometric tool used for the investigation and resolution of pure component contributions in unresolved mixtures. This method has been applied to a great variety of examples from different fields, such as hyphenated and multidimensional chromatographic systems, -omics data sets or spectroscopic images, among others.^15,26–29 In the case of untargeted LC-MS metabolomics studies, MCR-ALS can resolve a large number of overlapped chromatographic peaks, obtaining both the chromatogram and the mass spectrum of the components (ideally a single metabolite). It is especially helpful in untargeted metabolomics analysis, where there is no previous knowledge about the chemical compounds (metabolites) present in the analysed sample.^{12,17,18,30–32}

MCR-ALS decomposes experimental data sets arranged in a data matrix according to the following bilinear model:

$$D=C{S}^T+E$$

1

where D (size I × J) represents the experimental LC-MS data set (i.e. a single rice sample) organized in a data matrix in which its rows are the MS spectra at all retention times (i = 1, …, I), and columns are the chromatograms at all m/z values (j = 1, …, J). According to eqn (1), the D matrix is decomposed into the product of two-factor matrices, C and S^T, which are respectively the matrix of resolved elution profiles, C (size I × N), and the matrix of their corresponding mass spectra, S^T (size N × J). N represents the number of resolved components using the MCR-ALS method. The matrix E (size I × J) contains the residuals not explained by the model using the N considered components.^15,33 This data analysis strategy was extended to the simultaneous analysis of several samples. For instance, in the case of this work, different augmented data matrices (one for each metal treatment and selected chromatographic region), each with a total number of 40 data matrices (samples) were analysed. For instance, the augmented data matrix for the cadmium treated samples contained 8 control samples and 32 samples treated with Cd(ii) at four different concentration levels (10, 50, 100 and 1000 μM), and with 8 replicates for each level. Analogously, an augmented data matrix for copper treated samples contained 8 control samples and 32 samples treated with Cu(ii), at four different concentration levels (10, 50, 100 and 1000 μM). In both cases, these 40 samples were arranged in a single column-wise augmented data matrix (D_aug), which contains the individual data matrices (D_x where x = 1, …, 40), one for each rice sample, and settled one on the top of the other. This column-wise augmented data matrix (D_aug) can be decomposed by the MCR-ALS method using the same bilinear model:

$$D_{\text{aug}}=\left[\begin{array}{c}{D}_1\\ D_2\\ D_3\\ ⋮\\ D_{40}\end{array}\right]=\left[\begin{array}{c}{C}_1\\ C_2\\ C_3\\ ⋮\\ C_{40}\end{array}\right]S^T+\left[\begin{array}{c}{E}_1\\ E_2\\ E_3\\ ⋮\\ E_{40}\end{array}\right]=C_{\text{aug}}{S}^T+E_{\text{aug}}$$

2

According to eqn (2), resolved mass spectra (S^T) were forced to be the same for the common components (metabolites) in the different analysed rice samples (chromatographic runs). However, the elution profiles of the same component resolved in the column-wise augmented concentration matrix C_aug were allowed to be different for each one of the analysed rice samples (C_x (x = 1, …, 40)), as shown in eqn (2).^15,33 MCR-ALS solves eqn (1) (D, single data matrix case) or eqn (2) (D_aug, augmented data matrix case) starting with an estimation of the number of components of the considered matrix by the singular value decomposition (SVD) method.³⁴ This number is only an initial approximation, as the final number of components is decided by taking into account the data fitting results and the reliability of the resolved profiles. An initial estimate of S^T for the iterative optimization should be provided, which can be easily obtained by selection of the purest variables using a method based on the SIMPLISMA approach.^35,36 Then, the ALS constrained optimization procedure is carried out. The application of constraints provides chemical meaning to the pure mathematical solution. In the case of this work, the applied constraints were non-negativity (for chromatographic elution and spectra profiles of every component) and normalization of the spectra profile of every component (equal height).^33,37,38

Before applying MCR-ALS, individual data matrices of each rice sample (D_x) were normalized by dividing the individual values of each of them by the value of the chromatographic peak area of the internal standard (PIPES) of the considered sample. MCR-ALS analysis was carried out using the MCR-ALS toolbox freely available at www.mcrals.info.

Metabolite identification

After MCR-ALS analysis, matrices C_aug and S^T from MCR-ALS results were evaluated to detect the metabolites whose concentrations change due to the metal treatment of rice plants. The evaluation of the elution profiles resolved in C_aug for each component allowed the determination of the pure metabolites that showed a significant change in their concentrations between the different groups of samples (control and metal treated samples). For instance, Fig. 2A shows an example of the peak areas and the MS spectrum resolved for a component when control and Cd(ii) treated samples were considered. In this case, Cd(ii) produced a significant decrease in the peak areas of metabolites in comparison with the control samples. Therefore, from the corresponding resolved mass spectrum of this component (see Fig. 2B), it was possible to estimate the mass of the diagnostic ion associated to the metabolite causing the observed differences in the area of the resolved chromatogram when the control and the treated samples are compared. Since HESI is a soft ionization source, in most of the cases, the mass of the diagnostic ion corresponds to the mass of the deprotonated metabolite. However, it is now possible to deduce the accurate mass (four decimal figures) of the selected metabolite looking for the HRMS Q-Exactive Orbitrap raw mass spectrum of the considered metabolite. This allowed an initial identification (elemental composition) of the selected metabolite to be performed with a relative mass error lower than 5 ppm, which fulfils Directive 2002/657/CE mass spectrometric detection performance criteria and requirements.³⁹ This tentative identification of the metabolites changing concentration due to heavy metal stress was performed by comparing the accurate mass of the MCR-ALS resolved metabolites (after checked in the HRMS Q-Exactive Orbitrap raw data) with the theoretical exact mass values included in public databases such as MassBank,⁴⁰ Metlin,⁴¹ HMDB,⁴² and MetaCyc.⁴³

Fig. 2

An example of detection of pure metabolites showing a change in their concentration between control and Cd(ii) treated samples. (A) Area of the resolved elution profiles of a particular metabolite for each sample derived from c_aug,i. (B) Resolved mass spectrum of the particular metabolite (s^T_i in eqn (2)).

The last step in metabolite identification was the confirmation of the hypothesized chemical structure. This confirmation was performed by comparing the experimental AIF mass spectra with the corresponding theoretical MS/MS spectrum of the suspected metabolites (checked from MassBank⁴⁰ and the library of National Institute of Standards and Technology (NIST)⁴⁴ databases). In order to achieve the complete identification of the found metabolites, the requirements of the previously mentioned Directive 2002/657/CE³⁹ were also followed. To comply with this Directive, accurate mass measurements of diagnostic and product ions were required to have a relative error lower than 5 ppm. Furthermore, four identification points (IPs) were mandatory for positive identification. Since mass spectra were acquired with a resolution higher than 20 000 FWHM, the HRMS precursor ion and two products earn 2 and 2.5 IPs respectively, which made the achievement of the four mandatory IPs possible.

Candidate metabolites were searched on-line in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database,⁴⁵ specifically selecting O. sativa-annotated compounds when necessary.⁴⁶ The same strategy was followed to classify the identified metabolites into KEGG metabolic pathways.

Statistical assessment

The statistical significance of the concentration changes of the identified metabolites between sample groups was finally assessed. With this aim, the chromatographic peak areas of these metabolites were integrated for every sample. These peak areas can be obtained from MCR-ALS or using Xcalibur™ 2.2.44 software (Thermo Scientific, Hemel Hempstead, UK). In the case of this work, the original software was used and the accurate mass values obtained from MCR-ALS were extracted. The obtained values were normalized by diving by the value of the chromatographic peak area of the internal standard (PIPES) of the considered sample. One-way ANOVA was applied to the normalized areas for each metabolite with a level of significance equal to 0.05. The statistical significance of the enrichment in members from different KEGG metabolic pathways was tested using the hypergeometric distribution. In both cases, Bonferroni post hoc tests were applied to correct for multiple comparisons. IBM® SPSS® Statistics 22.0.0.1 software was used for ANOVA calculations. Network analyses were performed in R software using the igraph package.⁴⁷

Results and discussion

Detection and identification of rice metabolites from an untargeted study is a challenging task owing to the complexity of the data, the wide range of existing metabolites and the fact that there is no prior knowledge about the chemical compounds (metabolites) contained in the analysed samples. For this reason, we propose the combined use of powerful analytical approaches with chemometric strategies. In this work PLS-DA and MCR-ALS chemometric methods have been used to seek rice metabolites changing their concentrations due to metal treatment from an untargeted study.

Assessment of metal treatment in plants

Table 1 gives the Cd(ii) and Cu(ii) concentrations in the aerial part samples analysed using the ICP-MS method explained in the section entitled “Analysis of metal content in rice samples” (results for root samples are shown in Table S1 in the ESI†). Results show that Cu(ii) was already detected in the rice control samples, whereas the Cd(ii) concentration was negligible in these control samples (close to the limit of detection of the technique, 0.18 μg g⁻¹). The concentration of Cd(ii) in the aerial parts of the 1000 μM treated samples was 70 times bigger than in control samples, whereas the concentration of Cu(ii) in the aerial parts of the 1000 μM treated samples was only 1.9 times larger than in control samples. Therefore, it seems clear that Cd(ii) accumulation in the aerial parts of rice plants was much higher than Cu(ii) accumulation.

Table 1

Results obtained in the determination of Cd and Cu in the aerial parts of the analysed rice samples

Treatment	Content (μg g−1 dry weight)
	Cu(ii)	Cd(ii)
Control 1	18.1 ± 0.7	0.09 ± 0.02
Control 2	14.0 ± 0.6	<LOD
10 μM Cu	18.64 ± 0.04	0.27 ± 0.04
50 μM Cu	19 ± 1	<LOD
100 μM Cu	22.9 ± 0.1	<LOD
1000 μM Cu	26.0 ± 0.4	<LOD
10 μM Cd	23 ± 2	1.8 ± 0.1
50 μM Cd	17 ± 1	3.78 ± 0.09
100 μM Cd	22.0 ± 0.7	5.4 ± 0.2
1000 μM Cd	14.7 ± 0.8	6.3 ± 0.9

Limits of detection were 1.75 μg g⁻¹ for Cu and 0.18 μg g⁻¹ for Cd. Three replicates were analysed.

Analysis of chromatographic regions

Untargeted metabolomics chromatographic data can be extremely complex containing both useful and meaningless information for the purpose of the study. In order to simplify the analysis, a first selection of the chromatographic regions of interest (those regions with metabolites changing their concentration) is performed. As described in the Experimental section (see section entitled “Selection of relevant chromatographic regions”), PLS-DA of the TIC chromatograms and VIPs plot were used for this purpose. Five chromatographic regions were selected for Cd(ii) treatment (from 2.2 min to 7.2 min, from 12.2 min to 20.35 min, from 20.35 min to 22.3 min, from 22.3 min to 25.2 min and from 25.2 to 28.4 min) and six for Cu(ii) treatment (from 2.2 min to 4.8 min, from 4.8 min to 7.7 min, from 12.2 min to 19.5 min, from 23.1 min to 25.2 min and from 25.2 to 28.4 min). For each one of these chromatographic regions, a column-wise augmented data matrix was built for each metal treatment and analysed separately by MCR-ALS as detailed in the section entitled “Multivariate curve resolution by alternating least squares”.

Between 15 and 30 MCR-ALS components were resolved for each one of the total 11 chromatographic regions (5 for Cd(ii) and 6 for Cu(ii)) with explained variances (R²) larger of 98%. The total number of MCR-ALS components used to explain data variance and patterns of all these chromatographic regions were 115 for Cd(ii) treatment and 100 for Cu(ii) treatment. Nevertheless, not all of these MCR-ALS resolved components were assigned to individual metabolites. Some of them did not correspond to the true metabolite chromatographic peaks but to other noisy chromatographic contributions, such as background and solvent signals. Despite the high complexity of the untargeted LC-MS data set, due to the strong overlap among metabolite elution profile at the same retention times, MCR-ALS could properly resolve a large number of metabolites from the investigated rice samples.

Fig. 3 is an example of results of the application of MCR-ALS to the resolution of two strongly coeluted metabolites corresponding to the first selected chromatographic region (elution times between 2.2 and 7.2 min). Fig. 3A depicts the elution profiles of two strongly coeluted metabolites successfully resolved after MCR-ALS analysis. In the case of the upregulated example (red elution profile), there was an increase in the chromatographic peak heights from control samples to treated samples (with their maximum height at intermediate Cd(ii) concentration levels). In the case of the downregulated example (blue elution profile), there was a decrease in the chromatographic peak height when considering control and treated samples and the peak practically disappeared for the 1000 μM Cd(ii) treated samples. Fig. 3B gives the mass spectra of the two resolved metabolites and the m/z values of their diagnostic ions were used for their preliminary identification. In this case, these m/z values are not accurate mass values because MCR-ALS was applied to compressed data using the binning approach. The accurate mass of selected metabolites was checked afterwards in the experimental high-resolution LC-MS raw data, which conserved full precision and accuracy of four decimal points obtained by Q-Exactive (quadrupole-Orbitrap). For instance, in the case of upregulated metabolite (red mass spectrum), the diagnostic ion appeared at 319.00 m/z and had the accurate mass value of 319.0466 m/z in the original raw data. In the case of downregulated metabolite (blue mass spectrum), the MCR-ALS diagnostic ion at 193.05 m/z had an accurate mass value of 193.0509 m/z in the original raw data.

Fig. 3

An example of MCR-ALS results obtained in the analysis of LC-HRMS data from 40 rice samples under Cd(ii) treatment at 5 concentration levels. (A) Resolved elution profiles obtained for two coeluted metabolites. Red profile corresponds to a metabolite showing lower concentrations in the control than in the treated samples (upregulated metabolite), and blue profile is for one metabolite showing higher concentrations in the control than in the treated samples (downregulated metabolite). (B) Resolved mass spectra for the upregulated metabolite (in red) and the downregulated metabolite (in blue).

Once the 11 column-wise augmented data matrices (from the 11 different chromatographic regions) were analysed by MCR-ALS, the peak areas of the resolved elution profiles were statistically evaluated as described in the experimental section (see section entitled “Metabolite identification”). The main aim was to detect the possible metabolites whose concentrations change significantly due to the metal treatment. Fig. 4A depicts a Venn diagram showing the total number of metabolites detected by the MCR-ALS approach varying their concentration upon metal treatment. A total number of 77 candidates were exclusively found for Cd(ii) treatment, 20 solely for Cu(ii) treatment and 15 of them were common for both treatments.

Fig. 4

Venn diagram showing the number of metabolites selected for Cd(ii) treatment (purple), for Cu(ii) treatment (green) and commonly for both treatments. (A) Metabolites detected after MCR-ALS analysis. (B) Metabolites identified and confirmed after AIF analysis. (C) Metabolites showing a statistically significant change. (D) Statistically significant upregulated and downregulated metabolites for Cd(ii) and Cu(ii) treatment.

Metabolite identification

The accurate mass of the candidate metabolite components whose elution profile peak areas changed was compared with the exact mass values included in public databases such as MassBank,⁴⁰ Metlin,⁴¹ HMDB,⁴² and MetaCyc.⁴³ As an example, the tentative identification of the two MCR-ALS components shown in Fig. 3 is explained here in more detail. In the case of the upregulated metabolite (depicted in red in Fig. 3a and b), the accurate mass of the observed diagnostic ion was 319.0466 Da. This mass could be tentatively assigned with 2.1 ppm of relative error, to dihydromyricetin, whose deprotonated molecule ([M − H]⁻) has an exact mass of 319.0459 Da. In the case of the downregulated metabolite (depicted in blue in Fig. 3a and b), the accurate mass of the observed diagnostic ion was 193.0509 Da, which was tentatively assigned to ferulic acid (deprotonated molecule exact mass: 193.0501 Da) with a relative mass error of 4.1 ppm. A total number of 112 metabolites could be tentatively identified with a relative mass error lower than 5.0 ppm.

A further step in metabolite identification and confirmation was the comparison of the obtained experimental AIF mass spectra with their theoretical MS/MS spectra. As mentioned in the experimental section, entitled “LC-MS analysis”, the AIF scan mode of Q-Exactive (Orbitrap) allowed the product ion spectra of all detected ions to be obtained without the pre-selection of their precursor ions in the quadrupole. The AIF scan mode allowed metabolite identification through the product ion mass spectrum without the need to inject samples twice.

Fig. 5 depicts an example of how this AIF confirmation was performed. Fig. 5A represents the experimental AIF high-resolution mass spectrum corresponding to the chromatographic peak at 22.24 minutes. At this retention time, the corresponding full scan showed an ion at m/z 341.1092, which was tentatively identified as trehalose with 2.3 ppm of relative mass error. Fig. 5B corresponds to the MS/MS spectrum of trehalose obtained from a database (Massbank library), which shows the precursor ion at 341.1084 m/z and several product ions. Among these ions, the most intense were at 101.0245 and 179.0561 m/z. The experimental spectrum showed two major product ions at 101.0245 and 179.0564 m/z, which could be correlated with the corresponding theoretical product ions with less than 1.7 ppm of relative mass error. Taking all this information into account, trehalose was identified with 4.5 IPs (2 IPs earned for the HRMS precursor ion and 2.5 for the two product ions) and, therefore, the identification criteria recommended by Directive 2002/657/CE were fully accomplished.

Fig. 5

Identification of trehalose. (A) AIF mass spectrum at a retention time of 22.24 minutes. (B) MS/MS spectrum of trehalose obtained from the Massbank library (Compound ID PR100542).

Fig. 4B depicts a Venn diagram of finally identified and confirmed metabolites. As it can be seen, from the total number of 77 tentatively identified metabolites for Cd(ii) treatment, 66 were confirmed. From the total number of 20 metabolites detected for Cu(ii) treatment, 18 were verified. Moreover, from the total number of 15 metabolites present in both metal treatments, 13 were ratified. In summary, a total number of 97 metabolites were confirmed, which is 87% of the metabolites primarily selected by MCR-ALS analysis. Furthermore, most of the identified metabolites had 4.5 IPs of identification, accomplishing the requirements of Directive 2002/657/CE.³⁹ This is an important result of this study which represents an improvement in comparison to the results obtained in our previous work about Japanese rice under Cd(ii) and Cu(ii) stress,¹⁷ in which the confirmation of the metabolite identification was not possible due to instrumental limitations (no accurate measurement of product ions was available).

Tables S2–S4 in the ESI† show the final AIF identification results. On the one hand, Table S2 (ESI†) gives the common metabolites obtained by both metal (Cd(ii) and Cu(ii)) treatments. On the other hand, Tables S3 and S4 (ESI†) give the metabolites identified only for Cd(ii) and for Cu(ii) metal treatment, respectively. These tables show the accurate mass of the diagnostic ions of the selected metabolites, the name of the identified metabolites, the ions assigned to diagnostic ions, the relative mass errors and the product ions used for confirmation. As it was expected, diagnostic ions corresponding to most of the identified metabolites were the deprotonated ion ([M − H]⁻). Nevertheless, some metabolites also generated other ions frequently observed in LC-MS. For example, the loss of a water molecule ([M − H₂O − H]⁻) or the formation of adduct ions with mobile phase components ([M − H + HAc]⁻) or metal cations in the LC-MS system ([M − 2H + K]⁻) or ([M − 2H + Na]⁻). It should be highlighted that in all cases, the relative error in the accurate mass value was lower than 5 ppm, accomplishing Directive 2002/657/CE requirements.³⁹ Furthermore, a minimum of two product ions were formed and detected for the major part of these identified metabolites, achieving therefore the four mandatory IPs for positive identification. Only 27 of the total 112 metabolites did not attain this recommendation: for 15 of them fragment ions with an m/z value lower than 90 were predicted, but they could not be detected under the mass scanning range used in this study. Also, 12 others showed only one predicted product ion. These 27 metabolites have been marked on Tables S2–S4 (ESI†) with an asterisk, (*). All these results confirm the potential of the proposed untargeted approach as a tool to gather the list of the most interesting metabolites showing changes due to the stress conditions of the study, in this case, Cd(ii) and Cu(ii) of the rice samples.

Finally, one-way ANOVA was applied to the chromatographic peak areas of the identified metabolites in order to investigate which of these metabolites presented statistically significant differences between sample groups treated by metal ions. A total number of 48 metabolites showed a significant change only for Cd(ii) treatment, 17 exclusively for Cu(ii) treatment and 6 commonly for both metal treatments, with a p-value lower than 0.05. Tables 2–4 show the metabolites that showed significant changes. On the one hand, Table 2 displays the 6 common metabolites changing significantly with both metal (Cd(ii) and Cu(ii)) treatments. On the other hand, Tables 3 and 4 give the metabolites showing a significant change only for Cd(ii) and for Cu(ii) metal treatment, respectively. Tables 2–4 give p-values (obtained after a Bonferroni post hoc test) for all those metabolites that had significant variation in their concentrations (peak areas), the value of the fold change for 10, 50, 100 and 1000 μM treated samples and if the identified metabolite is up or down regulated. The major part of the metabolites presented a similar fold change for the four levels of treatment (10, 50, 100 and 1000 μM). This result probably reflects that at the lowest concentration of Cd(ii) and Cu(ii) tested in this study (10 μM) the metabolites were already significantly affected, without any further modification at higher metal concentrations (50, 100 and 1000 μM). However, in some cases, appreciable differences between treatment levels were observed. For instance, in the case of synaptic acid for copper treatment (Table 2) the fold change decreased gradually when the level of treatment increased, which means that the concentration of synaptic acid in rice is lower when the concentration of Cu(ii) increased. Another relevant example is the case of o-feruloylquinic acid for cadmium (Table 3), where the fold changes were similar for the three lower levels of treatment (10, 50 and 100 μM) but it was almost doubled for 1000 μM treated samples, indicating a large effect of this high metal concentration. Another interesting point is that 11 metabolites showing a significant change for Cd(ii) metal treatment (Table 3), presented a fold change at 100 μM level treatment lower than at the other three levels. For instance, this is observed for glutamine or 1-O-vanilloyl-beta-d-glucose. This different behaviour might indicate an adaptive plant response at this intermediate concentration for this reduced number of metabolites. Further work should be done to corroborate this hypothesis, focusing particularly on the pathways related to these compounds.

Table 2

Metabolites presenting statistically significant differences between sample groups treated with Cd(ii) and Cu(ii)

Exact mass	Compound name	KEGG	Ion assignment	Rel. mass error (ppm)	AIF (m/z)	Cd-treatment						Cu-treatment
						Corrected p-value	Up/down-regulated	Fold change control vs.				Corrected p-value	Up/down-regulated	Fold change control vs.
								10	50	100	1000			10	50	100	1000
305.0184	UMP	C00105	[M − H2O − H]−	3.1	111.0201/211.0014	5.461 × 10−3	Down	0.24	0.44	0.42	0.67	3.57 × 10−5	Down	0.26	0.23	0.19	0.16
223.0611	Synaptic acid	C00482	[M − H]−	0.4	208.0379/93.0347/121.0297/149.0246/164.0481/193.0146	5.115 × 10−3	Up	1.89	1.84	2.18	1.65	7.90 × 10−3	Down	0.79	0.54	0.41	0.38
328.0457	Cyclic-AMP	C00575	[M − H]−	1.5	134.0054/107.0504/192.9911	3.407 × 10−3	Down	0.24	0.44	0.35	0.61	3.15 × 10−5	Down	0.25	0.20	0.15	0.13
157.0508	2-Isopropylmaleate	C02631	[M − H]−	0.9	127.0039/112.0531/123.0453	5.97 × 10−3	Up	2.40	2.58	1.07	3.34	3.07 × 10−3	Up	1.45	1.66	3.18	1.52
327.2186	2,3-Dinor-8-iso prostaglandin F1alpha	C14795	[M − H]−	2.7	310.2150/293.2122/125.0972	1.49 × 10−4	Down	0.18	0.51	0.32	0.46	2.94 × 10−4	Down	0.20	0.13	0.18	0.09
481.2577	Stearyl citrate	—	[M − 2H + K]−	0.7	190.0149/174.0171/268.2727	1.430 × 10−3	Down	0.04	0.12	0.07	0.13	9.12 × 10−5	Down	0.52	0.29	0.18	0.17

Table 3

Metabolites presenting statistically significant differences between sample groups treated with Cd(ii)

Exact mass	Compound name	KEGG	Ion assignment	Rel. mass error (ppm)	AIF (m/z)	Corrected p-value	Up/down-regulated	Fold change control vs.
								10	50	100	1000
176.9363	Diphosphate	C00013	[M − H]−	2.1	133.0146/114.9489	2.21 × 10−4	Down	0.49	0.45	0.78	0.35
145.0143	2-Oxoglutaratea	C00026	[M − H]−	0.5	101.0246	8.11 × 10−3	Up	1.17	1.32	1.15	1.05
145.0619	Glutamine	C00064	[M − H]−	0.2	127.0513/128.0355/109.0409	3.603 × 10−3	Up	3.77	3.02	1.64	4.33
104.0355	Serinea	C00065	[M − H]−	1.9	—	1.90 × 10−2	Up	2.20	1.96	1.45	2.46
171.0068	Glycerol 3-phosphate	C00093	[M − H]−	2.1	96.9697/152.9597	9.85 × 10−6	Down	0.39	0.39	0.48	0.51
259.0225	Glucose 1-phosphate	C00103	[M − H]−	0.1	96.9155/138.9805/181.0511/240.9522	1.36 × 10−2	Down	0.65	0.59	0.66	0.59
118.0511	Threoninea	C00188	[M − H]−	1.0	—	6.53 × 10−3	Up	2.78	2.46	1.38	3.47
105.0195	Glyceratea	C00258	[M − H]−	1.5	—	3.31 × 10−2	Down	0.98	0.62	1.21	0.81
299.0990	d-Ribulose	C00309	[2M − H]−	2.1	149.0459/133.0509	6.11 × 10−3	Down	0.29	0.64	0.46	0.44
671.4685	PA(16:0/18:2(9Z,12Z))	C00416	[M − H]−	4.2	391.2256/255.2329/433.2364/409.2363/415.2257/279.2329	1.20 × 10−2	Down	0.32	0.43	0.67	0.29
173.0458	Shikimate	C00493	[M − H]−	1.2	93.0347/99.0453/111.0453/137.0245/155.0352	1.82 × 10−2	Down	0.87	0.61	1.14	0.60
133.0143	Malic acida	C00711	[M − H]−	0.6	115.0038	3.04 × 10−3	Down	0.77	0.88	0.93	0.47
213.1498	(−)-Menthone	C00843	[M − H + HAc]−	0.9	111.0453/152.9959/138.0199	1.72 × 10−2	Down	0.95	0.51	1.17	0.64
127.0517	5,6-Dihydrothyminea	C00906	[M − H]−	2.7	—	4.08 × 10−3	Up	3.36	2.66	1.58	3.69
344.0407	Cyclic-GMP	C00942	[M − H]−	1.7	133.0151/150.0422	5.01 × 10−4	Down	0.23	0.42	0.34	0.53
239.0773	Galactose	C00984	[M − H + HAc]−	0.1	179.0565/161.0458	1.32 × 10−2	Down	1.01	0.92	0.99	0.95
311.1353	4-Hydroxybutanoatea	C00989	[3M − H]−	1.7	—	2.86 × 10−3	Up	2.01	3.41	2.81	3.18
244.0227	3-Phosphoserine	C01005	[M − H + HAc]−	0.4	184.0021/96.9697	1.89 × 10−2	Up	1.95	1.86	1.65	1.51
99.0089	4-Fumaryl-acetoacetate	C01061	[M − 2H]2−	1.7	100.0123/154.0273/112.0123/98.0011	1.76 × 10−2	Down	1.06	1.06	1.59	0.98
341.1092	Trehalose	C01083	[M − H]−	2.3	101.0245/113.0246/119.0305/143.0349/161.0456	3.68 × 10−3	Down	0.87	0.77	0.81	0.69
113.0610	2-Hydroxycyclohexan-1-onea	C01147	[M − H]−	1.6	99.0452	5.40 × 10−3	Down	0.50	0.46	0.80	0.46
193.0509	Ferulic acid	C01494	[M − H]−	1.1	134.0374/149.0609/178.0274	9.92 × 10−4	Down	0.35	0.45	0.54	0.49
371.0990	Syringin	C01533	[M − H]−	1.8	209.0820/433.1508	3.81 × 10−2	Up	2.13	2.14	1.95	1.93
157.0369	Allantoin	C01551	[M − H]−	1.1	129.0196/114.0311	1.63 × 10−3	Up	4.65	4.17	2.56	5.78
337.0932	Columbamine	C01795	[M − H]−	0.8	275.0952/289.1108	1.97 × 10−2	Down	1.09	0.77	1.60	0.87
128.0356	5-Oxoprolinea	C01879	[M − H]−	2.1	—	3.89 × 10−3	Up	3.61	2.97	1.13	4.56
367.1041	O-Feruloylquinic acid	C02572	[M − H]−	1.9	202.1086/122.0375/199.0615/174.0561/129.0559/192.0432	1.29 × 10−2	Up	1.82	1.83	1.63	2.67
319.0466	Dihydromyricetin	C02906	[M − H]−	2.1	194.0284/124.0168/177.0196/160.0169/143.0139/107.0140	1.49 × 10−3	Down	0.59	0.58	1.10	0.65
365.1092	cis-3,4-Leucopelargonidin	C03648	[M − H + HAc]−	4.9	109.02960/137.02453/289.07368	1.46 × 10−4	Down	0.67	0.42	1.00	0.40
274.0120	6-Pyruvoyltetrahydropterin	C03684	[M − 2H + K]−	1.2	192.0903/162.0489	2.36 × 10−2	Down	0.66	0.66	0.34	0.60
623.1642	Apigenin 7-O-beta-d-glucoside	C04608	[M − H + HAc]−	3.9	432.1024/431.0988/270.0420/269.0454/239.0345/151.0040	9.80 × 10−4	Down	0.47	0.57	1.28	0.72
337.0570	5-Amino-1-(5-phospho-d-ribosyl)imidazole-4-carboxamide	C04677	[M − H]−	4.5	96.96978/125.04167	9.75 × 10−3	Down	0.73	0.74	1.19	0.68
563.1421	Apigenin 7-O-[beta-d-apiosyl-(1->2)-beta-d-glucoside]	C04858	[M − H]−	2.6	269.0666/431.0988	1.04 × 10−3	Down	0.52	0.49	1.12	0.67
315.0717	3′-Hydroxy-N-methyl-(S)-coclaurine	C05202	[M − H]−	1.4	108.02183/122.0375/191.0955	5.19 × 10−3	Up	1.32	0.87	1.12	1.18
313.1145	(2R)-1-O-beta-d-Galactopyranosylglycerol	C05401	[M − H + HAc]−	1.7	259.0945/236.0954/191.0509/162.0499	1.43 × 10−5	Down	0.25	0.56	0.37	0.36
191.0560	l-Quinic acida	C06746	[M − H]−	0.6	96.96971	3.96 × 10−3	Down	0.76	0.50	0.99	0.51
134.0374	Medicarpin	C10503	[M − 2H]2−	0.5	269.1031/254.0956	1.18 × 10−2	Down	0.49	0.45	0.66	0.48
815.2290	Cyanidin-3-O-rutinoside-5-O-β-d-glucoside	C12646	[M − H + HAc]−	4.8	93.0347/177.0197/163.0616/179.0565/149.0457/133.0508	1.80 × 10−6	Down	0.22	0.25	0.86	0.24
251.1037	2,6-Dihydroxy-N-methylmyosmine	C16151	[M − H + HAc]−	0.1	109.01709/157.03167/176.67197	3.29 × 10−3	Down	0.38	0.53	0.49	0.63
329.0874	1-O-Vanilloyl-beta-d-glucose	C20470	[M − H]−	1.1	179.0565/149.0457/165.0406/93.0347	2.92 × 10−2	Up	9.52	7.42	0.93	61.69
793.5193	1-18:1-2-16:0-monogalactosyldiacylglycerol	—	[M − 2H + K]−	4.9	255.2329/295.2643	4.63 × 10−4	Down	0.46	0.67	0.54	0.36
227.1290	6(E)-8-Oxogeraniol	—	[M − H + HAc]−	0.4	167.1079/152.0845/137.0609/151.1009	3.28 × 10−3	Down	0.44	0.50	0.66	0.42
799.2141	6′′′-O-Sinapoylsaponarin	—	[M − H]−	4.8	92.0256/131.0352/142.0053/756.1915	4.51 × 10−3	Down	0.28	0.30	1.68	0.66
277.0333	Caffeoylmalic acid	—	[M − H2O − H]−	4.9	134.0351/278.0457/227.0337/186.0163/141.0182/116.0107/178.0240	2.68 × 10−5	Down	0.36	0.43	0.51	0.34
785.2200	Kaempferol 3-(2G-apiosylrobinobioside)	—	[M − H + HAc]−	4.9	93.0347/177.0198/131.0351/147.0303/117.0194/101.0245	3.24 × 10−5	Down	0.35	0.28	0.85	0.28
347.0599	Maclurin 3-C-(2′′-galloyl-6′′-p-hydroxybenzoyl-glucoside)	—	[M − 2H]2−	2.5	695.1251/93.0347/109.0296/124.0166/125.0245/136.0166	2.50 × 10−8	Up	1000.0	1000.0	1000.0	1000.0
331.0554	S-7-Methylthioheptylhydroximoyl-l-cysteine	—	[M − 2H + K]−	1.1	293.0990/130.0876/277.0839/278.0763/219.0772	3.88 × 10−3	Up	9.23	7.69	3.23	9.55
385.0547	Shoyuflavone Aa	—	[M − H]−	4.7	133.0144	4.91 × 10−2	Down	0.72	0.43	0.82	0.39

a

Metabolites identified with less than four IPs.

Table 4

Metabolites presenting statistically significant differences between sample groups treated with Cu(ii)

Exact mass	Compound name	KEGG	Ion assignment	Rel. mass error (ppm)	AIF (m/z)	Corrected p-value	Up/down-regulated	Fold change control vs.
								10	50	100	1000
175.0251	Ascorbic acid	C00072	[M − H]−	1.8	113.0246/157.0367/115.07661/139.0880/130.28988	2.686 × 10−5	Up	5.62	7.40	36.71	22.14
134.0474	Adenine	C00147	[M − H]−	1.5	106.0664/107.0365	3.97 × 10−3	Down	0.52	0.52	0.43	0.42
266.0895	Adenosine	C00212	[M − H]−	0.0	134.0474/107.0365	1.37 × 10−3	Down	0.52	0.45	0.41	0.40
409.2365	1-Palmitoylglycerol 3-phosphate	C00416	[M − H]−	1.0	152.9959/255.2329/171.0065	1.35 × 10−4	Down	0.27	0.25	0.19	0.17
433.2365	LPA(0:0/18:2(9Z,12Z))	C00416	[M − H]−	1.1	278.2207/153.9993/169.9923/416.2294	2.20 × 10−5	Down	0.43	0.33	0.14	0.14
316.1173	Glycerophosphocholinea	C00670	[M − H + HAc]−	2.0	199.0341	4.36 × 10−3	Up	1.46	1.47	1.37	1.62
226.9966	3-Dehydroquinate	C00944	[M − 2H + K]−	1.0	127.0403/171.0302/109.0295/153.0194/125.0245/189.0199/117.0348	1.260 × 10−4	Up	1.80	1.65	1.62	2.51
267.1084	2′-Deoxyribose	C01801	[2M − H]−	0.5	133.0509/103.0403	5.405 × 10−4	Up	24.43	22.25	30.03	14.04
431.2209	6-Aminopenicillanate	C02954	[2M − H]−	4.4	215.1078/171.0066	4.893 × 10−5	Down	0.45	0.31	0.14	0.14
179.0387	Methyl-5-thio-d-ribose	C03089	[M − H]−	2.2	161.0403/144.0305/130.0229/116.0072/102.0279	7.104 × 10−3	Up	2.03	1.85	2.21	1.15
102.0563	GABAa	C03665	[M − H]−	2.4	—	2 × 10−3	Down	0.49	0.50	0.49	0.56
312.0954	5-(3′-Carboxy-3′-oxopropenyl)-4,6-dihydroxypicolinate	C05641	[M − H + HAc]−	1.3	99.0088/151.0262/219.0161	1.936 × 10−3	Down	0.48	0.42	0.34	0.34
299.0776	(1R,6R)-6-Hydroxy-2-succinylcyclohexa-2,4-diene-1-carboxylate	C05817	[M − H + Hac]−	1.1	239.0561/137.0246/101.0246/222.0532/194.0588	3.914 × 10−7	Up	0.80	1.36	5.65	5.64
242.0800	Pantothenol	C05944	[M − 2H + K]−	0.2	102.0562/126.0938	1.60 × 10−3	Up	1.44	1.45	1.34	1.62
336.0875	S-(Hydroxymethyl)glutathione	C14180	[M − H]−	1.3	320.0678/319.0844/305.0689/277.0738/185.0571	3.201 × 10−7	Up	3.65	5.67	13.53	9.87
503.2417	PG(18:4(6Z,9Z,12Z,15Z)/0:0)	—	[M − H]−	1.4	152.9959/90.0350	1.990 × 10−5	Down	0.70	0.41	0.24	0.24
459.2337	Fusicoplagin Aa	—	[M − 2H + Na]−	4.8	347.2593	4.90 × 10−5	Down	0.11	0.03	0.07	0.03

a

Metabolites identified with less than four IPs.

Analysis of biological effects

A subset of 71 metabolites was identified as significantly affected by at least one of the treatments (see Venn diagram in Fig. 4D), with 60 of them annotated as bona fide O. sativa metabolites in KEGG (Tables 2–4). KEGG pathway analysis combining the two treatments detected 11 O. sativa pathways with at least four metabolites significantly affected by either Cd(ii) or Cu(ii) exposure (Table 5). Hypergeometric distribution analyses indicated that the number of affected metabolites was higher than expected by a random distribution (metabolite enrichment analysis, asterisks in Table 5) for many of these pathways, although the analysis is limited by the difficulty in determining which metabolites could be detected in our analyses. In any case, secondary metabolism, glycerolipid and glycerophospholipid metabolism, and some components of the carbon and amino acid metabolism seem to be affected by at least one of the treatments (Table 5).

Table 5

KEGG pathway analysis of rice metabolites with significant changes in concentration upon Cd(ii) or Cu(ii) exposure

KEGG pathway	Description	Metabolites	Input	Hit	Totala	Total hits	p-valueb
Cd exposure
osa01100	Metabolic pathways	C00026, C00064, C00065, C00093, C00103, C00105, C00188, C00258, C00309, C00416, C00482, C00493, C00984, C00989, C01005, C01061, C01083, C01494, C03684, C03692, C04677, C05202	40	23	2558	1629	0.092
osa01110	Biosynthesis of secondary metabolites	C00026, C00065, C00093, C00103, C00188, C00258, C00416, C00482, C00493, C00843, C01083, C01494, C01795, C02631, C02906, C03648, C04677, C05202, C10503	40	19	2558	648	0.001**
osa01230	Biosynthesis of amino acids	C00026, C00064, C00065, C00188, C00493, C01005	40	6	2558	147	0.018*
osa00230	Purine metabolism	C00064, C00212, C00147, C00575, C00942, C04677	40	5	2558	105	0.017*
osa02010	ABC transporters	C00064, C00065, C00093, C00188, C01083	40	5	2558	95	0.012*
osa00561	Glycerolipid metabolism	C00093, C00258, C00416, C03692, C05401	40	5	2558	120	0.027*
osa01200	Carbon metabolism	C00026, C00065, C00258, C00989, C01005	40	4	2558	57	0.009**
osa00970	Aminoacyl-tRNA biosynthesis	C00064, C00065, C00188, C01005	40	4	2558	62	0.012*
osa00630	Glyoxylate and dicarboxylate metabolism	C00026, C00064, C00065, C00258	40	4	2558	62	0.012*
osa00260	Glycine, serine and threonine metabolism	C00065, C00188, C00258, C01005	40	4	2558	57	0.009**
Cu exposure
osa01100	Metabolic pathways	C00072, C00105, C00147, C00212, C00416, C00482, C00681, C00944, C03089, C05817	18	10	2558	1629	0.146
osa01110	Biosynthesis of secondary metabolites	C00072, C00416, C00482, C00681, C00944, C02631, C05817	18	7	2558	648	0.086
osa00230	Purine metabolism	C00416, C00670, C00681	18	3	2558	105	0.03*
osa00564	Glycerophospholipid metabolism	C00147, C00212, C00575	18	3	2558	43	0.003**

a

Total O. sativa KEGG annotations.

b

Hypergeometric distribution.

* Metabolite enrichment analysis.

Network analyses show the mutual correlations between affected metabolites by either treatment and the associated metabolic pathways (Fig. 6). The graph shows five distinct functional groups, related to amino acid, purine, and glycerolipid metabolism, together with a more general, secondary metabolism group (Fig. 6). The graph also shows that metabolites related to amino acid metabolism became mainly elevated in both Cd(ii)- and Cu(ii)-treated samples (marked in red and orange in Fig. 6, respectively), whereas the concentrations of purine- and glycerolipid-related metabolites (including glycerophospholipids) were mainly reduced by both treatments (blue, green and purple characters in Fig. 6). It also shows that, although the overlap between significantly changed metabolites by the two treatments is relatively small (Fig. 4D), both treatments affected similar pathways in essentially the same direction (red and orange on one side, blue, green and purple on the other, Fig. 6). Note that only one out of the 60 metabolites, synaptic acid, showed a divergent response to both exposures, increasing its concentration upon Cd exposure and decreasing it upon Cu exposure (Table 2, labelled in black in Fig. 6). The increase in amino acid metabolism (with some metabolites also annotated in the carbon metabolism pathway, Table 5) can be explained as a detoxification mechanism of plants and as a protection mechanism of cell constituents. On the other hand, the decrease in nucleotide and lipid concentrations may be related to a reduced growth rate and/or photosynthetic activity. Treated plants were significantly small and yellowish compared to controls (see Fig. S2 in ESI†), probably as a result of oxidative stress associated with metal poisoning. The same growth limitations may be related to the general decrease on secondary metabolism (Fig. 6). Finally, and considering metabolites not annotated in KEGG, the data shows changes in concentrations of glycosides (6′′′-O-sinapoylsaponarin, kaempferol 3-(2G-apiosylrobinobioside) and maclurin 3-C-(2′′′′-galloyl-6′′′′-p-hydroxybenzoyl-glucoside), Table 3), as well as of stearyl citrate (Table 2), related to their biosynthesis. Glycosides are very abundant in plants and they form strong complexes with metals (Cd(ii) and Cu(ii)), which explains the significant changes observed in their concentrations under metal treatment. This type of glycoside alteration has also been observed in our previous work¹⁷ and also in other studies and plants, such as in Arabidopsis thaliana.¹¹

Fig. 6

Network analysis showing the mutual correlations between affected metabolites by both treatments and the associated metabolic pathways. Compounds and pathways are represented by their KEGG IDs, CXXXXX numbers for compounds (Table 3), and osaXXXXX labels for metabolic pathways (Table 4). Red and orange labels correspond to compounds whose concentrations increased upon Cd and Cu treatments, respectively; blue and green ones represent those whose concentration decreased. Purple labels indicate compounds whose concentrations decreased in both treatments, whereas the single compound with a contradictory response is shown in black (see the text). Ovals approximately delimit clusters of compounds and pathways affecting specific biological functions, as depicted in the graph.

Conclusions

The identification and confirmation of metabolites whose concentrations change due to the effects of cadmium and copper treatment in Japanese rice were assessed. The MCR-ALS chemometric procedure was first applied to LC-HRMS data with the aim of resolving elution and spectra profiles of metabolites present in the analysed rice samples. Metabolites whose elution profile peak areas were changed by metal treatment were then further identified from the accurate mass values of the corresponding diagnostic ions. High-resolution mass spectrometric detection combined with MCR-ALS data treatment has proven to be a powerful tool for untargeted metabolomics studies where no previous knowledge about what metabolites are affected by the investigated treatment was available. High-resolution AIF data analysis allowed the confirmation of metabolite identity. Most of the identified metabolites had the mandatory score of 4.5 IPs of identification, accomplishing Directive 2002/657/CE requirements.

ANOVA statistical assessment of chromatographic peak areas revealed that concentration changes of 54 metabolites were statistically significant under Cd(ii) treatment and that concentration changes of 23 metabolites were statistically significant under Cu(ii) treatment. Cd(ii) treatment appeared more severe that Cu(ii) treatment as it caused more significant changes in rice plant metabolism. However, the affected pathways (secondary metabolism and amino acid-, purine-, carbon- and glycerolipid-metabolism) are essentially the same, suggesting an underlying similarity of the responses of the plant to both divalent cations. These responses are consistent with a reduction in plant growth and/or photosynthetic capacity and with the induction of defence mechanisms to reduce cell damage.

Acknowledgements

The research leading to these results has received funding from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 320737. The authors would like to thank Center for Research in Agricultural Genomics (CRAG) for kindly supplying Japanese rice seeds.

Footnotes

†

Electronic supplementary information (ESI) available. See DOI: 10.1039/c6mt00279j