Skip to main content

Gene identification and transcriptome analysis of low cadmium accumulation rice mutant (lcd1) in response to cadmium stress using MutMap and RNA-seq

Abstract

Background

Cadmium (Cd) is a widespread toxic heavy metal pollutant in agricultural soil, and Cd accumulation in rice grains is a major intake source of Cd for Asian populations that adversely affect human health. However, the molecular mechanism underlying Cd uptake, translocation and accumulation has not been fully understood in rice plants.

Results

In this study, a mutant displaying extremely low Cd accumulation (lcd1) in rice plant and grain was generated by EMS mutagenesis from indica rice cultivar 9311 seeds. The candidate SNPs associated with low Cd accumulation phenotype in the lcd1 mutant were identified by MutMap and the transcriptome changes between lcd1 and WT under Cd exposure were analyzed by RNA-seq. The lcd1 mutant had lower Cd uptake and accumulation in rice root and shoot, as well as less growth inhibition compared with WT in the presence of 5 μM Cd. Genetic analysis showed that lcd1 was a single locus recessive mutation. The SNP responsible for low Cd accumulation in the lcd1 mutant located at position 8,887,787 on chromosome 7, corresponding to the seventh exon of OsNRAMP5. This SNP led to a Pro236Leu amino acid substitution in the highly conserved region of OsNRAMP5 in the lcd1 mutant. A total of 1208 genes were differentially expressed between lcd1 and WT roots under Cd exposure, and DEGs were enriched in transmembrane transport process GO term. Increased OsHMA3 expression probably adds to the effect of OsNRAMP5 mutation to account for the significant decreases in Cd accumulation in rice plant and grain of the lcd1 mutant.

Conclusions

An extremely low Cd mutant lcd1 was isolated and identified using MutMap and RNA-seq. A Pro236Leu amino acid substitution in the highly conserved region of OsNRAMP5 is likely responsible for low Cd accumulation in the lcd1 mutant. This work provides more insight into the mechanism of Cd uptake and accumulation in rice, and will be helpful for developing low Cd accumulation rice by marker-assisted breeding.

Background

Cadmium (Cd) is a widespread toxic soil pollutant, which is easily taken up by crop plants and accumulates in the edible parts, thus posing potential health risks for humans [1]. Rice (Oryza sativa L.), as staple food for half of the world’s population, is a major source of Cd intake for Chinese and Japanese populations, accounting for 40–50% of the daily Cd intake [2,3,4]. Thus, it is urgent to develop strategies to reduce rice grain Cd accumulation and the associated health risk of Cd intake. Selection and breeding of low Cd-accumulation rice varieties is one of the most economical and effective approaches to alleviating Cd toxicity and minimize Cd accumulation in rice grains grown on Cd-contaminated soil.

There are many transport processes involved in Cd accumulation in rice grains, including root uptake, root-to-shoot translocation by xylem, transfer from xylem to phloem and remobilization of Cd from leaves to grains [5,6,7]. Previous studies have identified several loci controlling these processes in rice plants by using quantitative trait locus (QTL) analyses. Three QTLs on chromosomes 3, 6, and 8 were identified as influencing Cd concentration in brown rice [8]. A major QTL qlGCd3 for low grain Cd content was identified on chromosome 3 [9]. Two QTLs, qGCd7 and qCdT7, for grain Cd accumulation and Cd translocation, respectively, were detected on chromosome 7 [10, 11]. Recently, several transporters involved in these processes have been characterized in rice plants. Natural resistance-associated macrophage protein (NRAMP) family proteins play an important role in entry of Cd into rice root cells [12, 13]. OsNRAMP5 was mainly expressed in the plasma membrane of rice root and participated in manganese (Mn) and Cd uptake. Osnramp5 mutant resulted in significantly reduced Mn and Cd concentrations in rice straw and grain [14]. OsNRAMP1 showed transport activity for Fe and Cd uptake, and was more highly expressed in high Cd-accumulating rice than low Cd-accumulating one [15]. Zinc-regulated transporter/iron-regulated transporter-like protein (ZIP) proteins and pleiotropic drug resistance (PDR)-type ATP-binding cassette (ABC) transporters have been also reported to be involved in Cd uptake and efflux in plants. Iron (Fe) deficiency-inducible expressions of OsIRT1 and OsIRT2 enhanced Cd uptake and translocation in rice [16]. AtPDR8 conferred Cd tolerance in Arabidopsis by pumping Cd2+ out across the plasma membrane of root epidermal cells [17]. OsABCG36/OsPDR9 expression was significantly induced in rice roots by Cd exposure [18]. In addition, P-type ATPase heavy metal transporters (HMA) have been identified as the regulator for xylem Cd transport in rice by mediating vacuolar sequestration of Cd in rice roots [19]. Over-expression of OsHMA3 in rice plants enhanced Cd tolerance and reduced Cd accumulation in rice grain. These studies provide deep insight into the genetic basis of Cd uptake and translocation in rice. However, the comprehensive expressions of these genes in low Cd-accumulation rice and its relationship to Cd accumulation have been little investigated in rice, especially for indica rice genotypes.

In past decades, microarray analysis has been a powerful tool providing important information on plants response to heavy metal stresses at the transcriptomic level. The comprehensive transcriptomic changes have been reported by various microarray analyse in rice and barley in response to Cd, copper (Cu) and mercury (Hg) stress [20,21,22]. These microarrays analyses improved our understanding of heavy metal-responsive gene expression in plants. However, microarrays are only applicable for detecting known transcripts and always miss low abundance genes, leading to a loss of some important candidate gene information. Recently, high-throughput RNA sequencing (RNA-seq) approach has become increasingly popular in transcriptomics studies, with increased sequencing depth and high detection sensitivity. To date, this technique has been applied in identifying the genes related to heavy metal responses in rice plants. Liu et al. [23] and He et al. [24] detected arsenite-responsive miRNAs and Cd toxicity-responsive genes in rice roots using RNAseq-based approach, and found that miRNAs and post-transcriptional alternative splicing may be involved in rice response to As and Cd exposure, respectively. Tan et al. [18] revealed universal Cd-responsive genes in rice roots exposed to various Cd stresses. Differential long non-coding RNAs were detected at an early stage in the rice response to Cd stress, suggesting their role in response to Cd stress [25]. Although these studies have identified several genes and transcription factors involved in response to heavy metal stresses in rice, the molecular mechanism underlying Cd uptake, translocation and accumulation has not been fully understood in low Cd accumulation rice.

