- Research article
- Open access
- Published:
Transcriptome analysis of phosphorus stress responsiveness in the seedlings of Dongxiang wild rice (Oryza rufipogon Griff.)
Biological Research volume 51, Article number: 7 (2018)
Abstract
Background
Low phosphorus availability is a major factor restricting rice growth. Dongxiang wild rice (Oryza rufipogon Griff.) has many useful genes lacking in cultivated rice, including stress resistance to phosphorus deficiency, cold, salt and drought, which is considered to be a precious germplasm resource for rice breeding. However, the molecular mechanism of regulation of phosphorus deficiency tolerance is not clear.
Results
In this study, cDNA libraries were constructed from the leaf and root tissues of phosphorus stressed and untreated Dongxiang wild rice seedlings, and transcriptome sequencing was performed with the goal of elucidating the molecular mechanisms involved in phosphorus stress response. The results indicated that 1184 transcripts were differentially expressed in the leaves (323 up-regulated and 861 down-regulated) and 986 transcripts were differentially expressed in the roots (756 up-regulated and 230 down-regulated). 43 genes were up-regulated both in leaves and roots, 38 genes were up-regulated in roots but down-regulated in leaves, and only 2 genes were down-regulated in roots but up-regulated in leaves. Among these differentially expressed genes, the detection of many transcription factors and functional genes demonstrated that multiple regulatory pathways were involved in phosphorus deficiency tolerance. Meanwhile, the differentially expressed genes were also annotated with gene ontology terms and key pathways via functional classification and Kyoto Encyclopedia of Gene and Genomes pathway mapping, respectively. A set of the most important candidate genes was then identified by combining the differentially expressed genes found in the present study with previously identified phosphorus deficiency tolerance quantitative trait loci.
Conclusion
The present work provides abundant genomic information for functional dissection of the phosphorus deficiency resistance of Dongxiang wild rice, which will be help to understand the biological regulatory mechanisms of phosphorus deficiency tolerance in Dongxiang wild rice.
Background
Phosphorus (P) is an essential mineral element required for plant growth and development. It plays a crucial role in the processes of energy transfer, signal transduction, photosynthesis and respiration [1]. Rice is one of the most important food crops in the world, and it is also considered as a model organism for monocotyledon genomics research [2]. Although the total amount of P in the soil may be high, it is often present in unavailable forms or in forms that are only available outside the rhizosphere [3]. Thus, the deficiency of phosphorus in soil is a worldwide problem. According to the statistics, about 5.8 billion hm2 of arable land worldwide is deficient in P. And about 67 million hm2 of arable land in China is deficient in P, which resulted in yield reduction by 5–15% (about 25–75 billion Kg) [4]. Over the past several decades, P fertilizer application for crop growth has been increased rapidly, but P-use efficiency has decreased to a low level of 10–20% [5]. P fertilizer consumes non-renewable phosphate rock reserves, which are expected to be exhausted in the near future [6]. Furthermore, much of the applied P has caused serious environmental pollutions [7]. Therefore, it is a significant development direction to improve the status of P-deficiency in crops to explore the absorption and utilization of high P-use efficiency of crops.
Previous studies have found that there are significant genotypic differences in P-deficiency tolerance among different genotypes [8] which make it possible to breed varieties with improved P-deficiency tolerance. During the course of domestication from wild rice to cultivated rice, however, the number of alleles of cultivated rice reduced by 50–60% compared to wild rice [9]. Moreover, the wide application of commercial hybrids in recent decades has caused the loss of many excellent local varieties, and the genetic diversity of cultivated rice has become more and more narrow.
Dongxiang wild rice (O. rufipogon Griff., Hereinafter referred to as DXWR), a Chinese type of wild rice grown in Dongxiang County, Jiangxi Province (28°14′ N latitude and 116°30′ E longitude), is considered to be the northernmost region in the world where O. rufipogon is found, which is one of unique wild resources in Jiangxi province and one of the second class national protected wild plants in China [10]. Previous studies confirmed that DXWR has many useful genes lacking in cultivated rice, including P-deficiency tolerance [11], strong cold tolerance [12], high grain yield [13] and drought resistance [14]. Our previous results and relative reports indicated that DXWR was more resistant to low-P stress than low-P tolerant cultivated rice ‘Dalidao’ and ‘Liantangzao3’, suggesting that DXWR has strong resistance to low-phosphorus stress [15]. Therefore, it is considered to be a valuable resource for the exploitation and utilization of P-deficiency tolerance genes in rice. So far, there has been a great deal of quantitative trait loci (QTLs) for abiotic stress tolerance in DXWR, including P-deficiency tolerance [14, 16,17,18,19,20]. However, the regulatory mechanisms of P-deficiency tolerance have not been fully understood. Thus, a better understanding of P-deficiency tolerance mechanisms would be helpful for breeding P-deficiency tolerance rice cultivars.
In this study, the transcriptome of DXWR under P-deficiency stress was obtained by experiment. Then some candidate genes were found by combining the DEGs interval of this study with the previously identified QTLs interval associated with P-deficiency tolerance. The results will provide a basis for explaining the molecular mechanism of resistance, as well as cloning and utilizing the P-deficiency genes from wild rice.
Methods
Material and treatment
Dongxiang wild rice (O. rufipogon Griff., hereafter referred to as DXWR) was used as test material. Disinfection of rice seeds with 10% sodium hypochlorite, then soak the seeds with water at room temperature for 30 h. Selection of the uniformly germinated seeds of DXWR were grown in a plastic pot in a plant growth chamber at day/night temperature of 30 °C/26 °C (14 h day/10 h night) with relative humidity of 70% and 3000 lx of light intensity. Germinated seeds with coleoptiles 8–12 mm in length, adding Yoshida culture medium [21]. At two and half leaves stage (about 15 days), the seedlings were treated with P-deficiency. P-deficiency (0.016 mM NaH2PO4) and P-sufficiency (0.32 mM NaH2PO4) were as treatment and control, respectively. There were 10 seedlings per treatment with three replications. The growth culture solution was renewed every 3 days. After 9 days treatment under P-deficiency stress, leaves and roots under low P-stressed tissues (LLP and RLP) and control group stressed tissues (LCK and RCK) were collected and immediately frozen in liquid nitrogen, then stored at − 80 °C until RNA extraction.
RNA extraction, cDNA library preparation, and transcriptome sequencing
Total RNA was extracted for three biological replicates from the sampled leaf or root tissues using the TRIzol kit following the instructions (Invitrogen) of manufacturer. The quality and quantity of the resulting RNA were examined using agarose gel electrophoresis and an ND-1000 spectrophotometer (NanoDrop Technology, USA). Magnetic beads with Oligo (dT) were used to isolate mRNA from the total RNA. Mixed with the fragmentation buffer, the mRNA is fragmented into short fragments. Then cDNA was synthesized using the mRNA fragments as templates. Short fragments were purified and resolved with EB buffer for end reparation and poly (A) addition. After that, the short fragments were connected with adapters. After agarose gel electrophoresis, the suitable fragments (200Â bp) were selected for the PCR amplification as templates. The library was sequenced using the Illumina HiSeqâ„¢ 2000 platform.
Reads filtration and assessment of differential gene expression
Before assembly, adaptor sequences, empty reads, low quality sequences with ‘N’ percentage over 10% and those containing more than 50% bases with a Q < 20 were removed using the Perl program written according to the custom method of program editing. After filtering, the remaining reads were called clean reads and used for downstream bioinformatics analysis. The retained high-quality reads were mapped to the Nipponbare reference genome [22] by TopHat. And then the resulting aligned reads were used to create a RABT (Reference Annotation Based Transcript) assembly using Cufflinks [23]. Expression levels for each gene were calculated by quantifying the reads according to the RPKM (reads per kilobase per million reads) method [24]. We used ‘FDR (false discovery rate) ≤ 0.001 and the absolute value of log2 RPKM ratio ≥ 1’ as the threshold to judge the significance of gene expression difference [25].
Gene ontology (GO) term analysis
Blast2GO program was used to classify unigenes to GO terms based on molecular function, biological processes and cellular components [26] for leaf and root tissues, at p < 0.05.
Validation of transcriptome sequencing
qRT-PCR was performed to confirm 15 randomly selected differentially expressed genes (DEGs) (10 up-regulated and 5 down-regulated genes) from each group of P-deficiency-induced genes identified from RNA sequencing using the SYBR premix Ex Taq kit on a StepOnePlusâ„¢ Real-Time PCR System. Diluted cDNA was amplified using gene specific primers (Additional file 1: Table S1) and SYBR Green real-time PCR master mix (Toyobo). All reactions were performed using one biological sample and three technical replicates, and each sample was conducted at least in triplicate and normalized using OsActin1 as an internal control [27].
Results and discussion
Transcriptome sequencing statistics
To understand the molecular mechanism of DXWR under P-deficiency stress, we investigated the gene expression of DXWR in response to P-deficiency stress by transcriptome sequencing. The total RNA was extracted from the leaves of DXWR at seedling stage to construct cDNA library for transcriptome sequencing analysis. In total, 46.78, 48.82, 43.59, and 46.76 million high-quality paired end reads were generated by Illumina-sequencing the LCK, LLP, RCK, and RLP cDNA libraries, respectively (Table 1). The O. sativa ssp. japonica cv. Nipponbare genome has been completely sequenced through Sanger sequencing technology which is considered to be the best tool for assembling and annotating the rice genomes [28,29,30]. Therefore, in this study, we used the Nipponbare genome as reference for reads matching. The alignment results indicated that 73.37–75.20% (71.43–73.72% uniquely matched) of the total reads were mapped to the Nipponbare reference genome and 61.35–63.52% (38.08–39.57% uniquely matched) were mapped to the gene regions (Tables 1, 2). Meanwhile, there was a significant difference in the percentage of reads matched to the genome and the gene region, especially the unique matching reads number, which was similar to previous findings. The previous study reported that about 72% of the total reads mapped to the genome and 46% to the gene regions with 68 and 38% uniquely matched for deep transcriptome sequencing of rhizome and aerial-shoot in Sorghum propinquum [31]. This may be due to the fact that reads match the intergenic spacer or alternatively spliced region of mRNA. The 24.80–26.63% of reads remained unmapped, mainly attributable to gene intervals and the differences between DXWR and reference genome sequence (Table 1). Among more than 60% of the mapped genes, on average, at least 50% were covered by the uniquely mapped reads, and only approximately 15% of the genes had gene coverage of 20% or lower (Additional file 2: Figure S1), which suggested that a high quality of transcriptome data was obtained.
Previous studies have shown that Asian cultivated rice was domesticated from wild varieties (O. rufipogon) [32]. Moreover, the sequence of wild rice W1943 has a great similarity with that of Nipponbare [33]. However, some of the sequences in the W1943 cDNAs that could not be matched to the genome may be located in the gap of genomic sequences or may be related to the W1943 specific gene [33], which suggests that the use of the Nipponbare genome to match wild rice (O. rufipogon) is limited.
Analysis of differentially expressed genes
The raw data obtained from Illumina sequencing can be used to assess the level of gene expression [34]. Putative DEGs from the P-deficiency-stressed and control samples (LLP vs. LCK and RLP vs. RCK) were identified. In the LLP sample, 323 and 861 transcripts were up- (Additional file 3: Table S2) and down-regulated (Additional file 4: Table S3), respectively, when compared to the LCK sample. In the RLP sample, 756 and 230 transcripts were up- (Additional file 5: Table S4) and down-regulated (Additional file 6: Table S5), respectively, when compared to the RCK sample. Among these DEGs, 43 and 35 transcripts were up- (Additional file 7: Table S6) and down-regulated (Additional file 8: Table S7), respectively, in both the LLP and RLP samples, when compared to the control samples (LCK and RCK). 38 transcripts were up-regulated in the RLP sample but down-regulated in the LLP sample (Additional file 9: Table S8), and only 2 transcripts (LOC_Os07g08680.1, a putative protein and LOC_Os11g37950.1, a putative WIP3-Wound-induced protein precursor) were down-regulated in the RLP sample but up-regulated in the LLP sample. The overall distribution of the number of up and down regulated transcripts in roots and leaves can be seen in Fig. 1.
We detected that the expression of several genes were dramatically changed (the absolute value of log2 RPKM ratio > 3), including LOC_Os01g02460.1, LOC_Os02g48850.1, LOC_Os01g27590.1 and LOC_Os03g05334.1 (up regulated in both LLP and RLP), among them, LOC_Os01g02460.1, LOC_Os02g48850.1 may be related to the synthesis of the receptor kinase LRK10 and a specific domain of tigr01615 family proteins; LOC_Os01g72360.1, LOC_Os01g72370.3, LOC_Os02g43410.1, LOC_Os07g27030.1, LOC_Os10g11889.1, LOC_Os10g11889.2 and LOC_Os11g15624.1 (down regulated in both LLP and RLP), among them, LOC_Os01g72370.3 may be related to the protein binding domain of helix loop helix DNA; LOC_Os01g06900.1 and LOC_Os05g25650.1 (up regulated in RLP but down regulated in LLP), the former function is not clear so far, and the latter may be related to the synthesis of Verticillium wilt resistance protein VE.
Furthermore, many of the detected DEGs represented genes that have been previously identified and demonstrated to play roles in responses to P-deficiency stress in cultivated rice (Table 3).
Transcription factors (TFs) have been described as important regulators of transcription and have been reported to play essential roles in abiotic stress responses by regulating a large spectrum of downstream stress responsive genes. In previous studies, OsPHR2 (LOC_Os07g25710) [35], OsUPS (LOC_Os03g13740) [36] and OsPI1 (LOC_Os05g34940) [37] have been reported to be involved in response to P-deficiency stress. There were also some P-deficiency tolerance related genes had been found in DEGs, including OsSPX1 (LOC_Os06g40120) [38], OsWRKY74 (LOC_Os09g16510) [39], OsPupK20-2 (LOC_Os12g26380) [40], OsPupK05 (LOC_Os12g26290) [40], OsPT6 (LOC_Os08g45000) [41], OsBBX2 (LOC_Os02g07930) [42], OsBBX7 (LOC_Os02g49230) [42], OsRCI2-9 (LOC_Os06g44220) [43] and OsbHLH172 (LOC_Os06g12210) [44]. Among them, only OsSPX1 was up-regulated both in leaves and roots of P-deficiency treat, whereas, OsPupK20-2, OsPupK05, OsPT6, OsPT2 and OsBBX1 up-regulated only in roots; OsRCI2-9 and OsbHLH172 up-regulated only in leaves; OsWRKY74 and OsBBX2 down regulated in leaves; and OsBBX7 down regulated in roots.
In this study, we found that some genes have been shown to be associated with P-deficiency tolerance, such as MYB family TF (LOC_Os01g16810.1), WRKY TFs (LOC_Os01g60640.1, LOC_Os01g61080.1, LOC_Os01g53040.1), NAC TFs (LOC_Os01g66120.1 and LOC_Os03g60080.1), rice bHLH TF (LOC_Os01g72370.3), transcription elongation factor (LOC_Os02g04160.1), homeodomain leucine zipper TF (LOC_Os03g10210.1), NAC transcriptional activator (LOC_Os03g21060.1), heat shock TF (LOC_Os03g25120.1), Ethylene response TF (LOC_Os04g46220.1) and AP2 domain containing protein (LOC_Os04g52090.1). This study of gene promoter analysis showed that there were PHR1 binding sites in response to P stress [45]. PHR1, a MYB TF, involved in response to P-deficiency stress [46], which was not very sensitive to P-deficiency stress, and acted a role in the downstream of the phosphorus signal transduction pathway. PHR1 consists of a MYB domain and a coiled coil domain (perhaps to form a dimer with an imperfect palindromic sequence at specific promoters). It plays a role in maintaining the balance of P under the condition of adequate nutrition. It was also found that the increase of sulfur transporter and iron transporter was regulated by P stress [45]. This suggested that the transcription factors played an important role in the transcriptional regulation of downstream genes in plants under P-deficiency stress. Previous studies have found some zinc finger proteins responded to P-deficiency stress, for example, two C2H2- type zinc finger protein gene ZOS3-12 (LOC_Os03g32230) and ZOS5-08 (LOC_Os05g37190): the former was proved to be related to nitrogen stress in rice, and the latter was related to defoliation, which suggested that zinc finger proteins played an important role in the regulation of multiple stress tolerance.
Uptake, transport and translocation of P in plants are performed by P transporters (Pht). Under P-deficiency stress, lipid changed from phospholipid to galactose and sulfanilamide, phosphatase activity increased, and phosphate monoester decomposed into Pi and the related fatty acid and phosphate transporter gene will be induced [47]. In this study, we also found some transporter genes that have been shown to be associated with P-deficiency tolerance. In the LLP sample, we found three genes editing transmembrane amino acid transporter (LOC_Os01g41420.1, LOC_Os06g12320.1, LOC_Os01g61044.2), two genes editing inorganic phosphate transporter (LOC_Os03g05620.1 and LOC_Os06g21950.1) and eight editing transporter protein (LOC_Os01g50820.1, LOC_Os08g06010.1, LOC_Os01g36720.1, LOC_Os08g31670.1, LOC_Os11g42430.1, LOC_Os02g36450.1, LOC_Os03g03680.1 and LOC_Os01g07310.1). In the RLP sample, we found six up-regulated genes editing inorganic phosphate transporter (LOC_Os06g21950.1, LOC_Os08g45000.1, LOC_Os10g30770.1, LOC_Os03g05640.1, LOC_Os10g30790.1 and LOC_Os06g21920.1), among which, LOC_Os06g21950.1 was also up-regulated in the LLP sample. So far, there are 4 Pht families, namely Pht1, Pht2, Pht3 and Pht4. Pht1 gene family encodes a high affinity protein (Phosphate transporter) to regulate the rate of P crossing the plasma membrane. In the Arabidopsis genome, there are 9 members in the Pht1 gene family, namely Pht1;1 ~ Pht1;9, four members (Pht1;1, Pht1;2, Pht1;3 and Pht1;4) expressed promoters, which were induced by P-deficiency stress. What’s more, Pht1;1 ~ Pht1;4 were related to plant uptake of Pi from soil. The Pht1 gene family plays an important role in the process of P transport in plants, and express in flowers, cotyledons, pollen, leaf vascular tissues and shoots [48, 49]. Pht2 and Pht3 family encode Pi transporter associated with organelle. Pht2;1 is the only member of the Pht2 family and is a low affinity phosphate transporter, which is present in the chloroplast membrane [50, 51]. Pht3 family has 3 members, and the protein is located in the mitochondria. The Pht4 gene family consists of 6 members, and its structure is similar to that of P transporter SLC17/type1. The Pht4 gene expressed in both roots and leaves, 5 of which were present in plastids, and the other in Golgi, suggesting that the Pht4 gene family is associated with the transport of Pi in the cytosol, plastids and Golgi [52]. Although the function of a low affinity Pi transporter system is present in the root of the plant, the genes encoding these transporters remain to be identified.
The transcriptome changes of cultivated rice to P starvation have been reported in previous study [53]. In this study, we compared the DEGs between DXWR and cultivated rice. We found two genes encoding metallothioneins (LOC_Os12g38270.2 and LOC_Os12g38290.1) were up-redulated in roots, which was same with previous research [53]. Metallothioneins affect metal tolerance and homeostasis and scavenge reactive oxygen species [54], which could be a mechanism to overcome the increase in certain ion concentration, such as iron, upon Pi starvation. Under Pi starvation, plants overaccumulated some ions, including iron [45, 55, 56]. We also found one gene encoded iron transporter (LOC_Os09g23300.1) was up-regulated both in roots and leaves and two genes encoded vacuolar iron transporters (LOC_Os04g45520.1 and LOC_Os09g23300.1) were up-regulated in roots and leaves, respectively (LOC_Os04g45520.1 was up-regulated in roots and LOC_Os09g23300.1 was up-regulated in leaves). In addition, LOC_Os08g06010.1, a putative glycerol 3-phosphate permease, and LOC_Os03g40670.1, a putative glycerophosphoryl diester phosphodiesterase suggested to be involved in Pi remobilization [57], were also up-regulated both in roots and leaves after P-deficiency stress. One gene encoded MYB family TF (LOC_Os02g22020.1), one encoded transporter (LOC_Os08g31670.1) and one encoded glycerophosphoryl diester phosphodiesterase family protein (LOC_Os02g31030.1) which were up-regulated both in roots and leaves of DXWR but only up-regulated in leaves of cultivated rice.
To confirm the validity of the DEG data, quantitative RT-PCR was performed to investigate the expression patterns of 15 randomly selected genes under the same conditions. Expression trends were consistent for all transcripts in RNA-Seq analysis and quantitative RT-PCR analysis, with a correlation coefficient of (R2 = 0.9195) (Fig. 2). Thus, the DEGs detected in this study can be considered to be a high accuracy.
Functional classification by gene ontology
The GO gene ontology mapping software (WEGO) was used to classify the function and draw gene ontology tree, and the down and up regulated transcripts in roots and leaves were classified into the functional groups. In the present study, a total of 12,180 root transcripts in RLP vs. RCK and 6803 leaf transcripts in LLP vs. LCK were assigned GO terms. Among the 12,180 root transcripts, 1829 were annotated for their molecular function, 5865 transcripts were annotated for their cellular component, and 4486 were annotated for their biological process (Fig. 3). Among the 6803 leaf transcripts, 1241 were annotated for their molecular function, 2816 transcripts were annotated for their cellular component, and 2746 were annotated for their biological process (Fig. 4). In the classification of biological processes, cellular processes and metabolic processes were the most functional groups, which indicated that the DXWR had a wide range of metabolic activities under the P-deficiency stress. In the classification of cell components, the transcripts of cells, cell components and organelles were the most abundant. In the classification of molecular function, the transcription of immobilized and active were the most highly represented groups.
We also identified biological process, cellular component and molecular function GO terms that were over-represented (p < 0.05) among the DEGs of LLP vs. LCK and RLP vs. RCK, respectively (Additional file 10: Table S9, Additional file 11: Table S10). It was found that the expression of transcripts in the extracellular matrix increased in roots and leaves for the cellular component category (GO: 0005576, annotated as extracellular region), suggesting that the same cell components in the root and leaf are involved in the P-deficiency response.
Kyoto encyclopedia of genes and genomes (KEGG) pathway mapping
KEGG pathway analysis showed that 700 of the 1184 leaf DEGs and 622 of the 986 root DEGs could be classified into 20 functional categories and 114 and 117 subcategories, respectively. Furthermore, the over-represented KEGG Orthology (KO) terms (Q < 0.05) could be classified into 7 and 4 categories, respectively (Fig. 5). As shown in Fig. 5, the most common KO terms represented by both the leaf and root DEGs was the translation KEGG pathways. The over-represented KO terms for the leaf and root DEGs were further classified into 16 and 17 subcategories, respectively (Additional file 12: Table S11, Additional file 13: Table S12). Among these subcategories, 2 subcategories were over-represented among the leaf DEGs, which are RNA transporters and mRNA monitoring pathways, suggesting that these pathway may regulate the expression of P-deficiency inducible genes. Leaf and root did not have the same subcategory. Moreover, 12 KO terms were exclusively enriched among the leaf DEGs and 7 KO terms were exclusively enriched among the root DEGs. This finding suggests that there could be considerable differences in the biochemical and physiological processes involved in the P-deficiency responses of leaves and roots, and these annotations provide a valuable resource for investigating the specific processes, functions, and pathways involved in such differences.
DEGs mapped to previously identified P-deficiency responses related QTL intervals
Map based cloning is a method for identifying potential QTLs sites. Many QTLs have been identified in DXWR were associated with yield and multiple resistance. Based on the Gramene QTL database, a total of 57 QTLs related to P-deficiency stress have been found. We located 278 genes differentially regulated by P-deficiency stress on 10 of these identified QTL intervals in rice (Additional file 14: Table S13). Among them, the QTL Accession AQBD003, AQCI008, AQCI013, and AQCI011 had the greatest number of co-localized DEGs with 86, 33, 31 and 29 genes, respectively (Additional file 15: Table S14).
86 DEGs were co-localized within the AQBD003 interval, including two genes encoding PHD-finger family protein (LOC_Os12g06480.1 and LOC_Os12g06480.2), four encoding NB-ARC domain protein (LOC_Os12g31200.1, LOC_Os12g17410.1, LOC_Os12g10330.1 and LOC_Os12g10410.1), and one gene encoding helix-loop-helix DNA-binding domain containing protein (LOC_Os12g32400.1). Previous research has found that the NB-ARC domain of OsPDRH9Â N protein response to light regulation, biological and abiotic stress [58].
33 DEGs were co-localized within the AQCI008 interval, including one gene encoding AP2 TF (LOC_Os03g12950.1), one encoding DHHC zinc finger protein domain (LOC_Os03g16790.1), and two genes encoding Ribosomal protein gene (LOC_Os03g17580.1 and LOC_Os03g15870.1). In addition, this interval also included some energy related proteins, such as fructose-1,6-bisphosphatase, glyoxal oxidase, oxidoreductase and protein phosphatase 2C. Protein phosphatase 2C (PP2C) plays an important role in biological signal transduction. PP2C participates in stress response in higher plants. As a negative regulator of most signaling pathways, PP2C can directly bind to kinase or regulatory proteins, as well as directly bind to DNA regulating the expression of related genes [59].
31 DEGs were co-localized within the AQCI013 interval, including three genes encoding NB-ARC domain protein (LOC_Os12g17410.1, LOC_Os12g10330.1 and LOC_Os12g10410.1). 29 DEGs were co-localized within the AQCI011 interval, including one gene encoding AP2 TF (LOC_Os04g42570.1) and another one encoding ethylene-responsive TF (LOC_Os04g46220.1). This interval also included one YABBY domains protein gene (LOC_Os04g45330.1). Study has shown that YABBY protein plays an important role in plant abiotic stress [60].
Comparative analysis of the DEGs and previously identified P-deficiency-related QTL intervals in DXWR
Until now, a total of 23 P-deficiency-related QTLs on chromosomes 1, 2, 3, 7, 8, 9 and 11 were identified in DXWR [61]. Among the 23 QTLs, the QTL cluster (the region RM254-RM1233) on chromosome 11 was high significant for P-deficiency tolerant (LOD 3.31 and LOD 3.52 for qFWS-4 and qFWR-4, respectively) [61].
In the present study, 13 DEGs were co-localized within the qFWS-4 interval (Additional file 16: Table S15), including one gene encoding bZIP TF (LOC_Os11g05640.1), one encoding MYB-like DNA-binding domain protein (LOC_Os11g03440.1), one encoding calmodulin depedent protein kinases (LOC_Os11g02240.1), one encoding auxin-repressed protein (LOC_Os11g44810.2), regulating transcription of auxin-responsive genes [62]. And another one gene encoded homeobox domain containing protein (LOC_Os11g06020.1), whose products are transcription factors sharing a characteristic protein fold structure that bound DNA [63,64,65], regulating gene transcription. Furthermore, another gene, LOC_Os11g44310.1, which encoded calmodulin binding protein, regulated the expression of calmodulin to regulate the intracellular calcium concentration to control cells in many important biochemical reactions, which might play an important role in stress responses [66].
11 DEGs were co-localized within the qFWR-4 interval (Additional file 17: Table S16), including one gene encoding phosphoglycerate mutase (LOC_Os11g05260.1), which was the key enzyme in sugar metabolism process, regulating the adaptation of plants to environment. Another one encoded nucleoside-triphosphatase (LOC_Os11g03290.1) [67], which maintained the connection between the nucleus cytoplasm and cytoplasm, including transcription and regulating biosynthesis to adapt to the environment. The interval also included LOC_Os11g05400.1, which encoded Ser/Thr protein phosphatase family protein, modified other proteins by chemically adding phosphate groups to them and, thus, regulated cellular pathways, signal transduction, and responses to biotic and abiotic stresses [68]. Meanwhile, LOC_Os11g04300.1 and LOC_Os11g03940.1 encoded retrotransposon proteins, whose replicative mode of transposition by means of an RNA intermediate rapidly increased the copy numbers of elements and thereby could increase genome size, which could induce mutations by inserting near or within genes to response to biotic and abiotic stresses. LOC_Os11g04290.1 and LOC_Os11g05380.1, which encoded cytochrome P450, modified other proteins by chemically adding OH group to them or as an enzyme catalyzed reaction and [69], thus, regulated cellular pathways, signal transduction, and responses to biotic and abiotic stresses.
Conclusion
DXWR has a lot of useful agronomic traits we need, such as resistance of cold, drought, salt and P-deficiency. Therefore, it is considered to be an important resource for rice breeding. In this study, we analyzed the transcriptome of roots and leaves of DXWR seedlings under P-deficiency stress. A large number of DEGs and some critical paths were identified, such as RNA transport and mRNA monitoring path. By combining the DEGs identified in the present study with previously identified P-deficiency resistance QTLs from rice, important candidate genes were identified, including a variety of transcription factors and some functional protein genes. These findings will be useful in the future studies of molecular adaptations to P-deficiency stress and will facilitate the genetic manipulation of rice to improve its P-deficiency resistance.
Abbreviations
- P:
-
phosphorus
- DXWR:
-
Dongxiang wild rice
- DEGs:
-
differentially expressed genes
- FDR:
-
false discovery rate
- GO:
-
gene ontology
- KEGG:
-
Kyoto Encyclopedia of Genes and Genomes
- KO:
-
KEGG Orthology
- RPKM:
-
reads per kilobase per million reads
- QTLs:
-
quantitative trait loci
- TFs:
-
transcription factors
- qRT-PCR:
-
quantitative real-time PCR
- NAC:
-
NAM, ATAF, and CUC
- PP:
-
protein phosphatase
- AP2:
-
APETALA2
- bHLH:
-
basic helix-loop-helix
- WEGO:
-
Web Gene Ontology Annotation Plot
- bZIP:
-
basic leucine zipper
References
Abel S, Ticconi CA, Delatorre CA. Phosphate sensing in higher plants. Physiol Plant. 2002;115(1):1–8.
Zhu DF, Cheng SH, Zhang YP, Lin XQ, Chen HZ. Analysis of status and constraints of rice production in the world. Sci Agric Sin. 2010;43(3):474–9.
Wu ZH, He LY, Zuo XD, Yang JF, Men YY. Influence of low-P stress on the biological characteristics of different rice genotypes at different stages. J Plant Nutr Fertilizer. 2008;14(2):227–34.
Wang RC, Cheng SH, Cao LY. Advancements in phosphorus deficiency tolerance study in rice (oryza sativa L.). Chin Agric Sci Bull. 2009;25(6):77–83.
Li ZJ, Xie Y, Dai AY, Liu LF, Li ZC. Root and shoot traits responses to phosphorus deficiency and QTL analysis at seedling stage using introgression lines of rice. J Genet Genomics. 2009;36(3):173–83.
Li Y, Gu M, Zhang X, Zhang J, Fan H, Li P, Li Z, Xu G. Engineering a sensitive visual-tracking reporter system for, real-time monitoring phosphorus deficiency in tobacco. Plant Biotechnol J. 2014;12(6):674–84.
Cordell D, Drangert JO, White S. The story of phosphorus: global food security and food for thought. Glob Environ Change. 2009;19(2):292–305.
Wissuwa M, Ae N. Genotypic variation for tolerance to phosphorus deficiency in rice and the potential for its exploitation in rice improvement. Plant Breeding. 2001;120(1):43–8.
Sun CQ, Wang XK, Li ZC, Yoshimura A, Iwata N. Comparison of the genetic diversity of common wild rice (Oryza rufipogon Griff.) and cultivated rice (O. sativa L.) using RFLP markers. Theor Appl Genet. 2001;102:157–62.
Xie J, Agrama HA, Kong D, Zhuang J, Hu B, Wan Y, Yan W. Genetic diversity associated with conservation of endangered Dongxiang wild rice (Oryza rufipogon). Genet Resour Crop Evol. 2010;57(4):597–609.
Chen XR, Chen M, He HH, Zhu CL, Peng XS, He XP, Fu JR, Ouyang LJ. Low-phosphorus tolerance and related physiological mechanism of Xieqingzao B//Xieqingzao B/Dongxiang wild rice BC1F9 populations. Chin J Appl Ecol. 2011;22(5):1169–74.
Luo XD, Dai LF, Cao JF, Jian SR, Chen YL, Hu BL, Xie JK. Identification and molecular cytology analysis of cold tolerance introgression lines derived from Oryza sativa L. mating with O. rufipogon Griff. Euphytica. 2012;187(3):461–9.
Tan L, Liu F, Xue W, Wang GJ, Ye S, Zhu ZF, Fu YC, Wang XK, Sun CQ. Development of Oryza rufipogon and O. sativa Introgression Lines and Assessment for Yield-related Quantitative Trait Loci. J Integr Plant Biol. 2007;49(6):871–84.
Zhang FT, Cui FL, Zhang LX, Wen XF, Luo XD, Zhou Y, Li X, Wan Y, Zhang JE, Xie JK. Development and identification of a introgression line with strong drought resistance at seedling stage derived from Oryza sativa L. mating with Oryza rufipogon Griff. Euphytica. 2014;200(1):1–7.
Chen XR, Chen M, He HH, Zhu CL, Peng XS, He XP, Fu JR, Ouyang LJ. Identification of low-phosphorus-tolerance in Dongxiang wild rice (Oryza rufipogon Griff.). Acta Agriculturae Universitatis Jiangxiensis. 2011;33(3):0405–11.
Zhang X, Zhou S, Fu Y, Su Z, Wang X, Sun C. Identification of a drought tolerant introgression line derived from Dongxiang common wild rice (O. rufipogon Griff.). Plant Mol Biol. 2006;62(1–2):247–59.
Xia H, Lu BR, Su J, Chen R, Rong J, Song Z, Wang F. Normal expression of insect-resistant transgene in progeny of common wild rice crossed with genetically modified rice: its implication in ecological biosafety assessment. Theor Appl Genet. 2009;119(4):635–44.
Huang DR, Chen J, Lai FX, Liu GJ, Zhuang JY. Analysis of quantitative trait loci for resistance to brown planthopper in Dongxiang wild rice (Oryza rufipogon Griff.). Acta Agronomica Sinica. 2012;38(2):210–4.
Mao D, Yu L, Chen D, Li L, Zhu Y, Xiao Y, Zhang D, Chen C. Multiple cold resistance loci confer the high cold tolerance adaptation of Dongxiang wild rice (Oryza rufipogon) to its high-latitude habitat. Theor Appl Genet. 2015;128(7):1359–71.
Xiao N, Huang WN, Li AH, Gao Y, Li YH, Pan CH, Ji H, Zhang XX, Dai Y, Dai ZY, et al. Fine mapping of the qLOP2 and qPSR2-1 loci associated with chilling stress tolerance of wild rice seedlings. Theor Appl Genet. 2015;128(1):173–85.
Yoshida S, Forno DA, Cock JH. Laboratory manual for physiological studies of rice. Int Rice Res Inst. 1976.
Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12(1):323.
Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, Baren MJ, Salzberg SL, Wold BJ, Pachter L. Transcript assembly and abundance estimation from RNA-Seq reveals thousands of new transcripts and switching among isoforms. Nat Biotechnol. 2010;28(5):511.
Mortazavi A, Williams BA, Mccue K, Schaeffer L, Wold B. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat Methods. 2008;5(7):621–8.
Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001;29(4):1165–88.
Conesa A, Götz S. Blast2GO: a comprehensive suite for functional analysis in plant genomics. Int J Plant Genomics. 2008;2008:619832.
Zhang FT, Luo LD, Hu BL, Wan Y, Xie JK. YGL138(t), encoding a putative signal recognition particle 54 kDa protein, is involved in chloroplast development of rice. Rice. 2013;6(1):1–10.
Ouyang S, Liu J, Jones KM, Gansberger K, Moffat K, Hill J, Bera J, Fadrosh D, Jin SH, Johri S, et al. The map-based sequence of the rice genome. Nature. 2005;436:793–800.
Kawahara Y, de la Bastide M, Hamilton JP, Kanamori H, McCombie WR, Ouyang S, Schwartz DC, Tanaka T, Wu J, Zhou S, et al. Improvement of the Oryza sativa Nipponbare reference genome using next generation sequence and optical map data. Rice. 2013;6(1):4.
Zhang FT, Xu T, Mao LY, Yan SY, Chen XW, Wu ZF, Chen R, Luo XD, Xie JK, Gao S. Genome-wide analysis of Dongxiang wild rice (Oryza rufipogon Griff.) to investigate lost/acquired genes during rice domestication. BMC Plant Biol. 2016;16(1):1–11.
Zhang T, Zhao X, Wang W, Huang L, Liu X, Zong Y, Zhu L, Yang D, Fu B, Li Z. Deep transcriptome sequencing of rhizome and aerial-shoot in Sorghum propinquum. Plant Mol Biol. 2014;84(3):315–27.
Huang XH, Kurata N, Wei XH, Wang ZX, Wang AH, Zhao Q, Zhao Y, Liu KY, Lu HY, Li WJ, et al. A map of rice genome variation reveals the origin of cultivated rice. Nature. 2012;490(7421):497–501.
Lu T, Yu S, Fan D, Mu J, Shangguan Y, Wang Z, Minobe Y, Lin Z, Han B. Collection and comparative analysis of 1888 full-length cDNAs from wild rice Oryza rufipogon Griff. W1943. DNA Res. 2008;15(5):285–95.
Bloom JS, Khan Z, Kruglyak L, Singh M, Caudy AA. Measuring differential gene expression by short read sequencing: quantitative comparison to 2-channel gene expression microarrays. BMC Genomics. 2009;10(1):221.
Zhou J, Jiao FC, Wu ZC, Li YY, Wang XM, He XW, Zhong WQ, Wu P. OsPHR2 is involved in phosphate-starvation signaling and excessive phosphate accumulation in shoots of plants. Plant Physiol. 2008;146:1673–86.
Hur YJ, Yi YB, Lee JH, Chung YS, Jung HW, Yun DJ, Kim KM, Park DS, Kim DH. Molecular cloning and characterization of OsUPS, a U-box containing E3 ligase gene that respond to phosphate starvation in rice (Oryza sativa). Mol Biol Rep. 2012;39(5):5883–8.
Hou XL, Wu P, Jiao FC, Jia QJ, Chen HM, Yu J, Song XW, Yi K. Regulation of the expression of OsIPS1, and OsIPS2, in rice via systemic and local Pi signalling and hormones. Plant Cell Environ. 2010;28(3):353–64.
Wang C, Ying SH, Huang H, Li K, Wu P, Shou H. Involvement of OsSPX1 in phosphate homeostasis in rice. Plant J Cell Mol Biol. 2010;57(5):895–904.
Dai X, Wang Y, Zhang WH. OsWRKY74, a WRKY transcription factor, modulates tolerance to phosphate starvation in rice. J Exp Bot. 2016;67(3):947–60.
Chin JH, Gamuyao R, Dalid C, Bustamam M, Prasetiyono J, Moeljopawiro S, Wissuwa M, Heuer S. Developing rice with high yield under phosphorus deficiency: Pup1 sequence to application. Plant Physiol. 2011;156(3):1202–16.
Seo HM, Jung Y, Song S, Kim Y, Kwon T, Kim DH, Jeung SJ, Yi YB, Yi G, Nam MH, et al. Increased expression of OsPT1, a high-affinity phosphate transporter, enhances phosphate acquisition in rice. Biotechnol Lett. 2008;30(10):1833–8.
Huang J, Zhao X, Weng X, Wang L, Xie W. The rice B-Box zinc finger gene family: genomic identification, characterization, expression profiling and diurnal analysis. PLoS ONE. 2012;7(10):e48242.
Ren HY, Wei Q. Isolation and functional analysis of phosphate-responsive gene OsRCI2-9 in oryza sativa. Sci Agric Sin. 2015;48:831–40.
Yu F, Park S, Rodermel SR. Functional redundancy of AtFtsH metalloproteases in thylakoid membrane complexes. Plant Physiol. 2005;138(4):1957–66.
Misson J, Raghothama KG, Jain A, Jouhet J, Block MA, Bligny R, Ortet P, Creff A, Somerville S, Rolland N, et al. A genome-wide transcriptional analysis using Arabidopsis thaliana Affymetrix gene chips determined plant responses to phosphate deprivation. Proc Natl Acad Sci USA. 2005;102(33):11934–9.
Rubio V, Linhares F, Solano R, Martin AC, Iglesias J, Leyva A, Paz-Ares J. A conserved MYB transcription factor involved in phosphate starvation signaling both in vascular plants and in unicellular algae. Genes Dev. 2001;15(16):2122–33.
Ticconi CA, Abel S. Short on phosphate: plant surveillance and countermeasures. Trends Plant Sci. 2004;9(11):548–55.
Karthikeyan AS, Varadarajan DK, Mukatira UT, D’Urzo MP, Damsz B, Raghothama KG. Regulated expression of Arabidopsis phosphate transporters. Plant Physiol. 2002;130(1):221–33.
Mudge SR, Rae AL, Diatloff E, Smith FW. Expression analysis suggests novel roles for members of the Pht1 family of phosphate transporters in Arabidopsis. Plant J. 2002;31(3):341.
Daram P, Brunner S, Rausch C, Steiner C, Amrhein N, Bucher M. Pht2;1 encodes a low-affinity phosphate transporter from Arabidopsis. Plant Cell. 1999;11(11):2153.
Versaw WK, Harrison MJ. A chloroplast phosphate transporter, PHT2;1, influences allocation of phosphate within the plant and phosphate-starvation responses. Plant Cell. 2002;14(8):1751–66.
Guo B, Jin Y, Wussler C, Blancaflor EB, Motes CM, Versaw WK. Functional analysis of the Arabidopsis PHT4 family of intracellular phosphate transporters. New Phytol. 2008;177(4):889–98.
Secco D, Jabnoune M, Walker H, Shou H, Wu P, Poirier Y, Whelan J. Spatio-temporal transcript profiling of rice roots and shoots in response to phosphate starvation and recovery. Plant Cell. 2013;25(11):4285.
Cobbett C, Goldsbrough P. Phytochelatins and metallothioneins: roles in heavy metal detoxification and homeostasis. Annu Rev Plant Biol. 2002;53(1):159.
Hirsch J, Marin E, Floriani M, Chiarenza S, Richaud P, Nussaume L, Thibaud MC. Phosphate deficiency promotes modification of iron distribution in Arabidopsis plants. Biochimie. 2006;88(11):1767–71.
Lei M, Liu Y, Zhang B, Zhao Y, Wang X, Zhou Y, Raghothama KG, Liu D. Genetic and genomic evidence that sucrose is a global regulator of plant responses to phosphate starvation in Arabidopsis. Plant Physiol. 2011;156(3):1116–30.
Cheng Y, Zhou W, El Sheery NI, Peters C, Li M, Wang X, Huang J. Characterization of the Arabidopsis glycerophosphodiester phosphodiesterase (GDPD) family reveals a role of the plastid-localized AtGDPD1 in maintaining cellular phosphate homeostasis under phosphate starvation. Plant J. 2011;66(5):781–95.
Zhao TT. Functional identification and expression analysis in stress responses of a NB-ARC domain containing protein. M.S. dissertation. Xiamen Univ. 2012.
Hu XB, Song FM, Zheng Z. Structure and function of protein phosphatase 2C in higher plants. Chin J Cell Biol. 2005;27(1):29–34.
Li XB, Yang CC, Qiu NW. Bioinformatic analysis of YABBY protein family in Arabidopsis and Chinese cabbage. Shandong Agric Sci. 2012;44(12):1–6.
Luo XD, Liu J, Dai LF, Zhang FT, Wan Y, Xie JK. Linkage map construction and QTL identification of P-deficiency tolerance in Oryza rufipogon Griff. at early seedling stage. Euphytica. 2017;213(4):96.
Lee J, Han CT, Hur Y. Molecular characterization of the Brassica rapa auxin-repressed, superfamily genes, BrARP1 and BrDRM1. Mol Biol Rep. 2013;40(1):197–209.
Gehring WJ. The homeobox in perspective. Trends Biochem Sci. 1992;17(8):277.
Gehring WJ. Exploring the homeobox. Gene. 1993;135(1–2):215–21.
Bürglin TR, Affolter M. Homeodomain proteins: an update. Chromosoma. 2016;125(3):497.
Yong B, Wang XY, Xu P, Zheng HY, Fei XT, Hong ZX, Ma QQ, Miao YZ, Yuan XH, Jiang YS, et al. Isolation and abiotic stress resistance analyses of a catalase gene from Ipomoea batatas (L.) Lam. Biomed Res Int. 2017;2017:6847532.
Allen TD, Cronshaw JM, Bagley S, Kiseleva E, Goldberg MW. The nuclear pore complex: mediator of translocation between nucleus and cytoplasm. J Cell Sci. 2000;113(10):1651–9.
Chen YL, Hung MH, Chu PY, Chao TI, Tsai MH, Chen LJ, Hsiao YJ, Shih CT, Hsieh FS, Chen KF. Protein phosphatase 5 promotes hepatocarcinogenesis through interaction with AMP-activated protein kinase. Biochem Pharmacol. 2017;138:49–60.
Dang TT, Franke J, Tatsis EC, O’Connor SE. Dual catalytic activity of a cytochrome P450 controls bifurcation at a metabolic branch point of alkaloid biosynthesis in Rauwolfia serpentina. Angew Chem Int Ed Engl. 2017;56(32):9440–4.
Authors' contributions
LXD and XJK conceived and designed the experiments; DQW, LXD, CYL and ZFT performed the experiments; DQW, LXD, CYL and ZY analyzed the data; DQW and LXD wrote the manuscript; ZY and ZFT reviewed the manuscript. All authors read and approved the final manuscript.
Acknowledgements
Not applicable.
Competing interests
The authors declare that they have no competing interests.
Availability of data and materials
The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.
Consent for publication
Not applicable.
Ethics approval and consent to participate
Not applicable.
Funding
This research was partially supported by the National Natural Science Foundation of China (31260255, 31360147 and 31660384), the Natural Science Foundation of Jiangxi Province, China (20151BAB204008 and 20171BAB204021), and the Major Projects in Jiangxi Province (20161ACF60022).
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Author information
Authors and Affiliations
Corresponding authors
Additional files
Additional file 1: Table S1.
Primers used for real-time PCR in this study.
Additional file 2: Figure S1.
Distribution of genes coverage in the leaves and roots of Dongxiang wild rice (DXWR) seedlings with or without low phosphorus treatment, respectively. A leaves without low phosphorus treatment (LCK). B leaves with low phosphorus treatment (LLP). C roots without low phosphorus treatment (RCK). D roots with low phosphorus treatment (RLP). Gene coverage is the percentage of a gene that is covered by reads and defined as the ratio of the number of bases in a gene covered by uniquely mapped reads to the number of total bases in the gene. The pie graph demonstrates the detailed percentage of the different gene coverage listing on the left of the pie graph.
Additional file 3: Table S2.
List of up-regulated genes in LLP vs. LCK.
Additional file 4: Table S3.
List of down-regulated genes in LLP vs. LCK.
Additional file 5: Table S4.
List of up-regulated genes in RLP vs. RCK.
Additional file 6: Table S5.
List of down-regulated genes in RLP vs. RCK.
Additional file 7: Table S6.
List of up-regulated genes both in the LLP vs. LCK and RLP vs. RCK.
Additional file 8: Table S7.
List of down-regulated genes both in the LLP vs. LCK and RLP vs. RCK.
Additional file 9: Table S8.
List of down-regulated genes in the LLP vs. LCK but up-regulated in the RLP vs. RCK.
Additional file 10: Table S9.
Significant GO terms of DEGs in the biological process, cellular component and molecular function category for LLP vs. LCK.
Additional file 11: Table S10.
Significant GO terms of DEGs in the biological process, cellular component and molecular function category for RLP vs. RCK.
Additional file 12: Table S11.
Significant KO terms of DEGs in the LLP vs. LCK (Q-value < 0.05).
Additional file 13: Table S12.
Significant KO terms of DEGs in the RLP vs. RCK (Q-value < 0.05).
Additional file 14: Table S13.
Previously identified P-deficiency responses related QTL intervals.
Additional file 15: Table S14.
DEGs mapped to previously identified P-deficiency responses related QTL intervals.
Additional file 16: Table S15.
Co-localized DEGs within the qFWS-4 interval.
Additional file 17: Table S16.
Co-localized DEGs within the qFWR-4 interval.
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.
About this article
Cite this article
Deng, QW., Luo, XD., Chen, YL. et al. Transcriptome analysis of phosphorus stress responsiveness in the seedlings of Dongxiang wild rice (Oryza rufipogon Griff.). Biol Res 51, 7 (2018). https://doi.org/10.1186/s40659-018-0155-x
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s40659-018-0155-x