Insights into genes encoding respiratory burst oxidase homologs (RBOHs) in rubber tree (Hevea brasiliensis Muell. Arg.)

Industrial Crops & Products 128 (2019) 126–139

Contents lists available at ScienceDirect

Industrial Crops & Products journal homepage: www.elsevier.com/locate/indcrop

Insights into genes encoding respiratory burst oxidase homologs (RBOHs) in rubber tree (Hevea brasiliensis Muell. Arg.)

T

⁎

Zhi Zou , Jianghua Yang, Xicai Zhang Danzhou Investigation & Experiment Station of Tropical Crops, Ministry of Agriculture and Rural Affairs, Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences, Haikou, 571101, Hainan, P. R. China

A R T I C LE I N FO

A B S T R A C T

Keywords: Hevea brasiliensis Whole-genome duplication Tapping panel dryness Leaf senescence Respiratory burst oxidase Gene expression profiling

Superoxide-producing NADPH oxidases (NOXs), also known as respiratory burst oxidase homologs (RBOHs) in plants, are key players in the signaling network of reactive oxygen species (ROS) and involved in various plant processes. Despite the importance of RBOHs, little information is available in Hevea brasiliensis, a rubber-producing tree of the Euphorbiaceae family. This study presents the first genome-wide analysis of Rboh family genes in this special species. By using arabidopsis and rice RBOH proteins as queries, homology search of the rubber tree genome resulted in a relatively high number of ten Rboh genes. According to phylogenetic analysis, these genes were divided into seven groups named RBOHD, RBOHC, RBOHB, RBOHE, RBOHF, RBOHH, and RBOHN, which is consistent with best-reciprocal-hit-based sequence comparison with representative plant species. The former six groups are widely distributed in core eudicots, whereas RBOHN represents a novel but ancient group that has lost in many plants including poplar, grapevine, arabidopsis, and rice. Comparative genomics analysis showed that all rubber tree paralogs were derived from the recent whole-genome duplication (WGD) shared by cassava. Gene expression profiling revealed diverse patterns of different family members over various rubber tree tissues as observed in cassava, which also supports expression divergence between duplicates. Moreover, key members, which are related to tapping panel dryness, leaf senescence, hormone, and/or stress responses, were also identified in rubber tree. Results obtained from this study will facilitate further studies of Rboh genes in rubber tree as well as other species.

1. Introduction Hevea brasiliensis Muell. Arg. (2n = 4x = 36), commonly known as rubber tree or rubber, is a perennial tropical plant of the family Euphorbiaceae (spurge), which also includes several economically important species such as cassava (Manihot esculenta Crantz, 2n = 4x = 36), castor bean (Ricinus communis L., 2n = 20), and physic nut (Jatropha curcas L., 2n = 22) (Bredeson et al., 2016; Zou et al., 2016, 2018). Like cassava (Olsen and Schaal, 2001), rubber tree is also indigenous to the Amazon rainforest in South America (Lieberei, 2007; Zou et al., 2009). Recent whole genome sequencing and comparative genomics analyses revealed that rubber tree and cassava shared one recent whole-genome doubling event after the divergence with castor bean and physic nut, though the morphology of rubber tree and cassava is distinct and cassava is characterized as a perennial shrub (Zou et al.,

2015; Bredeson et al., 2016; Pootakham et al., 2017). According to available genome sequences, the rubber tree genome is relatively bigger (Reyan7-33-97, 1.37 Gb) than that of cassava (AM560-2, 582.28 Mb), and also contains more protein-coding genes (i.e. 43,792 vs 33,033) (Bredeson et al., 2016; Tang et al., 2016). In the past 100 years, rubber tree has been widely planted in Southeast Asia for commercial production of natural rubber (NR, cis-1,4-polyisoprene), an indispensable industrial raw material (Zou et al., 2009). NR is specifically synthesized and stored in the laticifer, a highly differentiated single-cell-type tissue that is mainly located in the soft inner bark of the tree trunk (Hao and Wu, 2000). The NR-containing latex which represents the laticifer cytoplasm is usually harvested through tapping the bark every 2–5 days (Zou et al., 2015). Ethephon (ETH), an ethylene generator, is also exploited for yield promotion (Liu, 2016). Over the past decades, the rubber yield has significantly increased for wide cultivation of high-

Abbreviations: BRH, best reciprocal hit; Chr, chromosome; DUOX, dual oxidase; EST, expressed sequence tag; ETH, ethephon; Ka, nonsynonymous substitution rate; Ks, synonymous substitution rate; Mya, million years ago; NOX, NADPH oxidase; NR, natural rubber; OG, orthologous group; Rboh, respiratory burst oxidase homolog; ROS, reactive oxygen species; RNA-seq, RNA sequencing; TPD, tapping panel dryness; WGD, whole-genome duplication ⁎ Corresponding author. E-mail addresses: [email protected], [email protected] (Z. Zou), [email protected] (J. Yang), [email protected] (X. Zhang). https://doi.org/10.1016/j.indcrop.2018.11.005 Received 4 August 2018; Received in revised form 1 November 2018; Accepted 2 November 2018 0926-6690/ © 2018 Elsevier B.V. All rights reserved.

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

2. Materials and methods

yielding clones and extensive utilization of ETH. However, great losses have also been caused by the frequent occurrence of TPD (tapping panel dryness), a complex physiological syndrome characterized as tapping incision blocked partly or entirely during latex exploiting (Chen et al., 2003; Zou et al., 2012, 2017a). Although a great deal of effort has been made, the mechanism underlying remains poorly understood (Venkatachalam et al., 2007; Li et al., 2010, 2016; Liu et al., 2015; Putranto et al., 2015; Montoro et al., 2018). The involvement of ROS (reactive oxygen species), mainly resulted from high-strength tapping and ETH over-stimulating, was proposed in 1980s (Chrestin et al., 1984). Thus far, several genes related to ROS-scavenging systems have been characterized, including superoxide dismutase (SOD), peroxidase (POD), catalase (CAT), ascorbate peroxidase (APX), dehydroascorbate reductase (DHAR), glutathione peroxidase (GPX), and glutathione reductase (GR) (Chao et al., 2015; Deng et al., 2015). By contrast, little is known about ROS generators. ROS are highly reactive molecules that are able to damage cellular components but also act as important signaling molecules (Suzuki et al., 2011; Dietz et al., 2016). Although there are various ROS generating systems, the superoxide-producing NADPH oxidase (NOX), first identified in human phagocytic cells, is most extensively studied (Kaur et al., 2014). This enzyme, also known as respiratory burst oxidase for the rapid increase in oxygen consumption during phagocytosis, is comprised of several subunits including the catalytic core gp91phox or NOX2. NOX2 is a membrane-integrated glycoprotein with a molecular mass of about 91 kDa, which contains NADPH and FAD-binding domains in the C-terminal cytoplasmic region, and heme-binding His pairs in the N-terminal transmembrane region (Sumimoto, 2008). NOX homologs are widely distributed in a broad range of organisms including vertebrates, invertebrates, fungi, and higher plants (Bedard et al., 2007). In animals, besides NOX2, six additional types have also been found, which are called NOX1, NOX3, NOX4, NOX5, DUOX1 (dual oxidase 1), and DUOX2. In addition to possessing the basic NOX structure, at the N-terminus, NOX5 also contains four calcium-binding EF-hand motifs while DUOX types contain two EF-hand motifs and a peroxidase-like domain (Sumimoto, 2008). NOXs in plants are known as respiratory burst oxidase homologs (RBOHs) which are the closest homologs of NOX5 (Sagi and Fluhr, 2006). Since the first plant homolog OsRBOHA was characterized in rice (Oryza sativa) (Groom et al., 1996), a high number of RBOHs have been found in several plant species (Lightfoot et al., 2008; Marino et al., 2011; Cheng et al., 2013; Wang et al., 2013a, Chang et al., 2016). They were proven to play key roles in various plant processes such as growth, development, hormone signaling, responses to abiotic stresses and biotic interactions (Torres et al., 2002; Marino et al., 2012; Kaur et al., 2014; Wang et al., 2018). Their transcript level was reported to be regulated by various hormones and environmental stresses (Sagi and Fluhr, 2006; Wang et al., 2013a; Chang et al., 2016), and the enzyme activity is tightly controlled by various post-translational modifications such as dimerization, calcium binding, phosphatidic acid binding, phosphorylation, and nitrosylation (Ogasawara et al., 2008; Zhang et al., 2009; Yun et al., 2011). The regulatory mechanism of EF-hand motifs has been investigated by crystallizing the N-terminal region of OsRBOHB (Oda et al., 2010). More recently, crystal structures of two key components of the basic NOX structure, i.e. the catalytic FAD/NADPH-binding region and hemebinding transmembrane region were also resolved by using Cylindrospermum stagnale NOX5, which was more likely acquired by cyanobacteria through gene transfer from a higher eukaryote (Magnani et al., 2017), highlighting key structural elements common to the entire NOX family. Despite the importance of RBOHs, little information is available in rubber tree. This paper presents a genome-wide identification of Rboh family genes in rubber tree, their evolution patterns, the global gene expression profiles, and also their putative roles in TPD, leaf senescence, hormone and stress responses. These findings will facilitate further studies of Rboh genes in rubber tree and other species.