In this study, 9311 and its mutant with extremely low Cd content in rice grains were used to investigate the differences in Cd uptake and accumulation in rice seedlings exposed to Cd stress. The candidate gene of low Cd accumulation mutant was identified by MutMap method and the comprehensive transcriptional changes related to Cd accumulation were determined using RNA-seq. Moreover, the tissue specificity of the candidate gene and its temporal response to Cd exposure in lcd1 and WT roots were examined by quantitative real time RT-PCR. The objective of this work is to clarify the physiological and molecular mechanism in response to Cd stress in low Cd accumulation rice, which will be helpful for developing low Cd accumulation rice by marker-assisted breeding.

Methods

Plant culture and treatment

Seeds of indica rice (cv. Yangdao 6, namely 9311), provided by the Jiangsu Lixiahe Agricultural Research Institute, were subjected to ethyl methanesulfonate (EMS) mutagenesis. Approximately 450 g seeds were immersed in distilled water for 12 h, and then incubated in 1% EMS for 12 h. After the EMS treatment, the seeds were washed in distilled water. The low Cd accumulation mutant (namely lcd1) in rice grains was identified according to the grain Cd concentration of 35,689 M2 plants determined by inductively coupled plasma-mass spectrometry (ICP-MS; X Series 2, Thermo Fisher Corp., Waltham, MA, USA), and was further self-pollinated to obtain M5 generation of lcd1. The mapping population was generated from the cross between the M5 generation of lcd1 and WT and was grown in Yoshida’s nutrient solution containing 0.1 μM Cd [26]. The F2 individuals with extremely high and low Cd accumulation in rice leaves were isolated at seedling stage and used as the mapping population for the MutMap method.

Field experiments were conducted during rice growing season from 2016 to 2017. Seeds of lcd1 and WT were sterilized in 10% NaClO for 20 min and washed with deionized water. Then seeds were germinated on moist filter paper at 37 °C and sown in sterilized moist quartz sand. Uniform 25-day-old seedlings were subsequently respectively transplanted into the Cd-contaminated experimental paddy field in Hangzhou (Field A, 1.5 mg kg− 1 of total Cd, pH 5.4), Quzhou (Field B, 2.6 mg kg− 1 of total Cd, pH 5.8), Yiyang (Field C, 4.5 mg kg− 1 of total Cd, pH 4.8) and Xiangtan (Field D, 0.35 mg kg− 1 of total Cd, pH 5.2). Rice plants grown in different experimental paddy fields were uniformly managed and harvested at tillering, heading, filling and maturity stages, respectively. All samples were oven dried to constant weight at 70 °C for mineral element concentration analysis.

Hydroponics experiments were also conducted to investigate the effects of Cd stress on growth and Cd accumulation of lcd1 and WT rice seedlings. Sterilized rice seeds were germinated on moist filter paper at 37 °C and grown on a plastic mesh floating on half-strength Yoshida’s nutrient solution. Uniform 15-day-old seedlings were transferred into fresh nutrient solutions for different Cd concentration treatments (0 and 5 μM). Each treatment was conducted in five replicated runs. After different times (0, 3, 6, 12 h and 14 days) of Cd treatment, rice plants were carefully washed with deionized water and separated into roots and shoots, immediately frozen in liquid nitrogen and then stored at − 80 °C until use.

Kinetic analyses of cd uptake in roots

The Cd uptake in roots was determined according to Uraguchi et al. [27] with some modifications. Uniform 15-day-old seedlings of lcd1 and WT were transferred to the uptake solution containing 0.5 mM CaCl2 and 2 mM MES (pH 5.6) for 24 h, and then carried out at room temperature for the uptake experiments.

For the dose-dependent Cd uptake experiment, rice seedlings were transferred to 1-L plastic container with the uptake solution containing different concentrations of CdCl2 (0, 0.1, 0.25, 0.5, 1, 5, 10, 15, 30 and 50 μM). Each treatment was replicated three times. After 1 h of uptake, plants were carefully washed with deionized water and separated into roots and shoots, oven dried and digested with HNO3 for Cd concentration analysis. Values of Km and Vmax for each genotype were calculated using software (GraphPadPRISM4; GraphPad Software Inc., CA, USA).

For the time-course Cd uptake experiment, rice seedlings were transferred to 1-L plastic container with the uptake solution containing 5 μM CdCl2 for different time treatments (0, 5, 10, 30, 60, 90, 120 and 180 min). Each treatment was replicated three times. The concentration of Cd in roots was determined as described below.

Determination of metal elements concentrations

The metal elements (Cd, Fe, Mn, Zn and Cu) concentrations were determined by ICP-MS according to Cao et al. [28]. Briefly, dried rice root, shoot and grain samples (0.25 g) were digested in concentrated HNO3, and diluted to 50 mL by deionized water for metal elements determination.

MutMap analysis

DNAs for MutMap analysis were extracted from rice leaves with DNeasy Plant Mini Kit (Qiagen, Hilden, Germany). Two parents and two bulked DNA pools were used for MutMap analysis. The bulked DNA were prepared by mixing DNA equally from 31 F2 individuals with extremely high and low Cd accumulation in rice leaves respectively. DNA quality and concentration were checked using the NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® DNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). Sequencing libraries were generated using Truseq Nano DNA HT Sample preparation Kit (Illumina, USA) and were sequenced by Illumina HiSeq4000 platform (Illumina, CA, USA).

The raw data (WT parent: SRR8695240, lcd1 parent: SRR8695241, extremely high Cd accumulation pool: SRR8695238, extremely low Cd accumulation pool: SRR8695239) were filtered by phred quality score to remove adapter sequences and low-quality bases. The clean reads were aligned to the 9311 reference genome (ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/oryza_indica/dna/) using BWA software. Alignment files were converted to BAM files through SAMtools software and applied to GATK to identify reliable SNPs. The SNP-index of two pools were calculated using sliding window methods and the difference between the SNP-index of two pools was calculated as ΔSNP-index.

RNA-seq analysis

Total RNAs for RNA-Seq analysis were extracted from rice roots of lcd1 and WT under Cd exposure with RNeasy Plant kit (Qiagen, Hilden, Germany). RNA quality and concentration were determined by NanoPhotometer® spectrophotometer (IMPLEN, CA, USA) and Qubit® RNA Assay Kit in Qubit® 2.0 Flurometer (Life Technologies, CA, USA). Sequencing libraries were generated using NEBNext® Ultra™ RNA Library Prep Kit for Illumina® (NEB, USA) and were sequenced on an Illumina Hiseq platform (Illumina, CA, USA).

The raw sequences for WT (SRR8718745, SRR8718748 and SRR8718749, three biological replicates) and lcd1 (SRR8718744, SRR8718746 and SRR8718747, three biological replicates) libraries were filtered using phred quality score to remove adapter sequences and low-quality bases. In high-quality reads, the sequences were trimmed and mapped to the 9311 reference genome using Hisat2. The best match was used to annotate genes and assign gene ontology information. FPKM values (fragments per kilobase of exon model per million mapped reads) were calculated to estimate gene expression levels. False positive and false negative errors were corrected by calculating the FDR (false discovery rate) adjusted q-values. The DEGs (differentially expressed genes) between lcd1 and WT under Cd exposure were selected based on normalized read count, and fold change (lcd1/WT) ≥ 1.5 and q-value < 0.05 were used as threshold values by using the DESeq R package (1.18.0). Gene ontology (GO) enrichment analysis of DEGs between lcd1 and WT under Cd exposure was performed using the GOseq R-package with q-value < 0.05.

Quantitative real-time RT-PCR

The tissue specificity of the candidate gene and its temporal expression under Cd exposure were examined by quantitative real time RT-PCR (qPCR). The mixture of first-strand cDNA from 1 μg of total RNA served as the templates for qPCR analysis using ReverTra Ace® qPCR RT Master Mix (Toyobo, Japan). The primers used for qPCR were designed by Primer Premier 5.0 (Premier, Palo Alto, California, USA) as listed in Additional file 1: Table S1. The qPCR reactions were run on a QuantStudio® 3 (Applied Biosystems, Thermo Fisher Scientific, Waltham, MA, USA), with a total reaction volume of 25 μL containing 2 μL of diluted cDNA template, 12.5 μL of SYBR Green Realtime PCR Master Mix, 1 μL of 10 μM forward and reverse primer, and 8.5 μL of nuclease-free water. The following reaction conditions were applied: incubation at 95 °C for 60 s, 40 cycles of denaturation at 95 °C for 15 s, annealing at 55 °C for 15 s and extension at 72 °C for 45 s, followed by the dissociation stage for melt curve analysis (from 60 to 94 °C). All of the samples were run with three replicates, and three no-template controls (NTC) were included in every run to monitor possible DNA contamination. ACTIN-1 was used as a control to calculate the relative gene expression levels.

To optimize qPCR conditions for each primer pair, PCR amplification specificity and efficiency were examined according to the MIQE guidelines [29]. The specificity of the PCR amplification was checked with a melt curve analysis. The result showed that the amplified product yielded a single peak at the melting temperature (Tm) for each gene, indicating that primer pairs for qPCR are highly specific. The PCR amplification efficiency was assessed using 5-fold serial dilution calibration-curve. The slopes of calibration curves range from − 3.53 to − 3.48, with corresponding R2 and efficiencies of 0.988–0.997 and 92.1–93.9%, respectively (Additional file 2: Table S2). conforming that the PCR amplification efficiency for each primer pair meets the MIQE recommended range (90.0–105.0%). The samples were normalized using ACTIN-1 expression and the relative gene expression levels were analyzed using the 2(−ΔΔCT) method due to that the amplification efficiencies of the target and reference were similar (Additional file 2: Table S2).

Statistics

Data were expressed as mean ± standard error of at least three biological replicates for each sample. Statistical analyses were performed using the statistical software graphpad prism (Graphpad Software, San Diego, CA, USA). Statistical significance was assessed using Duncan’s multiple comparison at P < 0.05 level.

Results

Isolation and characterization of the lcd1 mutant