2.1. Identification and manual curation of Rboh family genes Protein sequences of arabidopsis (Arabidopsis thaliana) and rice Rboh family genes described before (see Supplementary Table S1) were used as queries to search for rubber tree homologs from an in-house Reyan7-33-97 genome (Tang et al., 2016). Sequences with an E-value of less than 1e-5 in the TBLASTN search (Altschul et al., 1997) were collected, and positive genomic sequences were predicted as previously described (Zou et al., 2015). Subsequently, all gene models were further validated with rubber tree cDNAs, Sanger ESTs (expressed sequence tags) and RNA-seq (RNA sequencing) reads when available. The presence of the conserved NADPH_Ox domain (PF08414) was confirmed using SMART (http://smart.embl-heidelberg.de/). To facilitate analysis of the evolution pattern of HbRboh gene family, a similar approach was also used to identify homologs from several representative genomes such as cassava, castor bean, physic nut, poplar (Populus trichocarpa), and grapevine (Vitis vinifera), which were accessed from Phytozome v12 (https://phytozome.jgi.doe.gov/pz/portal.html) or NCBI (http://www. ncbi.nlm.nih.gov/, last accessed April 2018). They were selected for the following reasons: cassava, castor bean, and physic nut represent Euphorbiaceae plants in Malpighiales, where cassava shared the recent whole-genome duplication (WGD, named ρ in this study) with rubber tree and the latter two did not experience additional WGD after the socalled whole-genome triplication γ event shared by all core eudicots; poplar, a Salicaceae plant in Malpighiales, also experienced one lineage-specific WGD; and, grapevine represents an outgroup without any WGD (Tuskan et al., 2006; Jaillon et al., 2007; Chan et al., 2010; Wu et al., 2015; Bredeson et al., 2016). 2.2. Synteny analysis and definition of orthologous groups Synteny analysis was performed as described before (Zou et al., 2017b), where duplicate pairs were identified using the all-to-all BLASTP method and syntenic blocks were inferred using MCScanX (Wang et al., 2012). Duplicates derived from WGDs were defined when duplicated genes are located in syntenic blocks of duplicated chromosomes (Chrs). Tandem duplications were considered when two duplicated genes were consecutive in a genome. Additionally, BRH (best reciprocal hit)-based BLASTP analysis was used to define orthologs across different species, and the systematic name was assigned based on the best arabidopsis ortholog and phylogenetic analysis (see below). 2.3. Sequence alignment and phylogenetic analysis Multiple sequence alignments of deduced RBOH proteins were performed using MUSCLE (Edgar, 2004). Unrooted phylogenetic trees based on alignments were constructed using MEGA 6.0 (Tamura et al., 2013) with the maximum likelihood method (Jones-Taylor-Thornton (JTT) model, uniform rates and bootstrap of 1,000 replicates). Gene structures and sequence alignment were displayed using GSDS (http:// gsds.cbi.pku.edu.cn/) or Boxshade (http://www.ch.embnet.org/ software/ BOX_form.html), respectively. Ks (synonymous substitution rate) and Ka (nonsynonymous substitution rate) were calculated using codeml of the PAML package (Yang, 2007). 2.4. Analysis of protein properties and conserved motifs Theoretical molecular weight (MW), isoelectric point (pI), and grand average of hydropathicity (GRAVY) were calculated using ProtParam (http://web.expasy.org/protparam/). Protein subcellular localization was predicted using Cell-PLoc (Chou and Shen, 2010), and typical domains such as NADPH_Ox, EF-hand motif, heme-binding transmembrane helices, Ferric_reduct, FAD_binding_8, and NAD_binding_6 were inferred from SMART as well as sequence alignment. 127

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Interestingly, both MeRbohB and MeRbohN are located on Chromosome 8, which is similar to JcRbohB and JcRbohN located on Chromosome 6. By contrast, MeRbohB and MeRbohD are located on two different chromosomes, though JcRbohD stands very close to JcRbohB on Chromosome 6 (Fig. 1 and Supplementary Fig. S1). These results suggest possible chromosome rearrangement after their divergence. Macrosynteny was observed between rubber tree and cassava, which allows mapping rubber tree scaffolds to the corresponding cassava chromosomes based on synteny analysis. As a result, eight HbRboh genes could be anchored to seven cassava chromosomes, only excluding HbRbohD2 and HbRbohF2 (Fig. 1). The CDS of five gene pairs identified in rubber tree and cassava were shown to exhibit relatively high identity, varying from 87.6 to 92.2% (Table 2). The reciprocal BLASTP analysis revealed that HbRbohC1, HbRbohC2, HbRbohD1/HbRbohD2, HbRbohF1, and HbRbohF2 are the orthologs of MeRbohC1, MeRbohC2, MeRbohD, MeRbohF1, and MeRbohF2, respectively. MeRbohC1 and MeRbohC2 are more likely to result from the ρ WGD for their location in syntenic blocks. Exhibiting a similar Ks value to that of MeRbohC1/MeRbohC2, MeRbohF1, and MeRbohF2 may also be derived from the ρ WGD, though the chromosome location of MeRbohF2 has not been resolved. In fact, microsynteny was observed between HbRbohF1 and HbRbohF2-coding scaffolds, which is similar to HbRbohC1/HbRbohC2 and HbRbohD1/ HbRbohD2. Thereby, these gene pairs were also defined as WGD duplicates. Interestingly, cassava has completely lost the ortholog of HbRbohD2 after its divergence with rubber tree. Despite sharing the same recent WGD, the Ks value of duplicate pairs in rubber tree is relatively smaller than that in cassava, supporting a lower rate of evolution in perennial big trees. In fact, similar Ks values were also observed in another big tree, poplar. The Ka/Ks ratio of rubber tree and cassava duplicates varies from 0.0714 to 0.1998, which is similar to poplar duplicates (from 0.1424 to 0.2772) (Table 2).

Other conserved motifs were investigated using MEME (Bailey et al., 2009) with the optimized parameters: any number of repetitions; maximum number of motifs, 25; and, the optimum width of each motif, between 6 and 200 residues. The MAST program (http://meme-suite. org/tools/mast) was also used to search detected motifs in protein databases. 2.5. Gene expression analysis based on RNA-seq Detailed information of transcriptome data used in this study was shown in Supplementary Table S2, and all of them were generated via Illumina pair-end sequencing. Quality control was performed using fastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/), and filtered reads were mapped to coding sequences (CDS) of Rboh genes and other protein-coding genes using Bowtie 2 (Langmead and Salzberg, 2012). The relative transcript level was presented using FPKM (fragments per kilobase of exon per million fragments mapped) (Mortazavi et al., 2008). Parameters “FDR < 0.001” (Audic and Claverie, 1997) and “log2Ratio ≥ 1” were used as the threshold to judge the significance of gene expression difference. 2.6. Gene expression analysis based on qRT-PCR To learn more about the role of ethylene on the expression of HbRbohs in leaf, Reyan7-33-97 derived in vitro plantlets were used for the following experiments. The one-year-old bagged plantlets were grown in a greenhouse (Rubber Research Institute, Chinese Academy of Tropical Agricultural Sciences, Danzhou, China) illuminated with natural sunlight, and the chlorophyll content of leaves was monitored as described before (Zou et al., 2017b). For ETH treatment, six batches of at least three plantlets with similar growth performance were adopted. Mature leaves were treated with 50 μM ETH solution as described before (Zou et al., 2017a) and samples (at least three leaves per plantlet) were collected at six time points (i.e. 0, 2, 4, 6, 8, and 72 h) after the treatment. At the last time point, only leaves with the chlorophyll content of 45–55% relative to mature leaves were selected, which are defined as mid-senescent leaves (Zou et al., 2017a). Total RNA isolation, cDNA synthesis, and qRT-PCR analysis were performed as previously described (Zou et al., 2015, 2017a). Primers used in this study can be found in Supplementary Table S3, where HbYLS8 was used as the reference gene. All qRT-PCR assays were performed in triplicate for each biological sample. Estimation of relative gene abundance and statistical analysis were carried out as described before (Zou et al., 2015).

3.2. Phylogenetic analysis and definition of orthologous groups As shown in Fig. 2, an unrooted phylogenetic tree was constructed using full-length proteins of ten HbRbohs, nine MeRbohs, seven JcRbohs, seven RcRbohs, ten PtRbohs, seven VvRbohs, ten AtRbohs, and nine OsRbohs. According to the tree, seven phylogenetic groups named RBOHD, RBOHC, RBOHB, RBOHE, RBOHF, RBOHH, and RBOHN, were identified. The former six which also exist in arabidopsis were named based on their best arabidopsis counterparts, whereas RBOHN represents a novel group that is absent from rice, arabidopsis, grapevine as well as poplar. The result is highly consistent with the BRH-based homologous analysis, which revealed seven OGs (orthologous groups) across eight species examined (Table 3). RBOHD/OG1, RBOHB/OG3, RBOHE/OG4, RBOHF/OG5, and RBOHH/OG6 are widely distributed; RBOHC/OG2 exists in all examined eudicots but excludes from rice, a model monocotyledon plant; whereas, RBOHN/OG7 is limited to four Euphorbiaceae plant species examined (Fig. 2 and Table 3).

3. Results 3.1. Identification, chromosome location and synteny analysis of Rboh family genes The initiative TBLASTN search of the rubber tree genome resulted in 18 putative loci from 17 scaffolds. After discarding loci putatively encoding FROs (ferric chelate reductases) without the conserved NADPH_Ox domain, a total of ten RBOH-coding loci were retained and they are located on ten different scaffolds (Table 1). Although only two genes (i.e. HbRbohF2 and HbRbohB) were shown to have EST hits in NCBI GenBank, the expression of other family members was supported by available RNA-seq reads derived from transcriptomes of somatic embryogenesis, shoot apex, root, leaf, bark, laticifer, flower, and seed. Survey of cassava, castor bean, physic nut, poplar, and grapevine genomes resulted in nine, seven, seven, ten, or seven Rboh family members, respectively. Their expression was all supported by available ESTs and/or RNA-seq reads, which also allows the optimization of nine predicted gene models (Supplementary Table S1). In cassava, except for MeRbohF2, other family members were shown to distribute across seven out of the 18 cassava chromosomes (Fig. 1).

3.3. Exon-intron structures, sequence features and conserved motifs The exon-intron structure of Rboh family genes was analyzed on the basis of optimized gene models. Results showed that 69 examined genes contain at least seven introns (i.e. AtRbohD), where OsRbohD harbors the maximum of 14 introns (Table 1 and Supplementary Table S1). Generally, the exon-intron structure is highly conserved within a certain group, supporting the subclassification: RBOHD, RBOHC, and RBOHB usually contain 11 introns, whereas RBOHE, RBOHF, RBOHH, and RBOHN feature 13 introns (Fig. 3a). Nevertheless, gain or loss of introns was observed in several Rboh genes: JcRbohD has lost the fifth intron as found in its orthologs in rubber tree, cassava, castor bean, poplar, and grapevine, i.e. HbRbohD1/HbRbohD2, MeRbohD, RcRbohD, PtRbohD1/PtRbohD2, and VvRbohD1/VvRbohD2; AtRbohD has lost the fifth, eighth, ninth, and tenth introns as found in its orthologs in other 128

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Table 1 Rboh family genes in rubber tree. Gene name

HbRbohD1 HbRbohD2 HbRbohC1 HbRbohC2 HbRbohB HbRbohE HbRbohF1 HbRbohF2 HbRbohN HbRbohH

Scaffold position

scaffold0854:388025-383531 scaffold0294:998357-993858 scaffold0702:191765-197756 scaffold0169:1471360-1464635 scaffold0170:1829314-1834285 scaffold0272:200130-194199 scaffold0280:1267116-1259368 scaffold0319:27915-36806 scaffold0478:284762-7000 scaffold1323:156233-151848

Nucleotide length (bp, from start to stop codons) CDS

Gene

2733 2760 2742 2736 2652 2778 2850 2853 2205 2661

3888 4085 5202 5713 4142 5438 6842 7563 5836 4243

Intron no.

EST hits

AA

MW (kDa)

pI

GRAVY

11 11 11 11 11 13 13 13 13 13

0 0 0 0 1 0 0 2 0 0

910 919 913 911 883 925 949 950 734 886

102.62 103.51 103.74 102.65 100.35 104.84 107.79 107.96 84.64 101.08

9.13 9.23 9.04 8.98 9.15 8.91 9.28 9.27 9.24 9.37

−0.271 −0.283 −0.301 −0.221 −0.286 −0.210 −0.287 −0.238 −0.117 −0.178

Abbreviations: amino acid, AA; base pair, bp; coding sequence, CDS; expressed sequence tag, EST; grand average of hydropathicity, GRAVY; isoelectric point, pI; kilodalton, kDa; molecular weight, MW.

species; AtRbohC and AtRbohG have lost the fifth intron as found in its paralog AtRbohA and orthologs in rubber tree, cassava, physic nut, castor bean, poplar, and grapevine, i.e. HbRbohC1/HbRbohC2, MeRbohC1/MeRbohC2, JcRbohC, RcRbohC, PtRbohC1/PtRbohC2, and VvRbohC; JcRbohB was shown to have lost the first intron as found in its orthologs in rubber tree, castor bean, poplar, grapevine, and arabidopsis, i.e. HbRbohB, RcRbohB, PtRbohB, VvRbohB, and AtRbohB; MeRbohB has lost the fifth intron as found in its orthologs in other species; OsRbohB and OsRbohH have lost the eleventh intron as found in their orthologs in other species; OsRbohA has lost the first intron as found in its paralog OsRbohC and orthologs in rubber tree, cassava, physic nut, castor bean, poplar, grapevine, and arabidopsis, i.e. HbRbohF1/ HbRbohF2, MeRbohF1/MeRbohF2, JcRbohF, RcRbohF, PtRbohF1/

PtRbohF2, VvRbohF, and AtRbohF/AtRbohI; AtRbohH and AtRbohJ have lost the third intron as found in its orthologs in rubber tree, cassava, physic nut, castor bean, poplar, grapevine, and arabidopsis, i.e. HbRbohH, MeRbohH, JcRbohH, RcRbohH, PtRbohH1/PtRbohH2, and VvRbohH; OsRbohD has lost the fourth intron as found in its paralog OsRbohE and orthologs in other species; OsRbohD has gain one more intron following the fifth exon as found in its orthologs in other species; RcRbohN has lost the second intron as found in its orthologs in rubber tree, cassava, and physic nut, i.e. HbRbohN, MeRbohN, and RcRbohN; MeRbohN also has lost the seventh intron as found in its orthologs in rubber tree, castor bean, and physic nut (Fig. 3a). The deduced HbRBOH proteins consist of 734–950 amino acid residues, and the theoretical MW and pI value were predicted to range

Fig. 1. Chromosomal locations of nine MeRboh genes and their collinear genes in rubber tree. Shown are seven cassava chromosomes encoding MeRbohs and the chromosome serial number is indicated at the top of each chromosome. The line connects the duplicate pair located in syntenic blocks. Eight HbRbohs shown just behind their collinear genes in cassava were marked in orange. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 129

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

from 84.64 to 107.96 kDa or from 8.91 to 9.37, respectively. The average protein length and MW are similar across eight examined species, i.e., 898 residues or 101.92 kDa in rubber tree, 897 residues or 101.90 kDa in cassava, 887 residues or 100.72 kDa in physic nut, 887 residues or 100.86 kDa in castor bean, 907 residues or 103.01 kDa in poplar, 897 residues or 101.82 kDa in grapevine, 906 residues or 102.93 kDa in arabidopsis and 903 residues or 101.55 kDa in rice. The average pI value is approximately 9.16 in rubber tree and most examined species, however, grapevine RBOHs harbor a relatively small value of 8.80. The GRAVY value was all shown to be less than 0 (varying from -0.081 to -0.362), indicating their hydrophilic feature. Since four Euphorbiaceae plant species contain all seven RBOH groups, a focused analysis of their conserved motifs was performed by using MEME. Among 25 motifs identified, Motifs 1–10, 12, 13, 15, 17, 21, and 22 are broadly distributed; Motifs 11 and 14 are present in most

Table 2 Rboh duplicates derived from recent WGD in rubber tree, cassava and poplar. Ks and Ka were calculated using PAML. Duplicate pair

Identity (%)

Ks

Ka/Ks

HbRbohD1/HbRbohD2 HbRbohC1/HbRbohC2 HbRbohF1/HbRbohF2 MeRbohC1/MeRbohC2 MeRBohF1/MeRBohF2 PtRbohC1/PtRbohC2 PtRbohD1/PtRbohD2 PtRbohF1/PtRbohF2 PtRbohH1/PtRbohH2

91.2 92.2 91.5 87.6 89.4 91.2 91.0 92.4 89.9

0.2732 0.2071 0.2472 0.3810 0.3491 0.2466 0.2504 0.2363 0.2314

0.0714 0.1998 0.1334 0.1629 0.1265 0.1629 0.2090 0.1424 0.2772

Abbreviations: nonsynonymous substitution rate, Ka; synonymous substitution rate, Ks.

Fig. 2. Phylogenetic analysis of RBOHs in rubber tree, cassava, castor bean, physic nut, poplar, grapevine, arabidopsis, and rice. Sequence alignment was performed using MUSCLE and the unrooted phylogenetic tree was constructed using bootstrap maximum likelihood tree (1,000 replicates) method of MEGA6. The distance scale denotes the number of amino acid substitutions per site. The name of each group/OG is indicated next to the corresponding group. (Abbreviations: orthologous group, OG). 130

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Table 3 Species-specific distribution of seven OGs identified in the Rboh gene family. OGs that are limited to rubber tree and cassava were shown in bold, and systematic names were assigned to groups only when at least one member was found in poplar, grapevine or arabidopsis. Group

OG

Rubber

Cassava

Physic nut

Castor bean

Poplar

Grapevine

Arabidopsis

Rice

RBOHD

OG1

HbRbohD1 HbRbohD2

MeRbohD

JcRbohD

RcRbohD

PtRbohD1 PtRbohD2

VvRbohD1 VvRbohD2

OsRbohI

RBOHC

OG2

RcRbohC

JcRbohB

RcRbohB

PtRbohC1 PtRbohC2 PtRbohB

VvRbohC

OG3

MeRbohC1 MeRbohC2 MeRbohB

JcRbohC

RBOHB

HbRbohC1 HbRbohC2 HbRbohB

AtRbohD AtRbohA AtRbohG AtRbohC

VvRbohB

AtRbohB

RBOHE

OG4

HbRbohE

MeRbohE

JcRbohE

RcRbohE

PtRbohE

VvRbohE

AtRbohE

RBOHF

OG5

RcRbohF

JcRbohN JcRbohH

RcRbohN RcRbohH

PtRbohF1 PtRbohF2 ND PtRbohH1 PtRbohH2

VvRbohF

OG7 OG6

MeRbohF1 MeRbohF2 MeRbohN MeRbohH

JcRbohF

RBOHN RBOHH

HbRbohF1 HbRbohF2 HbRbohN HbRbohH

AtRbohF AtRbohI ND AtRbohH AtRbohJ

OsRbohB OsRbohH OsRbohF OsRbohG OsRbohA OsRbohC ND OsRbohE OsRbohD

ND VvRbohH

ND

Abbreviations: not detected, ND; orthologous group, OG.

binding and Ca2+-dependent phosphorylation; Motifs 3, 8, 9, and 4 belong to the Ferric_reduct domain (PF01794) which is characterized as a transmembrane component similar to ferric reductase; Motifs 13, 1, 16, and 6 belong to the FAD_binding_8 domain (PF08022) bearing the FAD-binding site; whereas, Motifs 10, 2, and 7 belong to the NAD_binding_6 domain (PF08030) bearing the NADPH-binding site (Fig. 4 and Supplementary Fig. S2). Compared with the C-terminal, the Nterminal region of RBOHs is relatively variable and RBOHN was shown

groups with the exception of RBOHD; Motif 19 is present in most groups except for RBOHH and RBOHN; Motif 24 is present in most groups except for RBOHE and RBOHN; Motifs 16, 18, and 20 are limited to RBOHD, RBOHC, RBOHB, and RBOHF; Motifs 23 and 25 are limited to RBOHB, RBOHC, and RBOHD (Fig. 3b and Supplementary Fig. S2). Motifs 15, 21, 11, 5, and 14 belong to the NADPH_Ox domain which catalyzes the production of superoxides; Motifs 14 and 17 are characterized as the EF-hand motif (PF13499), which is involved in Ca2+-

Fig. 3. Structural analysis of rubber tree, cassava, castor bean, and physic nut Rboh genes. a. Shown is the graphic representation of exon-intron structures displayed using GSDS. b. Shown is the distribution of conserved motifs among RBOH proteins, where different motifs are represented by different color blocks as indicated at the bottom of the figure and the same color block in different proteins indicates a certain motif. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article). 131

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Fig. 4. Sequence alignment of HbRBOH proteins. Identical and similar amino acids are highlighted in black or dark grey, respectively. Histidine residues involved in heme binding are indicated by arrows. Conserved domains such as NADPH_Ox, EF-hand motif, Ferric_reduct, FAD_binding_8 and NAD_binding_6 are indicated by single underlines, whereas putative transmembrane α-helices (TMs) are indicated by double underlines. (Abbreviations: transmembrane, TM) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

132

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Fig. 5. Tissue-specific expression profiles of rubber tree and cassava Rboh genes. Color scale represents FPKM normalized log10 transformed counts where green indicates low expression and red indicates high expression. (Abbreviations: friable embryogenic callus, FEC; shoot apical meristem, SAM; somatic organized embryogenic structure, OES; root apical meristem, RAM) (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article).

transcript level was not very high (Fig. 5a). According to their expression patterns over various tissues, ten HbRbohs were grouped into three main clusters: Cluster I is preferentially expressed in bark and female flower, including two RbohFs (i.e. HbRbohF1 and HbRbohF2), two RbohCs (i.e. HbRbohC1 and HbRbohC2) and one RbohD (i.e. HbRbohD1), where HbRbohD1 and HbRbohF2 were also highly expressed in male flower; Cluster II is typically expressed in root, including HbRbohB and HbRbohE; Cluster III is predominantly expressed in male flower, including HbRbohN, HbRbohD2, and HbRbohH, where HbRbohN was also highly expressed in leaf (Fig. 5a). To learn more about the expression evolution of HbRboh genes, tissue-specific expression profiling of Rboh family genes was also performed in cassava, which shared the ρ WGD with rubber. The examined 11 tissues, i.e. shoot apical meristem (SAM), lateral bud, leaf, mid-vein, petiole, stem, fibrous root, storage root, root apical meristem (RAM), friable embryogenic callus (FEC) and somatic organized embryogenic structure (OES), were previously generated by Wilson et al. (2017) with at least two replicates. Interestingly, the expression of MeRbohH was only detected in FEC, and the transcript level is extremely low. Both MeRbohC1 and MeRbohF2 were constitutively expressed, where MeRbohF2 represents the most expressed gene in FEC, SAM and stem and MeRbohC1 represents the most expressed gene in other examined tissues. According to the total transcripts of the gene family, MeRbohs were most abundant in RAM, stem, fibrous root, and petiole (Class I); moderate in mid-vein, storage root, lateral bud, and leaf (Class II, counting 34–58% of Class I); and, considerably low in FEC, OES, and SAM (Class III, counting 16–25% of Class I). Based on tissue-specific expression patterns, MeRbohs can be divided into three main clusters: Cluster I is typically expressed in fibrous root, including MeRbohC2, MeRbohF1, MeRbohE, MeRbohD, and MeRbohN; Cluster II is predominantly expressed in stem and petiole, including MeRbohC1 and MeRbohF2; and, Cluster III that includes MeRbohB is preferentially expressed in FEC and OES (Fig. 5b).

to lack a large part of the N-terminal present in other groups (Fig. 4).

3.4. Tissue-specific expression profiles of Rboh genes in rubber tree and cassava Along with whole genome sequencing, several tissue transcriptomes of Reyan7-33-97 have also been sequenced by using the Illumina platform as described before (Tang et al., 2016), i.e. root, leaf, bark, laticifer, male flower, female flower, and seed. As shown in Fig. 5a, transcriptional profiling showed that ten HbRbohs were expressed in at least one of the examined tissues, i.e., ten in male flower, nine in female flower, nine in seed, nine in root, nine in leaf, nine in bark, and two in laticifer. Based on the FPKM value, the transcripts of the total gene family were shown to be most abundant in root and bark (denoted Class I for convenience, same as follows); moderate in female flower, male flower, and leaf (Class II, counting 17–44% of Class I); and, considerably low in seed and laticifer (Class III, counting 1–3% of Class I). Members of RBOHC and RBOHB contributed the major transcripts in most examined tissues: HbRbohC1 and HbRbohC2 occupied 71%, 43%, or 36% of total transcripts in female flower, leaf and bark respectively, and HbRbohB occupied 87%, 83% or 34% in root, seed and bark respectively. In laticifer, RBOHF occupied 95% of the total transcripts, whereas RBOHD and RBOHN occupied 60% of the total transcripts in male flower. RBOHN and RBOHD may also play an important role in leaf, counting about 28% or 11% of the total transcripts respectively. By contrast, members in RBOHE and RBOHH were considerably less expressed (Fig. 5a). Several key genes were also identified in a certain tissue and expression divergence was observed between duplicates. HbRbohB represents the most expressed gene in root and seed, but it was not or lowly expressed in other examined tissues. HbRbohC1 represents the most expressed gene in bark, leaf, and female flower, in contrast, the transcript level of its paralog HbRbohC2 was relatively low in all examined tissues. Similar expression pattern was also observed between other two duplicates, where the transcript levels of HbRbohD1 and HbRbohF2 were relatively higher than that of their paralogs respectively. HbRbohF2 represents the most and second expressed gene in laticifer or female flower, respectively. HbRbohN represents the most and second expressed gene in male flower or leaf, respectively. The expression of HbRbohH was shown to be male flower-specific, though its

3.5. Transcriptional profiling of HbRboh genes during leaf development Expression profiles of HbRbohs during leaf development were further analyzed over five typical stages as described before (Fang et al., 2016; Zou et al., 2017b), i.e. bronze, color-change, pale-green, mature, 133

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Fig. 6. Expression profiles of HbRboh genes during leaf development and in response to various hormones and stresses. a. Profiles during leaf development. b. Responses to ethephon treatment. c. Responses to jasmonic acid or coronatine treatment. d. Responses to cold stress in two varieties with different sensitivity. e. Responses to M. ulei and C. mycelium infection. f. Profiles in bark or laticifer of healthy and TPD-affected rubber trees. (Abbreviations: control, CK; coronatine, COR; jasmonic acid, JA; tapping panel dryness, TPD).

and senescent. As shown in Fig. 6a, the total transcripts were similar during early leaf development, however, dramatic decline was observed when leaves have matured and especially senesced. Based on their

expression patterns, eight HbRbohs were grouped into two main clusters: Cluster I is predominantly expressed in color-change leaves, including two RbohCs (i.e. HbRbohC1 and HbRbohC2), HbRbohD2 and 134

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

mycelium infection (Hurtado Páez et al., 2015) was also investigated. As shown in Fig. 6e, after the treatment of M. ulei for 48 h, HbRbohC2 and HbRbohN were upregulated by 2.0 or 3.3 folds in FX3864 leaves respectively. As for C. mycelium, young emerging leaves of the sensitive RRII105 and tolerant GT1 were inoculated and collected as a mixed sample from 6, 12, and 24 h after the treatment. Compared with the noninoculated control, three or four genes were significantly regulated in RRII105 and GT1 respectively: In RRII105, HbRbohC2, HbRbohD1, and HbRbohF1 were upregulated by 2.4, 3.4 or 5.7 folds respectively; in GT1, HbRbohD1, HbRbohD2, HbRbohF1, and HbRbohN were upregulated by 11.4, 13.5, 5.1 or 3.2 folds respectively, though their transcript levels were not very high (Fig. 6e).

HbRbohN; Cluster II is preferentially expressed in pale-green leaves, including two RbohFs (i.e. HbRbohF1 and HbRbohF2), HbRbohD1 and HbRbohE. Compared with mature leaf, two genes were shown to be significantly regulated in senescent leaf, i.e. about three-folds increase for HbRbohD1 and about five-folds decrease for HbRbohF1 (Fig. 6a). 3.6. The response of HbRboh genes to ethephon treatment in laticifer and bark Given the importance of ETH stimulation on rubber yield promotion, Reyan7-33-97 latex samples collected at four time points (i.e. 0, 3, 12, and 24 h) after ETH treatment were selected for RNA-seq as previously described (Tang et al., 2016). Expression profiling showed that the transcript level of HbRbohs was extremely low in laticifer (Fig. 6b), which is consistent with tissue-specific expression analysis as described above. HbRbohF2, a relatively more expressed gene in this special tissue, was significantly upregulated at the time point of 12 h (about three folds). In PR107 bark (Liu et al., 2016), ETH application for 24 h could upregulate HbRbohF1 by about 2.1 folds, by contrast, transcripts of HbRbohD2 and HbRbohN were shown to significantly decrease at both time points examined (i.e. 37.5 or 19.3 for 8 and 24 h respectively) (Fig. 6b).

3.9. HbRboh genes associated with TPD in bark and laticifer The frequent occurrence of TPD makes it to be one of the most important factors that affect rubber production (Zou et al., 2017a). In China, TPD trees were classed into five main grades based on the severity, where the second and fourth grades representing slight and severe TPD in this study are defined by a quarter or from a half to three quarters of the tapping cut affected respectively (Zou et al., 2012). HbRbohs associated with TPD in PR107 bark were examined based on transcriptomes of healthy and severe TPD trees (without ETH stimulation) (Liu et al., 2016). Interestingly, the profiling revealed that most HbRbohs were downregulated in TPD trees and five of them were identified as differentially expressed genes, i.e. HbRbohB with about 345 folds, HbRbohC1 with about 2.6 folds, HbRbohD1 with about 53.8 folds, HbRbohF1 with about 3.8 folds, and HbRbohF2 with about 4.9 folds (Fig. 6f). However, the effect seems to be tissue-specific. In laticifer, only HbRbohF2 and HbRbohC1 transcripts were detected, where HbRbohC1 was excluded from differential analysis for its low expression. The profiling showed that HbRbohF2 was not significantly regulated, though slight and severe TPD trees were both examined (Fig. 6f). In the case of ETH stimulation (Montoro et al., 2018), about two-folds upregulation of HbRbohF2 was observed in laticifer of ETH-treated trees in related to the control. It’s worth noting that no significant change was observed between severe (ETH-induced) and slight (with and without ETH-treated) TPD trees (Fig. 6f).

3.7. The response of HbRboh genes to jasmonic acid or coronatine treatment in bark Jasmonic acid (JA) is a phytohomorne that can not only induce secondary laticifer differentiation but also promote rubber productivity (Hao and Wu, 2000; Tian et al., 2015; Liu et al., 2018). Transcriptional profiling showed that, after treating Reyan88-13 bark with 0.02% or 0.04% w/w JA for 24 h, three HbRbohs (i.e. HbRbohB, HbRbohC2, and HbRbohF2) were significantly regulated: the transcript level of HbRbohB increased by about 4.2 or 2.2 folds under two concentrations, respectively; HbRbohC2 increased by about two folds and HbRbohF2 decreased by about 4.4 folds under the high concentration (Fig. 6c). Coronatine (COR), an effective mimic of the active form of JA, can also induce the differentiation of secondary laticifers (Zhang et al., 2015; Wu et al., 2016). Expression profiling showed that, after treating Reyan7-33-97 epicormic shoots with 20 μM COR for a period of time, HbRbohC1 and HbRbohN were significantly regulated: HbRbohC1 was upregulated by 2.1 folds in the mixed early sample collected at time points of 1 h, 2 h, and 4 h whereas HbRbohN was downregulated by 15.2 folds in the mixed late sample collected at time points of 8 h, 1 d, 2 d, and 3 d (Fig. 6c).

3.10. HbRboh genes associated with ethylene response in leaf and ethyleneinduced leaf senescence It’s well established that ethylene acts as an important plant hormone regulating fruit ripening, abscission and leaf senescence (Iqbal et al., 2017). Thereby, ETH solution was used to artificially induce leaf senescence in rubber tree, and leaves of mid-senescent as well as early samples after the treatment were subjected for expression profiling of five relatively more expressed genes (i.e. HbRbohC1, HbRbohC2, HbRbohD1, HbRbohF1, and HbRbohF2). As shown in Fig. 7, qRT-PCR analysis revealed that four genes were significantly regulated at least one time point. Among them, HbRbohC2 was initially inhibited by ETH but upregulated in senescent leaves; conversely, HbRbohF2 was initially induced by ETH but downregulated in senescent leaves; the transcripts of HbRbohD1 and HbRbohC1 gradually increased with the extension of time after the treatment, and the difference is that HbRbohC1 ultimately returned to the control level in senescent leaves (Fig. 7).

3.8. The response of HbRboh genes to cold or pathogenic infection in leaf As a tropical plant, rubber tree is highly susceptible to low temperature and thus selecting plants resistant to cold is an important objective in the breeding of rubber tree. During the past decades, several cold-resistant clones have been developed in China, e.g. CATAS93114. In a previous study, plantlets of CATAS93-114 and the sensitive Reken501 were subjected to 4 °C chilling, and their leaf transcriptome response was characterized at 0, 2, 8, and 24 h after the treatment (Cheng et al., 2018). As shown in Fig. 6d, the profiling revealed that the response pattern of two clones is distinct, where three genes (i.e. HbRbohC1, HbRbohC2, and HbRbohN) in Reken501 were significantly regulated at least one time point in contrast to none in CATAS93-114: the transcripts of HbRbohC1 increased by 2.6 or 2.5 folds at time points of 2 and 8 h, respectively; HbRbohC2 was upregulated by 2.2 folds at the time point of 8 h; and, HbRbohN was upregulated by 3.6, 3.7, or 9.3 folds at three time points respectively. It is worth noting that the transcript level of HbRbohN was shown to be extremely low in unstressed CATAS93-114 leaves but considerably high in Reken501 (Fig. 6d). The response of HbRbohs to Microcyclus ulei and Corynespora

4. Discussion 4.1. The rubber tree genome encodes a relatively high number of Rboh family genes Despite the importance of RBOHs in ROS-mediated cellular signaling and various plant processes, to date, research is still focused on model plants such as arabidopsis and rice (Sagi and Fluhr, 2006; Wang et al., 2013a; Chang et al., 2016; Kaur and Pati, 2016, 2018; Kaur et al., 135

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Fig. 7. Expression profiles of HbRboh genes in ethephon-treated leaves as well as ethephon-induced senescent leaves. Significant difference was determined under the parameters “FDR < 0.001” and “log2Ratio ≥ 1”.

RBOHB/OG3, RBOHE/OG4, RBOHF/OG5, RBOHH/OG6, and RBOHN/ OG7. RBOHB, RBOHE, RBOHF, and RBOHH correspond to Group I, Group III, Group IV, and Group V as described by Marino et al. (2011). By contrast, the proposed Group II is clearly split into two phylogenetic groups, i.e. RBOHD and RBOHC, whereas RBOHN has not been reported ever before. When taking advantage of sequenced plant genomes to trace their origin, RBOHE was shown to be a primitive group that can be traced back to moss (Physcomitrella patens) (Rensing et al., 2008); RBOHD is more likely to appear in the ancient seed plants, since it can be found in gymnosperms (e.g. Picea abies) but not in spikemoss (Selaginella moellendorffii) (Banks et al., 2011); RBOHB and RBOHF may appear some time before the divergence of monocots and eudicots, since they are widely present in angiosperms but not Amborella trichopoda, the only living species on the sister lineage to all other flowering plants (Amborella Genome Project, 2013); similar to RBOHH, RBOHN can also be traced back to A. trichopoda, though it has lost in a high number of plant species including poplar, arabidopsis, grapevine, and rice; RBOHC may be generated along with the γ WGD, since it can be found in most eudicots but not in Aquilegia coerulea, a member of the basal-most eudicot clade. This means that core eudicots initially contain all seven groups including RBOHN.

2018). Both of them were shown to have experienced massive gene loss and chromosomal rearrangement after several recent WGDs (Bowers et al., 2003; Wang et al., 2005). In a previous study, five groups were proposed based on the phylogenetic analysis of ten AtRBOHs and nine OsRBOHs (Marino et al., 2011), however, the evolutionary relationships of RBOHs were not well resolved and novel groups seems to be present in other plant species. In the present study, genome-wide identification and manual curation of Rboh family genes were performed in rubber tree as well as five representative core eudicots. These species belong to different lineages of core eudicots (see Supplementary Fig. S3): rubber tree, cassava, castor bean, and physic nut are representatives of the Euphorbiaceae family in Malpighiales; poplar is a well-known model woody plant in the Salicaceae family of Malpighiales; grapevine is a Vitaceae plant in Parietales (Tuskan et al., 2006; Jaillon et al., 2007; Zou et al., 2016, 2018). The gene family of ten members in rubber tree is comparative to ten, ten, nine or nine present in poplar, arabidopsis, cassava, and rice respectively, but relatively more than seven found in castor bean, physic nut as well as grapevine. In contrast to no recent Rboh duplicate present in both castor bean and physic nut, one local duplicate was found in grapevine and a relatively high number of duplicates were identified in rubber tree, cassava, poplar, arabidopsis as well as rice (Table 2 and Supplementary Tables S1 and S4), reflecting the occurrence of one or more recent WGDs (Bowers et al., 2003; Wang et al., 2005; Tuskan et al., 2006; Bredeson et al., 2016; Pootakham et al., 2017). According to comparative genomics analyses, after the divergence of monocots and eudicots, the ancestor of rice experienced three WGDs named τ, σ, and ρ, whereas the last common ancestor of core eudicots experienced the γ WGD that was estimated to occur at approximately 117 million years ago (Mya) (Jiao et al., 2012, 2014). Moreover, the ancestral arabidopsis was proven to experience two more WGDs known as β and α, which was estimated to occur within a window of 61–65 or 23–50 Mya respectively (Bowers et al., 2003; Vanneste et al., 2014); poplar was proven to experience one recent WGD at 60–65 Mya, which is specific to Salicaceae (Tuskan et al., 2006); and, rubber tree and cassava were proven to share one recent doubling event at 39–47 Mya (Bredeson et al., 2016; Pootakham et al., 2017). By contrast, no recent WGD was found in castor bean, physic nut as well as grapevine (Jaillon et al., 2007; Chan et al., 2010; Wu et al., 2015). Phylogenetic and BRH-based homologous analyses divided 69 Rboh genes into seven groups or OGs, i.e. RBOHD/OG1, RBOHC/OG2,

4.2. WGDs act as the main force for the expansion of HbRboh gene family Although more than one genome assemblies have been available in rubber tree, the most complete one is derived from Reyan7-33-97, which consists of 7,453 scaffolds (Tang et al., 2016). Nevertheless, due to lack of a high-density genetic map in this species, it is not easy to determine whether a duplicate is derived from WGD or segmental duplications from the view of synteny analysis. In cassava, a genetic map with 22,403 markers has been established, which allows anchoring 89.0% of the AM560-2 assembly to 18 chromosomes (Bredeson et al., 2016). Since cassava shared the ρ WGD and harbors similar chromosome structures with rubber tree (Bredeson et al., 2016), it is possible to define WGD duplicates based on synteny analysis as well as calculating the Ks value of duplicate pairs. Duplicates in rubber tree and cassava were all shown to result from the recent WGD, which is similar to poplar and rice (Table 2, Supplementary Tables S1 and S4). By contrast, four arabidopsis duplicates were derived from the β WGD (1), α WGD (2) as well as proximal duplication (1) (Wang et al., 2013b). Moreover, species-specific expansion of different groups was observed in rubber tree as well as other examined species: the expansion of RBOHB and 136

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

be variety-specific.

RBOHE was only observed in rice; RBOHD has expanded in rubber tree, poplar, and grapevine; RBOHC has expanded in rubber tree, cassava, poplar, and arabidopsis; RBOHF has expanded in rubber tree, cassava, poplar, arabidopsis, and rice; RBOHH has expanded in poplar, arabidopsis, and rice (Table 3 and Supplementary Fig. S3). In addition to gene copies, structural divergence of different group members was also observed. Compared with other examined species, the exon-intron structure is relatively more variable in arabidopsis and rice due to chromosomal rearrangement after WGDs (Bowers et al., 2003; Wang et al., 2005). The result is consistent with the fast evolution of annual herbs and slow evolution of perennial trees and shrubs (Luo et al., 2015).

4.5. HbRboh gene expression is modified in TPD and leaf senescence Despite the involvement of RBOHs in TPD occurrence was proposed 30 years ago (Chrestin et al., 1984), the molecular mechanism underlying is yet to be resolved. In severe TPD trees without ETH stimulation, downregulation of HbRbohB, HbRbohC1, HbRbohD1, HbRbohF1, and HbRbohF2 in bark was observed. Nevertheless, in the laticifer of severe TPD trees with or without ETH stimulation, no significant changes were observed for both HbRbohF2 and HbRbohC1 that are expressed in this special tissue. By contrast, upregulation of HbRbohF2 in the laticifer of ETH-treated trees (Montoro et al., 2018) was observed, indicating a putative role of HbRboh genes on TPD may be due to transcription induction by ETH. However, possible post-translational modifications cannot be ruled out. According to recent evidence from cytology and molecular biology, TPD is more likely to be a type of programmed cell death (PCD) that was characterized by nuclear DNA fragmentation and significant regulation of senescence-associated genes (Chen et al., 2003; Zou et al., 2012). Indeed, overexpression of a Myb transcription factor HbMyb1 in tobacco can suppress stress-induced cell death (Peng et al., 2011). In barley, RBOHF2 was shown to participate in leaf age-dependent accelerated senescence (Torres et al., 2017). In rubber tree, significant regulation of HbRboh gene expression is not only found in natural senescent leaves (i.e. HbRbohD1 and HbRbohF1) but also in ETH-induced senescent leaves (i.e. HbRbohC2, HbRbohD1, and HbRbohF2). Among them, upregulation of HbRbohD1 is shared by two types of senescent leaves, implying its key role in leaf senescence. Moreover, HbRbohD1 as well as HbRbohF2 are induced by ETH, whereas HbRbohC2 is inhibited by ETH, implying their different regulatory mechanism. In the end, HbRboh gene expression is modified in TPD, leaf senescence and ETH treatment. Moreover, the fact that these genes are differentially regulated by different types of tissue/stress could indicate their roles in the redox status of a cell in diverse physiological traits.

4.3. Expression divergence of duplicated HbRboh genes and tissue-specific isoforms As shown in Figs. 3–4, several typical domains or motifs are highly conserved between duplicates in rubber tree as well as cassava. Thereby, expression divergence may play a more important role in functional diversification of duplicated HbRboh genes. Indeed, the transcript level of HbRbohC1, HbRbohD1, and HbRbohF2 is relatively higher than their paralogs in all examined tissues, which is similar to MeRbohC1 and MeRbohF2 acting as predominant isoforms in cassava. Similar expression patterns were also reported in arabidopsis and rice (Sagi and Fluhr, 2006; Wang et al., 2013a). In arabidopsis, AtRbohD and AtRbohF represent the predominant isoforms and are ubiquitously expressed, whereas AtRbohH and AtRbohJ are preferentially expressed in pollen and stamen (Sagi and Fluhr, 2006). HbRbohH, the ortholog of AtRbohH and AtRbohJ (Kaya et al., 2014) in rubber tree, was also shown to exhibit a male flower-specific expression pattern, implying their similar functions. Nevertheless, HbRbohC1 and HbRbohF2 represent the ubiquitously expressed genes in rubber tree, which is similar to MeRbohC1 and MeRbohF2 in cassava, implying species-specific evolution. Similar thing is also found for HbRbohB and MeRbohB: HbRbohB is highly expressed in root, in contrast, the transcript level of MeRbohB is not very high in fibrous root as well as starchy-enriched storage root. It’s worth noting that the expression of only two HbRboh genes (HbRbohF2 and HbRbohC1) was detected in the rubber-producing laticifer, where HbRbohF2 represents the predominant isoform.

5. Conclusions This paper presents the first genome-wide study of the Rboh gene family in rubber tree, which includes the comparative evolutionary analysis with seven representative plant species, i.e. cassava, physic nut, castor bean, poplar, grapevine, arabidopsis, and rice. Seven OGs corresponding to seven phylogenetic groups, i.e. RBOHD/OG1, RBOHC/OG2, RBOHB/OG3, RBOHE/OG4, RBOHF/OG5, RBOHH/OG6, and RBOHN/OG7, were proposed for core eudicots, and rubber paralogs were shown to be derived from the recent WGD shared by cassava. Tissue-specific expression of certain HbRboh genes and expression divergence of HbRboh duplicates were observed. Moreover, gene expression patterns seem to be characteristic of TPD, leaf senescence, hormone and stress responses in rubber tree. These findings will provide valuable information for further studies of Rboh genes in rubber tree as well as other species.

4.4. HbRboh gene expression is modified in hormone and stress induction Secondary laticifers, which are located in the secondary phloem and closely related to rubber productivity, periodically differentiate from vascular cambia. The progress was shown to be induced by mechanical wounding (e.g. tapping), exogenous JA and linolenic acid (Hao and Wu, 1982, 2000). Moreover, wounding-induced laticifer differentiation is dependent on the activity of RBOHs (Tian et al., 2015). Results showed that HbRbohB, HbRbohC1, HbRbohC2, HbRbohF2, and HbRbohN are significantly regulated by JA or its mimic COR. Among them, induced expression of HbRbohC1 by COR in early stages suggests its potential role in laticifer differentiation. Four genes (i.e. HbRbohF2, HbRbohF1, HbRbohD2, and HbRbohN) were shown to be significantly regulated by ETH, where induced expression of HbRbohF2 in laticifer and HbRbohF1 in bark implies their possible roles in ETH-stimulated yield promotion. It’s well documented that RBOHs are involved in plant-pathogen interaction and responses to abiotic stresses (Lightfoot et al., 2008; Marino et al., 2011, 2012; Chaouch et al., 2012; Cheng et al., 2013; Wang et al., 2013a; Torres et al., 2017). The present study showed that HbRbohC2, HbRbohN, HbRbohD1, HbRbohD2, and HbRbohF1 are related to the response to M. ulei and C. mycelium. In another study, induced expression of HbRboh genes upon Pseudocercospora ulei infection was also reported for a partially resistant genotype MDF180 (HbRbohB and HbRbohF2) and a totally resistant genotype FX2784 (HbRbohB) (Koop et al., 2016). Additionally, HbRbohC1, HbRbohC2, and HbRbohN are associated with cold stress, though the response pattern was shown to

Competing interests The authors declare that they have no competing interests.

Authors’ contributions The study was conceived and directed by ZZ. All the experiments and analysis were directed by ZZ and carried out by ZZ, JY and XZ. ZZ wrote the paper. All the authors read and approved the final manuscript. 137

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

Consent for publication

Fang, Y., Mei, H., Zhou, B., et al., 2016. De novo transcriptome analysis reveals distinct defense mechanisms by young and mature leaves of Hevea brasiliensis (para rubber tree). Sci. Rep. 6, 33151. https://doi.org/10.1038/srep33151. Groom, Q.J., Torres, M.A., Fordham-Skelton, A.P., Hammond-Kosack, K.E., Robinson, N.J., Jones, J.D., 1996. rbohA, a Rice homologue of the mammalian gp91phox respiratory burst oxidase gene. Plant J. 10, 515–522. Hao, B.Z., Wu, J.L., 1982. Effects of wound (tapping) on laticifer differentiation in Hevea brasiliensis. Acta Bot. Sin. 24, 388–391. Hao, B.Z., Wu, J.L., 2000. Laticifer differentiation in Hevea brasiliensis: induction by exogenous jasmonic acid and linolenic acid. Ann. Bot. 85, 37–43. Iqbal, N., Khan, N.A., Ferrante, A., Trivellini, A., Francini, A., Khan, M.I.R., 2017. Ethylene role in plant growth, development and senescence: interaction with other phytohormones. Front. Plant Sci. 8, 475. https://doi.org/10.3389/fpls.2017.00475. Hurtado Páez, U.A., García Romero, I.A., Restrepo Restrepo, S., Aristizábal Gutiérrez, F.A., Montoya Castaño, D., 2015. Assembly and analysis of differential transcriptome responses of Hevea brasiliensis on interaction with Microcyclus ulei. PLoS One 10, e0134837. https://doi.org/10.1371/journal.pone.0134837. Jaillon, O., Aury, J.M., Noel, B., et al., 2007. The grapevine genome sequence suggests ancestral hexaploidization in major angiosperm phyla. Nature 449, 463–467. Jiao, Y., Leebens-Mack, J., Ayyampalayam, S., et al., 2012. A genome triplication associated with early diversification of the core eudicots. Genome Biol. 13, R3. https:// doi.org/10.1186/gb-2012-13-1-r3. Jiao, Y., Li, J., Tang, H., Paterson, A.H., 2014. Integrated syntenic and phylogenomic analyses reveal an ancient genome duplication in monocots. Plant Cell 26, 2792–2802. https://doi.org/10.1105/tpc.114.127597. Kaur, G., Sharma, A., Guruprasad, K., Pati, P.K., 2014. Versatile roles of plant NADPH oxidases and emerging concepts. Biotechnol. Adv. 32, 551–563. https://doi.org/10. 1016/j.biotechadv.2014.02.002. Kaur, G., Pati, P.K., 2016. Analysis of cis-acting regulatory elements of Respiratory burst oxidase homolog (Rboh) gene families in Arabidopsis and rice provides clues for their diverse functions. Comput. Biol. Chem. 62, 104–118. https://doi.org/10.1016/j. compbiolchem.2016.04.002. Kaur, G., Pati, P.K., 2018. In silico physicochemical characterization and topology analysis of Respiratory burst oxidase homolog (Rboh) proteins from Arabidopsis and rice. Bioinformation 14, 93–100. https://doi.org/10.6026/97320630014093. Kaur, G., Guruprasad, K., Temple, B.R.S., Shirvanyants, D.G., Dokholyan, N.V., Pati, P.K., 2018. Structural complexity and functional diversity of plant NADPH oxidases. Amino Acids 50, 79–94. https://doi.org/10.1007/s00726-017-2491-5. Kaya, H., Nakajima, R., Iwano, M., et al., 2014. Ca2+-activated reactive oxygen species production by Arabidopsis RbohH and RbohJ is essential for proper pollen tube tip growth. Plant Cell 26, 1069–1080. https://doi.org/10.1105/tpc.113.120642. Koop, D.M., Rio, M., Sabau, X., et al., 2016. Expression analysis of ROS producing and scavenging enzyme-encoding genes in rubber tree infected by Pseudocercospora ulei. Plant Physiol. Biochem. 104, 188–199. https://doi.org/10.1016/j.plaphy.2016.03. 022. Langmead, B., Salzberg, S.L., 2012. Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. Li, D., Deng, Z., Chen, C., et al., 2010. Identification and characterization of genes associated with tapping panel dryness from Hevea brasiliensis latex using suppression subtractive hybridization. BMC Plant Biol. 10, 140. https://doi.org/10.1186/14712229-10-140. Li, D., Wang, X., Deng, Z., Liu, H., Yang, H., He, G., 2016. Transcriptome analyses reveal molecular mechanism underlying tapping panel dryness of rubber tree (Hevea brasiliensis). Sci. Rep. 6, 23540. https://doi.org/10.1038/srep23540. Lieberei, R., 2007. South American leaf blight of the rubber tree (Hevea spp.): new steps in plant domestication using physiological features and molecular markers. Ann. Bot. 100, 1125–1142. Lightfoot, D.J., Boettcher, A., Little, A., Shirley, N., Able, A.J., 2008. Identification and characterisation of barley (Hordeum vulgare) respiratory burst oxidase homologue family members. Funct. Plant Biol. 35, 347–359. https://doi.org/10.1071/FP08109. Liu, J.P., 2016. Molecular mechanism underlying ethylene stimulation of latex production in rubber tree (Hevea brasiliensis). Trees 30, 1913–1921. https://doi.org/10.1007/ s00468-016-1455-9. Liu, J.P., Xia, Z.Q., Tian, X.Y., Li, Y.J., 2015. Transcriptome sequencing and analysis of rubber tree (Hevea brasiliensis Muell.) to discover putative genes associated with tapping panel dryness (TPD). BMC Genomics 16, 398. https://doi.org/10.1186/ s12864-015-1562-9. Liu, J.P., Zhuang, Y.F., Guo, X.L., Li, Y.J., 2016. Molecular mechanism of ethylene stimulation of latex yield in rubber tree (Hevea brasiliensis) revealed by de novo sequencing and transcriptome analysis. BMC Genomics 17, 257. https://doi.org/10. 1186/s12864-016-2587-4. Liu, J.P., Hu, J., Liu, Y.H., et al., 2018. Transcriptome analysis of Hevea brasiliensis in response to exogenous methyl jasmonate provides novel insights into regulation of jasmonate-elicited rubber biosynthesis. Physiol. Mol. Biol. Plants 24, 349–358. https://doi.org/10.1007/s12298-018-0529-0. Luo, M.C., You, F.M., Li, P., et al., 2015. Synteny analysis in Rosids with a walnut physical map reveals slow genome evolution in long-lived woody perennials. BMC Genomics 16, 707. https://doi.org/10.1186/s12864-015-1906-5. Magnani, F., Nenci, S., Millana Fananas, E., et al., 2017. Crystal structures and atomic model of NADPH oxidase. Proc. Natl. Acad. Sci. U. S. A. 114, 6764–6769. https://doi. org/10.1073/pnas.1702293114. Marino, D., Andrio, E., Danchin, E.G., et al., 2011. A Medicago truncatula NADPH oxidase is involved in symbiotic nodule functioning. New Phytol. 189, 580–592. https://doi. org/10.1111/j.1469-8137.2010.03509.x. Marino, D., Dunand, C., Puppo, A., Pauly, N., 2012. A burst of plant NADPH oxidases. Trends Plant Sci. 17, 9–15. https://doi.org/10.1016/j.tplants.2011.10.001.

All authors consent for publication. Acknowledgements This work was supported by the National Natural Science Foundation of China (31700580 and 31371556), the Natural Science Foundation of Hainan province (31100460) and the Central Public-interest Scientific Institution Basal Research Fund for Chinese Academy of Tropical Agricultural Sciences (1630022011014 and 1630022017030). The authors would like to thank three reviewers for their careful reading and valuable suggestions. The authors also appreciate those contributors who make related genome and transcriptome data accessible in public databases. Appendix A. Supplementary data Supplementary material related to this article can be found, in the online version, at doi:https://doi.org/10.1016/j.indcrop.2018.11.005. References Altschul, S.F., Madden, T.L., Schaffer, A.A., et al., 1997. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. 25, 3389–3402. Amborella Genome Project, 2013. The Amborella genome and the evolution of flowering plants. Science 342, 1241089. https://doi.org/10.1126/science.1241089. Audic, S., Claverie, J.M., 1997. The significance of digital gene expression profiles. Genome Res. 7, 986–995. Bailey, T.L., Boden, M., Buske, F.A., et al., 2009. MEME SUITE: tools for motif discovery and searching. Nucleic Acids Res. 37, W202–208. Banks, J.A., Nishiyama, T., Hasebe, M., et al., 2011. The Selaginella genome identifies genetic changes associated with the evolution of vascular plants. Science 332, 960–963. https://doi.org/10.1126/science.1203810. Bedard, K., Lardy, B., Krause, K.H., 2007. NOX family NADPH oxidases: not just in mammals. Biochimie 89, 1107–1112. Bowers, J.E., Chapman, B.A., Rong, J., Paterson, A.H., 2003. Unravelling angiosperm genome evolution by phylogenetic analysis of chromosomal duplication events. Nature 422, 433–438. Bredeson, J.V., Lyons, J.B., Prochnik, S.E., et al., 2016. Sequencing wild and cultivated cassava and related species reveals extensive interspecific hybridization and genetic diversity. Nat. Biotechnol. 34, 562–570. https://doi.org/10.1038/nbt.3535. Chan, A.P., Crabtree, J., Zhao, Q., et al., 2010. Draft genome sequence of the oilseed species Ricinus communis. Nat. Biotechnol. 28, 951–956. Chang, Y.L., Li, W.Y., Miao, H., et al., 2016. Comprehensive genomic analysis and expression profiling of the NOX gene families under abiotic stresses and hormones in plants. Genome Biol. Evol. 8, 791–810. https://doi.org/10.1093/gbe/evw035. Chao, J., Zhang, S., Chen, Y., Tian, W.M., 2015. Cloning, heterologous expression and characterization of ascorbate peroxidase (APX) gene in laticifer cells of rubber tree (Hevea brasiliensis Muell. Arg.). Plant Physiol. Biochem. 97, 331–338. https://doi.org/ 10.1016/j.plaphy.2015.10.023. Chaouch, S., Queval, G., Noctor, G., 2012. AtRbohF is a crucial modulator of defenceassociated metabolism and a key actor in the interplay between intracellular oxidative stress and pathogenesis responses in Arabidopsis. Plant J. 69, 613–627. https:// doi.org/10.1111/j.1365-313X.2011.04816.x. Chen, S.C., Peng, S.Q., Huang, G.X., Wu, K., Fu, X., Chen, Z., 2003. Association of decreased expression of a Myb transcription factor with the TPD (tapping panel dryness) syndrome in Hevea brasiliensis. Plant Mol. Biol. 51, 51–58. Cheng, C., Xu, X., Gao, M., et al., 2013. Genome-wide analysis of respiratory burst oxidase homologs in grape (Vitis vinifera L.). Int. J. Mol. Sci. 14, 24169–24186. https://doi. org/10.3390/ijms141224169. Cheng, H., Chen, X., Fang, J., An, Z., Hu, Y., Huang, H., 2018. Comparative transcriptome analysis reveals an early gene expression profile that contributes to cold resistance in Hevea brasiliensis (the Para rubber tree). Tree Physiol. 38, 1409–1423. https://doi. org/10.1093/treephys/tpy014. Chou, K.C., Shen, H.B., 2010. Cell-PLoc: a package of web-servers for predicting subcellular localization of proteins in various organisms. Nat. Protoc. 3, 153–162. Chrestin, H., Bangratz, J., d’Auzac, J., Jacob, J., 1984. Role of the lutoidic tonoplast in the senescence and degeneration of the laticifers of Hevea brasiliensis. Zeitschrift für Pflanzenphysio 114, 261–268. Deng, Z., Zhao, M., Liu, H., Wang, Y., Li, D., 2015. Molecular cloning, expression profiles and characterization of a glutathione reductase in Hevea brasiliensis. Plant Physiol. Biochem. 96, 53–63. https://doi.org/10.1016/j.plaphy.2015.07.022. Dietz, K.J., Mittler, R., Noctor, G., 2016. Recent progress in understanding the role of reactive oxygen species in plant cell signaling. Plant Physiol. 171, 1535–1539. https://doi.org/10.1104/pp.16.00938. Edgar, R.C., 2004. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 32, 1792–1797.

138

Industrial Crops & Products 128 (2019) 126–139

Z. Zou et al.

tree (Hevea brasiliensis Muell. Arg.). Planta 226, 499–515. Wang, Y., Tang, H., Debarry, J.D., et al., 2012. MCScanX: a toolkit for detection and evolutionary analysis of gene synteny and collinearity. Nucleic Acids Res. 40, e49. https://doi.org/10.1093/nar/gkr1293. Wang, G.F., Li, W.Q., Li, W.Y., Wu, G.L., Zhou, C.Y., Chen, K.M., 2013a. Characterization of Rice NADPH oxidase genes and their expression under various environmental conditions. Int. J. Mol. Sci. 14, 9440–9458. https://doi.org/10.3390/ijms14059440. Wang, W., Chen, D., Zhang, X., Liu, D., Cheng, Y., Shen, F., 2018. Role of plant respiratory burst oxidase homologs in stress responses. Free Radic. Res. 52, 826–839. https://doi. org/10.1080/10715762.2018.1473572. Wang, X., Shi, X., Hao, B., Ge, S., Luo, J., 2005. Duplication and DNA segmental loss in the rice genome: implications for diploidization. New Phytol. 165, 937–946. Wang, Y., Tan, X., Paterson, A.H., 2013b. Different patterns of gene structure divergence following gene duplication in Arabidopsis. BMC Genomics 14, 652. https://doi.org/ 10.1186/1471-2164-14-652. Wilson, M.C., Mutka, A.M., Hummel, A.W., et al., 2017. Gene expression atlas for the food security crop cassava. New Phytol. 213, 1632–1641. https://doi.org/10.1111/nph. 14443. Wu, P., Zhou, C., Cheng, S., et al., 2015. Integrated genome sequence and linkage map of physic nut (Jatropha curcas L.), a biodiesel plant. Plant J. 81, 810–821. Wu, S., Zhang, S., Chao, J., et al., 2016. Transcriptome analysis of the signalling networks in coronatine-induced secondary laticifer differentiation from vascular cambia in rubber trees. Sci. Rep. 6, 36384. https://doi.org/10.1038/srep36384. Yang, Z., 2007. PAML 4: phylogenetic analysis by maximum likelihood. Mol. Biol. Evol. 24, 1586–1591. Yun, B.W., Feechan, A., Yin, M., et al., 2011. S-nitrosylation of NADPH oxidase regulates cell death in plant immunity. Nature 478, 264–268. https://doi.org/10.1038/ nature10427. Zhang, Y., Zhu, H., Zhang, Q., et al., 2009. Phospholipase dalpha1 and phosphatidic acid regulate NADPH oxidase activity and production of reactive oxygen species in ABAmediated stomatal closure in Arabidopsis. Plant Cell 21, 2357–2377. https://doi.org/ 10.1105/tpc.108.062992. Zhang, S.X., Wu, S.H., Chen, Y.Y., Tian, W.M., 2015. Analysis of differentially expressed genes associated with coronatine-induced laticifer differentiation in the rubber tree by subtractive hybridization suppression. PLoS One 10, e0132070. https://doi.org/ 10.1371/journal.pone.0132070. Zou, Z., Huang, Q.X., Xie, G.S., Yang, L.F., 2018. Genome-wide comparative analysis of papain-like cysteine protease family genes in castor bean and physic nut. Sci. Rep. 8, 331. https://doi.org/10.1038/s41598-017-18760-6. Zou, Z., Gong, J., An, F., et al., 2015. Genome-wide identification of rubber tree (Hevea brasiliensis Muell. Arg.) aquaporin genes and their response to ethephon stimulation in the laticifer, a rubber-producing tissue. BMC Genomics 16, 1001. https://doi.org/ 10.1186/s12864-015-2152-6. Zou, Z., Liu, J., Yang, L., Xie, G., 2017a. Survey of the rubber tree genome reveals a high number of cysteine protease-encoding genes homologous to Arabidopsis SAG12. PLoS One 12, e0171725,. https://doi.org/10.1371/journal.pone.0171725. Zou, Z., Xie, G.S., Yang, L.F., 2017b. Papain-like cysteine protease encoding genes in rubber (Hevea brasiliensis): comparative genomics, phylogenetic and transcriptional profiling analysis. Planta 246, 999–1018. https://doi.org/10.1007/s00425-0172739-z. Zou, Z., Yang, L.F., Gong, J., et al., 2016. Genome-wide identification of Jatropha curcas aquaporin genes and the comparative analysis provides insights into the gene family expansion and evolution in Hevea brasiliensis. Front.. Plant Sci. 7, 395. https://doi. org/10.3389/fpls.2016.00395. Zou, Z., Yang, L.F., Wang, Z.H., Yuan, K., 2009. Biosynthesis and regulation of natural rubber in Hevea. Plant Physiol. J. 45, 1231–1238. Zou, Z., Yang, L.F., Wang, Z.H., Yuan, K., 2012. Strategies on prevention and treatment of tapping panel dryness syndrome of rubber (Hevea brasiliensis Muell. Arg.). J. Chin. Biotechnol. 9, 8–15.

Montoro, P., Wu, S., Favreau, B., Herlinawati, E., et al., 2018. Transcriptome analysis in Hevea brasiliensis latex revealed changes in hormone signalling pathways during ethephon stimulation and consequent Tapping Panel Dryness. Sci. Rep. 8, 8483. https://doi.org/10.1038/s41598-018-26854-y. Mortazavi, A., Williams, B.A., McCue, K., Schaeffer, L., Wold, B., 2008. Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. https://doi.org/10.1038/nmeth.1226. Oda, T., Hashimoto, H., Kuwabara, N., et al., 2010. Structure of the N-terminal regulatory domain of a plant NADPH oxidase and its functional implications. J. Biol. Chem. 285, 1435–1445. https://doi.org/10.1074/jbc.M109.058909. Ogasawara, Y., Kaya, H., Hiraoka, G., et al., 2008. Synergistic activation of the Arabidopsis NADPH oxidase AtRbohD by Ca2+ and phosphorylation. J. Biol. Chem. 283, 8885–8892. https://doi.org/10.1074/jbc.M708106200. Olsen, K., Schaal, B., 2001. Microsatellite variation in cassava (Manihot esculenta, Euphorbiaceae) and its wild relatives: further evidence for a southern Amazonian origin of domestication. Am. J. Bot. 88, 131–142. Peng, S.Q., Wu, K.X., Huang, G.X., Chen, S.C., 2011. HbMyb1, a Myb transcription factor from Hevea brasiliensis, suppresses stress induced cell death in transgenic tobacco. Plant Physiol. Biochem. 49, 1429–1435. https://doi.org/10.1016/j.plaphy.2011.09. 007. Pootakham, W., Sonthirod, C., Naktang, C., et al., 2017. De novo hybrid assembly of the rubber tree genome reveals evidence of paleotetraploidy in Hevea species. Sci. Rep. 7, 41457. https://doi.org/10.1038/srep41457. Putranto, R.A., Herlinawati, E., Rio, M., et al., 2015. Involvement of ethylene in the latex metabolism and tapping panel dryness of Hevea brasiliensis. Int. J. Mol. Sci. 16, 17885–17908. https://doi.org/10.3390/ijms160817885. Rensing, S.A., Lang, D., Zimmer, A.D., 2008. The Physcomitrella genome reveals evolutionary insights into the conquest of land by plants. Science 319, 64–69. Sagi, M., Fluhr, R., 2006. Production of reactive oxygen species by plant NADPH oxidases. Plant Physiol. 141, 336–340. Sumimoto, H., 2008. Structure, regulation and evolution of Nox-family NADPH oxidases that produce reactive oxygen species. FEBS J. 275, 3249–3277. https://doi.org/10. 1111/j.1742-4658.2008.06488.x. Suzuki, N., Miller, G., Morales, J., Shulaev, V., Torres, M.A., Mittler, R., 2011. Respiratory burst oxidases: the engines of ROS signaling. Curr. Opin. Plant Biol. 14, 691–699. https://doi.org/10.1016/j.pbi.2011.07.014. Tamura, K., Stecher, G., Peterson, D., Filipski, A., Kumar, S., 2013. MEGA6: Molecular Evolutionary Genetics Analysis version 6.0. Mol. Biol. Evol. 30, 2725–2729. Tang, C., Yang, M., Fang, Y., et al., 2016. The rubber tree genome reveals new insights into rubber production and species adaptation. Nat. Plants 2, 16073. https://doi.org/ 10.1038/nplants.2016.73. Tian, W.M., Yang, S.G., Shi, M.J., Zhang, S.X., Wu, J.L., 2015. Mechanical woundinginduced laticifer differentiation in rubber tree: an indicative role of dehydration, hydrogen peroxide, and jasmonates. J. Plant Physiol. 182, 95–103. https://doi.org/ 10.1016/j.jplph.2015.04.010. Torres, M.A., Dangl, J.L., Jones, J.D., 2002. Arabidopsis gp91phox homologues AtRbohD and AtRbohF are required for accumulation of reactive oxygen intermediates in the plant defense response. Proc. Natl. Acad. Sci. U. S. A. 99, 517–522. Torres, D.P., Proels, R.K., Schempp, H., Hückelhoven, R., 2017. Silencing of RBOHF2 causes leaf age-dependent accelerated senescence, salicylic acid accumulation, and powdery mildew resistance in barley. Mol. Plant Microbe Interact. 30, 906–918. https://doi.org/10.1094/MPMI-04-17-0088-R. MPMI04170088R. Tuskan, G.A., Difazio, S., Jansson, S., et al., 2006. The genome of black cottonwood, Populus trichocarpa (Torr. & Gray). Science 313, 1596–1604. Vanneste, K., Baele, G., Maere, S., Van de Peer, Y., 2014. Analysis of 41 plant genomes supports a wave of successful genome duplications in association with the Cretaceous-Paleogene boundary. Genome Res. 24, 1334–1347. https://doi.org/10. 1101/gr.168997.113. Venkatachalam, P., Thulaseedharan, A., Raghothama, K., 2007. Identification of expression profiles of tapping panel dryness (TPD) associated genes from the latex of rubber

139