We isolated a mutant (lcd1) displaying extremely low Cd content in rice grain from an EMS-mutagenized indica rice 9311. The lcd1 mutant was identified in screening for grains Cd concentration from among 35,689 M2 plants grown on Cd-contaminated paddy field and was further self-pollinated to obtain M5 generation of lcd1. At maturity, there were no obvious differences in plant height, tiller number and panicle morphology between lcd1 and WT grown in different paddy fields (Fig. 1a, b, Additional file 3: Table S3). However, the concentrations of Cd in lcd1 grains from all paddy fields were extremely low, ranging from 0.02 to 0.13 mg·kg− 1, which was much lower than the maximum limit of 0.2 mg·kg− 1 for Cd in rice; whereas the corresponding Cd concentrations in WT grains were 1.02–4.44 mg·kg− 1 (Fig. 1c). The grain Cd concentrations in lcd1 decreased by 96.2–98.0% compared with those in WT. Similar results were observed for root and shoot Cd concentrations of lcd1 at tillering, heading and filling stages (Fig. 1d, e). The other metal concentrations, such as Zn, Cu and Fe in grains were not significantly different between lcd1 and WT grown in most paddy fields. However, (Mn concentrations in lcd1 grains from all paddy fields were significantly lower than those in WT (Additional file 4: Figure S1).

Fig. 1
figure 1

Plant morphology and Cd concentrations in lcd1 and WT (cv.9311) plants and grains. a, lcd1 and WT plants after heading, bar = 10 cm. b, lcd1 and WT panicles, bar = 5 cm. c, Cd concentrations in lcd1 and WT grains grown on Cd-contaminated fields. d and e, Cd concentrations in shoots and roots of lcd1 and WT at tillering, heading and filling stage. * and ** indicate significant difference at p < 0.05 and p < 0.01 level between lcd1 and WT, respectively. The same as below

Rice seedlings of lcd1 and WT were exposed to 5 μM Cd in hydroponics. Under control condition, there were no visible differences in plant height, root length, root and shoot weights between lcd1 and WT. When grown in the presence of 5 μM Cd, both lcd1 and WT showed decreases in plant height and shoot weight compared with control, however, lcd1 shoots were much longer than WT in the presence of 5 μM Cd (Fig. 2a-d). Moreover, lcd1 accumulated less Cd in both shoots and roots compared with WT, irrespective of Cd treatment (Fig. 2e, f). The concentrations of Cu and Zn in both shoots and roots did not differ significantly between lcd1 and WT, however, the concentrations of Mn in both shoots and roots were significantly lower in lcd1 than in WT, while the shoot Fe concentration was relatively lower in lcd1 than in WT (Additional file 5: Figure S2).

Fig. 2
figure 2

Plant biomass and Cd concentrations in lcd1 and WT seedlings under Cd exposure. a-d, plant phenotype and biomass of lcd1 and WT seedling exposed to 5 μM CdCl2. E and F, Cd concentrations of the shoots and roots in lcd1 and WT seedling exposed to 5 μM CdCl2. Letters a-d indicate difference is significant at 0.05 level

The Cd uptake abilities of lcd1 and WT were further determined by a dose- and time-dependent Cd uptake in roots. The Cd uptake by roots was lower in lcd1 than in WT under all Cd treatment concentrations (Fig. 3a). The Vmax of lcd1 (162.00 ± 9.11 nmol g− 1 DW h− 1) was lower than that of WT (523.70 ± 33.27 nmol g− 1 DW h− 1), and the value of Km was greater in lcd1 (23.20 ± 2.83 μM) than in WT (9.26 ± 1.67 μM) (Table 1). The time-course of Cd uptake was 1.7–2.4 times higher in WT than in lcd1 through 180 min of Cd exposure (Fig. 3b). These results show that the Cd uptake ability of lcd1 is much lower than that of WT.

Fig. 3
figure 3

The dose- (a) and time-dependent (b) Cd uptake in roots of lcd1 and WT seedlings

Table 1 The Michaelis-Menten equation parameters of root Cd uptake in lcd1 and WT under different Cd exposure. Data are presented as means ± SE (n = 3)

Genetic analysis and identification of the candidate region of the lcd1 mutant by MutMap analysis

For rapid identification of the gene controlling the low Cd accumulation phenotype of the lcd1 mutant, we crossed the lcd1 mutant to WT 9311 in 2016 to obtain F1 progeny. F1 plants were self-pollinated to generate F2 progeny in 2017. 439 F2 individuals were grown in Yoshida’s nutrient solution containing 0.1 μM Cd. The progeny segregated in a 3:1 ratio (327:112) for WT phenotype and low Cd accumulation phenotype, respectively, indicating that the low Cd accumulation phenotype of the lcd1 mutant is caused by a single recessive mutation. We combined DNA form 31 F2 progeny that displayed extremely high (namely H2_H31 pool) and low Cd accumulation phenotype (namely H2_L31 pool) and subjected them to whole-genome sequencing using Illumina HiSeq4000 platform. We obtained 22,199,781,900 and 20,232,345,600 cleaned bases for Pool H2_H31 and H2_L31, respectively. These reads were aligned to the 9311 reference genome using the BWA software resulting in the identification of 1055 SNP positions.

For each SNP position, the value of SNP-index was obtained and a graph relating SNP positions and ΔSNP-index (the difference of the SNP-index of two pools) were generated for all 12 chromosomes of rice. As we expected, ΔSNP-index was distributed randomly around 0 for most parts of the genome (Additional file 6: Figure S3). Finally, we obtained the candidate region containing a cluster of SNPs with ΔSNP-index values ranging from 0.6 to 1 on chromosome 7 (Fig. 4a).

Fig. 4
figure 4

Identification of the lcd1 mutant by MutMap. a, Δ SNP-index plots of two pools on chromosome 7 generated by MutMap analysis. b, Structure of the OsNRAMP5 gene in lcd1 and WT. Gray boxes and black lines represent exons and introns in the OsNRAMP5 gene, respectively. SNP location is indicated below the gene structure with red font. c, Confirmation of SNP in the lcd1 mutant by Sanger sequencing. The SNP is indicated with red arrow. d, Homology comparison of the OsNRAMP5 protein in different crop plants. * indicates the identical amino acids, red box indicates the amino acid mutation location

Identification of the causal SNPs of the lcd1 mutant

There were nine SNPs with ΔSNP-index ranging from 0.6 to 1 identified in the candidate region, of which most SNPs located in non-coding region or had no gene annotation information (Table 2). A SNP at nucleotide position 8,887,787 corresponded to the seventh exon of the gene OsNRAMP5 (BGIOSGA024510, Os07g0257200), previously identified as an important transporter responsible for Mn and Cd uptake. This SNP led to a C to T transition at codon 236, resulting in replacement of a Pro (CCC) residue by Leu (CTC) in OsNRAMP5 (Fig. 4b). Sequence analysis confirmed the mutation site in the lcd1 mutant and multiple protein sequence comparison indicated that the residue affected was highly conserved in different crop plants (Fig. 4c, d), suggesting that the Pro236Leu substitution is likely responsible for the low Cd accumulation phenotype of the lcd1 mutant.

Table 2 Summary of SNPs showing ΔSNP-index from 0.6 to 1 within the candidate region

The tissue specificity and temporal expression level of OsNRAMP5 were further detected in lcd1 and WT under Cd exposure. OsNRAMP5 was expressed exclusively in WT root and little of its mRNA was detected in WT shoot (Fig. 5a). Moreover, OsNRAMP5 expression markedly increased with enhancing Cd treatment time in both WT and lcd1 roots compared with the control. However, lcd1 root showed a higher OsNRAMP5 expression at 3 and 6 h after Cd treatment, but significant lower expression at 12 h after Cd treatment than WT root (Fig. 5b), indicating a difference in Cd response pattern between lcd1 and WT root.

Fig. 5
figure 5

The tissue specificity of the OsNRAMP5 gene and its temporal response to Cd exposure in lcd1 and WT roots. a, Relative expression of OsNRAMP5 in WT root and shoot. b, Temporal responses of OsNRAMP5 in lcd1 and WT roots to Cd exposure

Transcriptome analysis of the lcd1 mutant under cd exposure by RNA-seq

To clarify how the mutation in OsNRAMP5 affects transcript levels of other genes involved in Cd accumulation in rice roots. Transcriptome of lcd1 and WT roots under Cd exposure were investigated by RNA-seq analysis. 15-day-old rice seedlings of lcd1 and WT were exposed to 5 μM Cd for 6 h, and root samples were used for sequencing on an Illumina platform. An average of ∼51.5 million clean reads per sample was attained from lcd1 and WT root cDNA libraries, among which more than 85% was uniquely mapped to the 9311 reference Genome (ftp://ftp.ensemblgenomes.org/pub/plants/release-36/fasta/oryza_indica/dna/). 25,090 expressed genes were detected in Cd-exposed WT root and 24,808 expressed genes in Cd-exposed lcd1 root, among which 24,102 genes were common ones, indicating similar genetic background between the two genotypes (Fig. 6a). The list of detected genes from RNA-Seq for lcd1 and WT under Cd exposure were provided in supplementary material (Additional file 7). A total of 1208 genes were differentially expressed between lcd1 and WT roots using threshold values of fold change (lcd1/WT) ≥ 1.5 and q-value < 0.05, including 577 up-regulated and 631 down-regulated ones (Fig. 6b). GO classification of biological process, cellular component and molecular function of the DEGs were all significantly enriched in transmembrane transport (Fig. 6c), indicating that transmembrane transport process plays an important role in differential Cd accumulation between lcd1 and WT roots.

Fig. 6
figure 6

Transcriptional changes between lcd1 and WT roots responsive to Cd exposure. a, Venn diagram of common Cd responsive genes in lcd1 and WT. b, Volcano map of differentially expressed genes (DEGs) between lcd1 and WT roots responsive to Cd exposure. c, Significantly enriched GO terms of DEGs between lcd1 and WT roots responsive to Cd exposure. GO terms belong to biological processes, molecular functions, and cellular components were shown in green, red, and blue, respectively. GO terms were sorted based on q-values (q-value < 0.05)

The changes in transcriptional levels of NRAMP, ZIP, HMA and ABC family genes, previously identified as important genes involved in Cd transmembrane transport, were further analyzed between lcd1 and WT (Table 3). In general, lcd1 root showed relatively lower expressions in NRAMP and ABC family genes, and higher expressions in ZIP and HMA family genes than WT root under Cd exposure. Among them, four genes (OsNRAMP1, OsABCG36/OsPDR9, OsZIP9 and OsHMA3) displayed larger relative changes in their expression levels between lcd1 and WT roots under Cd exposure compared with the other genes, indicating that they are involved in differential Cd accumulation between lcd1 and WT roots.

Table 3 Comparison of gene expressions related to Cd transmembrane transport between lcd1 and WT roots under Cd exposure using RNA-seq analysis

Discussion

Rice is a major dietary intake source of Cd for Asian populations, therefore, reducing rice grain Cd accumulation is very important for human health [2, 4]. Selection and breeding of low Cd-accumulation rice is recognized as the most economical and effective approach to reduce rice grain Cd concentration and associated Cd intake risk for humans [30]. In this study, we succeeded to obtain a non-transgenic low Cd accumulating rice mutant, lcd1, by EMS mutagenesis from 9311, a well-known commercial indica rice cultivar with its full genome sequence available. The lcd1 mutant had the similar agronomic traits to WT 9311 when grown on moderately or severely Cd-contaminated fields, but accumulated extremely low Cd in grains and plants (Additional file 3: Table S3, Fig. 1). Physiological studies in hydroponic culture further showed that lcd1 shoots were much longer than those of WT in the presence of 5 μM Cd, and had less Cd accumulation in shoots and roots compared with WT (Fig. 2). Moreover, root Cd uptake ability of lcd1 was much lower than that of WT (Fig. 3, Table 1). These results indicate that the lcd1 mutant has lower Cd-accumulating ability and higher tolerance to Cd stress than WT. It is noteworthy that the lcd1 mutant do not show any inhibition in rice growth compared with WT, although the concentrations of Mn were significantly lower in lcd1 rice seedling and grain than those in WT (Additional file 4: Figure S1, Additional file 5: Figure S2), suggesting the potential use of the lcd1 mutant in reducing intake risk of Cd from rice consumption without yield penalty.

For rapid identification of the gene controlling the low Cd accumulation phenotype of the lcd1 mutant, MutMap, a method based on whole-genome sequencing of bulked DNA of F2 segregating population was used in the present study. Compared with conventional map-based cloning, MutMap method can rapidly identify mutant genes and QTLs by using F2 progeny derived from the cross of the mutant and its wild type [31, 32]. Previous studies have identified some genes using MutMap method in rice. OsRR22, responsible for the salinity-tolerance phenotype for the hst1 mutant, has been identified by a MutMap method [31]. WB1, a regulator of endosperm development in rice, has been identified by a modified MutMap method [33]. In the present study, MutMap analysis, sequencing verification and homologous comparison results showed that a Pro236Leu amino acid substitution in a highly conserved region of OsNRAMP5 is likely responsible for the low Cd uptake and accumulation phenotype of the lcd1 mutant (Fig. 4). In line with our study, Ishikawa et al. [13] isolated three low Cd accumulation mutants from ion-beam irradiated seeds of japonica rice cultivar Koshihikari, namely lcd-kmt1, lcd-kmt2 and lcd-kmt3. Gene identification results showed that low Cd accumulation phenotype of three mutants were caused by 1-bp deletion, 433-bp insertion and ~ 277 kb deletion in OsNRAMP5 respectively. In the present study, we isolated the indica rice lcd1 mutant with the nonsynonymous mutation at codon 236 (CCC → CTC) in OsNRAMP5. The lcd1 mutant had similar extremely low Cd accumulation phenotype to the reported OsNRAMP5 mutants [13], but differed in the mutation site of the OsNRAMP5 gene, indicating that it may be a new OsNRAMP5 allelic mutant in rice. Furthermore, Tang et al. [34] and Sasaki et al. [14] reported that knockout of OsNRAMP5 using the CRISPR/Cas9 system and T-DNA insertion produced low Cd-accumulating indica rice (cv Huazhan) and japonica rice (cv Zhonghua11), respectively, indicating that OsNRAMP5 is a major transporter responsible for Cd uptake in rice. In the present study, a Pro236Leu amino acid substitution likely results in a defective OsNRAMP5 in the lcd1 mutant. Our results supported previous findings that OsNRAMP5 is a major transporter responsible for Cd uptake in rice. However, in the present study, the OsNRAMP5 expressions in lcd1 root were increased at 3 and 6 h after Cd treatment, but decreased at 12 h after Cd treatment than those in WT root (Fig. 5b). This might be due to the fact that low Cd concentration in the lcd1 mutant promotes Cd absorption by inducing OsNRAMP5 expression during the early stage of Cd exposure.

To clarify the molecular mechanism underlying Cd uptake, translocation and accumulation in the lcd1 mutant, transcriptome changes between lcd1 and WT under Cd exposure were investigated by RNA-seq analysis. RNA-seq provides an effective approach to identify the whole transcriptome changes in rice responsive to heavy metal stresses. Previous studies have demonstrated the transcriptional changes in rice plants responsive to Cd using high-throughput RNA sequencing [18, 24, 25]. Similar transcriptome investigations were performed in maize [35], barley [36] and bean (Phaseolus vulgaris L.) [37] under Cd, Zn and Cu stress. In this study, a total of 1208 DEGs were identified between lcd1 and WT roots under Cd exposure, and most were significantly enriched in transmembrane transport (Fig. 6), indicating that the transmembrane transport process contributes to differential Cd response and accumulation between lcd1 and WT roots. Hence, the changes in transcriptional levels of NRAMP, ZIP, HMA and ABC family genes, previously identified as important genes involved in Cd transmembrane transport, were further examined between lcd1 and WT (Table 3). Four genes, OsNRAMP1, OsABCG36/OsPDR9, OsZIP9 and OsHMA3 displayed larger relative changes in their expression levels between lcd1 and WT roots under Cd exposure compared with the other genes. OsNRAMP1 was involved in uptake and accumulation of Fe and Cd in rice. Over-expression of OsNRAMP1 in rice enhanced Cd accumulation in the leaves [15]. In this study, OsNRAMP1 expression was markedly decreased by 3.1-fold in lcd1 root compared with WT root at 6 h after Cd exposure, indicating that the low Cd uptake and accumulation in lcd1 root could be partly explained by the decreased expression of OsNRAMP1. Moreover, OsHMA3 has been identified as major regulator for xylem Cd transport in rice. OsHMA3 participates in the sequestration of Cd into vacuoles and affects root-to-shoot Cd translocation in rice. Overexpression of OsHMA3 enhanced Cd tolerance in rice [19], while loss-of-function allele of OsHMA3 resulted in high Cd accumulation in rice shoots and grain in some Indica rice cultivars [38]. In this study, OsHMA3 protein coding sequence in indica rice 9311 was identical to the type V alleles of OsHMA3 as reported by Yan et al. [38], which was proved to be functional by heterologous expression in yeast. There was no SNP and InDel detected in lcd1 OsHMA3 compared with WT, indicating no loss of function mutation of OsHMA3 in the lcd1 mutant. Furthermore, OsHMA3 expression was increased by 1.3-fold in lcd1 root compared with WT root at 6 h after Cd exposure, suggesting that increased OsHMA3 expression probably adds to the effect of OsNRAMP5 mutation to account for the significant decreases in Cd accumulation in rice plant and grain of the lcd1 mutant. In addition to OsHMA3, some PDR-type ABC transporters have been found to participate in Cd detoxification in plants. AtPDR8 conferred Cd tolerance in Arabidopsis by pumping Cd2+ out across the plasma membrane of root epidermal cells [17]. In this study, the expressions of four ABC transporter encoding genes were all significantly down-regulated in lcd1 root compared with WT root under Cd exposure, of which OsABCG36/OsPDR9 showed a 7.9-fold decreased expression. This might be due to that the lower Cd uptake and intracellular Cd level in the root cell of lcd1 caused weaker induction of OsABCG36/OsPDR9 expression in lcd1 root compared with in WT root. Actually, previous studies have demonstrated that OsPDR9 was differentially induced by Cd stress. Moons [39] suggested that OsPDR9 was rapidly induced by Cd stress in roots of rice seedlings. Tan et al. [18] also found that OsABCG36/OsPDR9 expression was remarkably induced by Cd in rice roots, which is basically consistent with our results, indicating the potential role of ABC transporters in regulating Cd tolerance in rice plants.

Availability of data and materials

All the data is contained within the manuscript.

Abbreviations

ABC:

ATP-binding cassette transporter

DEGs:

differential expressed genes

EMS:

ethyl methanesulfonate

GO:

gene ontology

HMA:

heavy-metal ATPase

NRAMP:

natural resistance-associated macrophage protein

PC:

phytochelatin

QTL:

quantitative trait locus

RNA-seq:

RNA sequencing

ZIP:

zinc-regulated transporter/iron-regulated transporter-like protein

References

  1. Nawrot TS, Staessen JA, Roels HA, Munters E, Cuypers A, Richart T, Ruttens A, Smeets K, Clijsters H, Vangronsveld J. Cadmium exposure in the population: from health risks to strategies of prevention. Biometals. 2010;23(5):769–82.

    Article  CAS  Google Scholar 

  2. Arunakumara KKIU, Walpola BC, Yoon MH. Current status of heavy metal contamination in Asia’s rice lands. Rev Environ Sci Bio. 2013;12(4):355–77.

    Article  CAS  Google Scholar 

  3. Wang L, Cui X, Cheng H, Chen F, Wang J, Zhao X, Lin C, Pu X. A review of soil cadmium contamination in China including a health risk assessment. Environ Sci Pollut R. 2015;22(21):16441–52.

    Article  CAS  Google Scholar 

  4. Ke S, Cheng XY, Zhang N, Hu HG, Yan Q, Hou LL, Sun X, Chen ZN. Cadmium contamination of rice from various polluted areas of China and its potential risks to human health. Environ Monit Assess. 2015;187(7):408.

    Article  Google Scholar 

  5. Uraguchi S, Fujiwara T. Rice breaks ground for cadmium-free cereals. Curr Opin Plant Biol. 2013;16(3):328–34.

    Article  CAS  Google Scholar 

  6. Clemens S, Aarts MG, Thomine S, Verbruggen N. Plant science: the key to preventing slow cadmium poisoning. Trends Plant Sci. 2013;18(2):92–9.

    Article  CAS  Google Scholar 

  7. Clemens S, Ma JF. Toxic heavy metal and metalloid accumulation in crop plants and foods. Annu Rev Plant Biol. 2016;67:489–512.

    Article  CAS  Google Scholar 

  8. Ishikawa S, Ae N, Yano M. Chromosomal regions with quantitative trait loci controlling cadmium concentration in brown rice (Oryza sativa). New Phytol. 2005;168(2):345–50.

    Article  CAS  Google Scholar 

  9. Abe T, Nonoue Y, Ono N, Omoteno M, Kuramata M, Fukuoka S, Yamamoto T, Yano M, Ishikawa S. Detection of QTLs to reduce cadmium content in rice grains using LAC23/Koshihikari chromosome segment substitution lines. Breeding Sci. 2013;63(3):284–91.

    Article  CAS  Google Scholar 

  10. Ishikawa S, Abe T, Kuramata M, Yamaguchi M, Ando T, Yamamoto T, Yano M. A major quantitative trait locus for increasing cadmium-specific concentration in rice grain is located on the short arm of chromosome 7. J Exp Bot. 2009;61(3):923–34.

    Article  Google Scholar 

  11. Tezuka K, Miyadate H, Katou K, Kodama I, Matsumoto S, Kawamoto T, Masaki S, Satoh H, Yamaguchi M, Sakurai K, Takahashi H, Satoh-Nagasawa N, Watanabe A, Fujimura T, Akagi H. A single recessive gene controls cadmium translocation in the cadmium hyperaccumulating rice cultivar Cho-Ko-koku. Theor Appl Genet. 2010;120(6):1175–82.

    Article  CAS  Google Scholar 

  12. Sasaki A, Yamaji N, Ma JF. Transporters involved in mineral nutrient uptake in rice. J Exp Bot. 2016;67(12):3645–53.

    Article  CAS  Google Scholar 

  13. Ishikawa S, Ishimaru Y, Igura M, Kuramata M, Abe T, Senoura T, Hase Y, Arao T, Nishizawa NK, Nakanishi H. Ion-beam irradiation, gene identification, and marker-assisted breeding in the development of low-cadmium rice. P Natl Acad Sci USA. 2012;109(47):19166–71.

    Article  CAS  Google Scholar 

  14. Sasaki A, Yamaji N, Yokosho K, Ma JF. Nramp5 is a major transporter responsible for manganese and cadmium uptake in rice. Plant Cell. 2012;24(5):2155–67.

    Article  CAS  Google Scholar 

  15. Takahashi R, Ishimaru Y, Senoura T, Shimo H, Ishikawa S, Arao T, Nakanishi H, Nishizawa NK. The OsNRAMP1 iron transporter is involved in cd accumulation in rice. J Exp Bot. 2011;62(14):4843–50.

    Article  CAS  Google Scholar 

  16. Nakanishi H, Ogawa I, Ishimaru Y, Mori S, Nishizawa NK. Iron deficiency enhances cadmium uptake and translocation mediated by the Fe2+ transporters OsIRT1 and OsIRT2 in rice. Soil Sci Plant Nutr. 2006;52(4):464–9.

    Article  CAS  Google Scholar 

  17. Kim DY, Bovet L, Maeshima M, Martinoia E, Lee Y. The ABC transporter AtPDR8 is a cadmium extrusion pump conferring heavy metal resistance. Plant J. 2007;50(2):207–18.

    Article  CAS  Google Scholar 

  18. Tan M, Cheng D, Yang Y, Zhang G, Qin M, Chen J, Chen Y, Jiang M. Co-expression network analysis of the transcriptomes of rice roots exposed to various cadmium stresses reveals universal cadmium-responsive genes. BMC Plant Biol. 2017;17(1):194.

    Article  Google Scholar 

  19. Sasaki A, Yamaji N, Ma JF. Overexpression of OsHMA3 enhances cd tolerance and expression of Zn transporter genes in rice. J Exp Bot. 2014;65(20):6013–21.

    Article  CAS  Google Scholar 

  20. Ding Y, Chen Z, Zhu C. Microarray-based analysis of cadmium-responsive microRNAs in rice (Oryza sativa). J Exp Bot. 2011;62(10):3563–73.

    Article  CAS  Google Scholar 

  21. Lopes MS, Iglesia-Turino S, Cabrera-Bosquet L, Serret MD, Bort J, Febrero A, Araus JL. Molecular and physiological mechanisms associated with root exposure to mercury in barley. Metallomics. 2013;5(9):1305–15.

    Article  CAS  Google Scholar 

  22. Lin CY, Trinh NN, Fu SF, Hsiung YC, Chia LC, Lin CW, Huang HJ. Comparison of early transcriptome responses to copper and cadmium in rice roots. Plant Mol Biol. 2013;81(4–5):507–22.

    Article  CAS  Google Scholar 

  23. Liu Q, Zhang H. Molecular identification and analysis of arsenite stress-responsive miRNAs in rice. J Agr Food Chem. 2012;60(26):6524–36.

    Article  CAS  Google Scholar 

  24. He F, Liu Q, Zheng L, Cui Y, Shen Z, Zheng L. RNA-Seq analysis of rice roots reveals the involvement of post-transcriptional regulation in response to cadmium stress. Front Plant Sci. 2015;6:1136.

    PubMed  PubMed Central  Google Scholar 

  25. Chen L, Shi S, Jiang N, Khanzada H, Wassan GM, Zhu C, Peng X, Xu J, Chen Y, Yu Q, He X. Genome-wide analysis of long non-coding RNAs affecting roots development at an early stage in the rice response to cadmium stress. BMC Genomics. 2018;19(1):460.

    Article  Google Scholar 

  26. Yoshida S, Forno DA, Cock JH, Gomez KA. Laboratory manual for physiological studies of rice. 3rd ed. Manila: International Rice Research Institute; 1976.

    Google Scholar 

  27. Uraguchi S, Mori S, Kuramata M, Kawasaki A, Arao T, Ishikawa S. Root-to-shoot cd translocation via the xylem is the major process determining shoot and grain cadmium accumulation in rice. J Exp Bot. 2009;60(9):2677–88.

    Article  CAS  Google Scholar 

  28. Cao ZZ, Mou RX, Cao ZY, Lin XY, Xu P, Chen ZJ, Zhu ZW, Chen MX. Nickel in milled rice (Oryza sativa L.) from the three main rice-producing regions in China. Food Addit Contam Part B. 2017;10(1):69–77.

    Article  CAS  Google Scholar 

  29. Bustin SA, Benes V, Garson JA, Hellemans J, Huggett J, Kubista M, Mueller R, Nolan T, Pfaffl MW, Shipley GL, Vandesompele J, Wittwer CT. The MIQE guidelines: minimum information for publication of quantitative real-time PCR experiments. Clin Chem. 2009;55(4):611.

    Article  CAS  Google Scholar 

  30. Ashrafzadeh S, Leung DW. Development of cadmium-safe crop cultivars: a mini review. J Crop Improv. 2016;30(2):107–17.

    Article  CAS  Google Scholar 

  31. Takagi H, Tamiru M, Abe A, Yoshida K, Uemura A, Yaegashi H, Obara T, Oikawa K, Utsushi H, Kanzaki E, Mitsuoka C, Natsume S, Kosugi S, Kanzaki H, Matsumura H, Urasaki N, Kamoun S, Terauchi R. MutMap accelerates breeding of a salt-tolerant rice cultivar. Nat Biotechnol. 2015;33(5):445.

    Article  CAS  Google Scholar 

  32. Abe A, Kosugi S, Yoshida K, Natsume S, Takagi H, Kanzaki H, Matsumura H, Yoshida K, Mitsuoka C, Tamiru M, Innan H, Cano L, Kamoun S, Terauchi R. Genome sequencing reveals agronomically important loci in rice using MutMap. Nat Biotechnol. 2012;30(2):174.

    Article  CAS  Google Scholar 

  33. Wang H, Zhang YX, Sun LP, Xu P, Tu RR, Meng S, Wu WX, Anis GB, Hussain K, Riaz A, Chen DB, Cao LY, Cheng SH, Shen XH. WB1, a regulator of endosperm development in Rice, is identified by a modified MutMap method. Int J Mol Sci. 2018;19(8):2159.

    Article  Google Scholar 

  34. Tang L, Mao BG, Li YK, Lv QM, Zhang LP, Chen CY, He HJ, Wang WP, Zeng XF, Shao Y, Pan YL, Hu YY, Peng Y, Fu XQ, Li HQ, Xia ST, Zhao BR. Knockout of OsNramp5 using the CRISPR/Cas9 system produces low cd-accumulating indica rice without compromising yield. Sci Rep. 2017;7(1):14438.

    Article  Google Scholar 

  35. Peng H, He X, Gao J, Ma H, Zhang Z, Shen Y, Pan G, Lin H. Transcriptomic changes during maize roots development responsive to cadmium (cd) pollution using comparative RNAseq-based approach. Biochem Bioph Res Co. 2015;464(4):1040–7.

    Article  CAS  Google Scholar 

  36. Kintlová M, Blavet N, Cegan R, Hobza R. Transcriptome of barley under three different heavy metal stress reaction. Genomics Data. 2017;13:15–7.

    Article  Google Scholar 

  37. Astudillo-Reyes C, Fernandez AC, Cichy KA. Transcriptome characterization of developing bean (Phaseolus vulgaris L.) pods from two genotypes with contrasting seed zinc concentrations. PLoS One. 2015;10(9):e0137157.

    Article  Google Scholar 

  38. Yan JL, Wang PT, Wang P, Yang M, Lian XM, Tang Z, Huang CF, Salt DE, Zhao FJ. A loss-of-function allele of OsHMA3 associated with high cadmium accumulation in shoots and grain of japonica rice cultivars. Plant Cell Environ. 2016;39(9):1941–54.

    Article  CAS  Google Scholar 

  39. Moons A. Ospdr9, which encodes a PDR-type ABC transporter, is induced by heavy metals, hypoxic stress and redox perturbations in rice roots 1. FEBS Lett. 2003;553(3):370–6.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We would like to thank Pro. Weixing Zhang (China National Rice Research Institute, Hangzhou, China) for proving the technical assistance.

Funding

This research was supported by the National Natural Science Foundation of China (No. 31701408, No. 31801304), the Zhejiang Provincial Natural Science Foundation of China (No. LQ17C130002) and the earmarked fund for China Agriculture Research System (CARS-01-47). These funding bodies played no role in the design of the study and collection, analysis, and interpretation of data and in writing the manuscript.

Author information

Authors and Affiliations

Authors

Contributions

ZC designed the research, performed the experiments and drafted the manuscript. XL and YY analyzed the data and helped in the preparation \of various parts of the manuscript. MG and PX provided support with instrumental analysis. MC supervised the study and contributed in writing the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Zhen Zhen Cao or Ming Xue Chen.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1:

Table S1. The primer sequences for quantitative real-time RT-PCR (DOCX 19 kb)

Additional file 2:

Table S2. The PCR amplification efficiency for each primer pair (DOCX 19 kb)

Additional file 3:

Table S3. Comparison of the main agronomic traits between lcd1 and WT plants (DOCX 15 kb)

Additional file 4:

Figure S1. Metal concentrations in lcd1 and WT grains grown on Cd-contaminated fields. (TIF 158 kb)

Additional file 5:

Figure S2. Metal concentrations in shoots and roots of lcd1 and WT seedlings exposed to 5 μM CdCl2. (TIF 86 kb)

Additional file 6:

Figure S3. Δ SNP-index plots of two pools on all 12 chromosomes generated by MutMap analysis. (TIF 494 kb)

Additional file 7:

The list of detected genes from RNA-Seq for lcd1 and WT under Cd exposure. co_exp: genes with the common expression for lcd1 and WT; specific_exp: genes specific for lcd1 and WT. (XLS 3494 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Cao, Z.Z., Lin, X.Y., Yang, Y.J. et al. Gene identification and transcriptome analysis of low cadmium accumulation rice mutant (lcd1) in response to cadmium stress using MutMap and RNA-seq. BMC Plant Biol 19, 250 (2019). https://doi.org/10.1186/s12870-019-1867-y

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12870-019-1867-y

Keywords