Transcriptome Analysis of Distinct Cold Tolerance Strategies in the Rubber Tree (Hevea brasiliensis)

Natural rubber is an indispensable commodity used in approximately 40,000 products and is fundamental to the tire industry. Among the species that produce latex, the rubber tree [Hevea brasiliensis (Willd. ex Adr. de Juss.) Muell-Arg.], a species native to the Amazon rainforest, is the major producer of latex used worldwide. The Amazon Basin presents optimal conditions for rubber tree growth, but the occurrence of South American leaf blight, which is caused by the fungus Microcyclus ulei (P. Henn) v. Arx, limits rubber tree production. Currently, rubber tree plantations are located in scape regions that exhibit suboptimal conditions such as high winds and cold temperatures. Rubber tree breeding programs aim to identify clones that are adapted to these stress conditions. However, rubber tree breeding is time-consuming, taking more than 20 years to develop a new variety. It is also expensive and requires large field areas. Thus, genetic studies could optimize field evaluations, thereby reducing the time and area required for these experiments. Transcriptome sequencing using next-generation sequencing (RNA-seq) is a powerful tool to identify a full set of transcripts and for evaluating gene expression in model and non-model species. In this study, we constructed a comprehensive transcriptome to evaluate the cold response strategies of the RRIM600 (cold-resistant) and GT1 (cold-tolerant) genotypes. Furthermore, we identified putative microsatellite (SSR) and single-nucleotide polymorphism (SNP) markers. Alternative splicing, which is an important mechanism for plant adaptation under abiotic stress, was further identified, providing an important database for further studies of cold tolerance.

Natural rubber is an indispensable commodity used in approximately 40,000 products and is 21 fundamental to the tire industry. Among the species that produce latex, the rubber tree [Hevea 22 brasiliensis (Willd. ex Adr. de Juss.) Muell-Arg.], a species native to the Amazon rainforest, is 23 the major producer of latex used worldwide. The Amazon Basin presents optimal conditions for 24 rubber tree growth, but the occurrence of South American leaf blight, which is caused by the 25 fungus Microcyclus ulei (P. Henn) v. Arx, limits rubber tree production. Currently, rubber tree 26 plantations are located in scape regions that exhibit suboptimal conditions such as high winds 27 and cold temperatures. Rubber tree breeding programs aim to identify clones that are adapted to 28 these stress conditions. However, rubber tree breeding is time-consuming, taking more than 20 29 years to develop a new variety. It is also expensive and requires large field areas. Thus, genetic 30 studies could optimize field evaluations, thereby reducing the time and area required for these 31 experiments. Transcriptome sequencing using next-generation sequencing (RNA-seq) is a 32 powerful tool to identify a full set of transcripts and for evaluating gene expression in model and 33 non-model species. In this study, we constructed a comprehensive transcriptome to evaluate the 34 cold response strategies of the RRIM600 (cold-resistant) and GT1 (cold-tolerant) genotypes. 35 Furthermore, we identified putative microsatellite (SSR) and single-nucleotide polymorphism 36 (SNP) markers. Alternative splicing, which is an important mechanism for plant adaptation under 37 abiotic stress, was further identified, providing an important database for further studies of cold 38 tolerance. 39

1
Introduction 40 Cold stress, which can be classified as chilling (0 to 15°C) and/or freezing (< 0°C) temperatures, 41 affects plant growth and development, limiting spatial distribution and yields (Chinnusamy et al.,42 2007; Yang et al., 2015). Furthermore, cold stress prevents plants from achieving their full 43 genetic potential, inhibiting metabolic reactions, reducing photosynthetic capacity and altering 44 membrane permeability (Chinnusamy et al., 2007;Sevillano et al., 2009). 45 Temperate plants can generally achieve cold acclimation and acquire tolerance to extracellular 46 ice formation in their vegetative tissues. However, tropical crops, such as maize and rice, lack 47 cold acclimation ability and are sensitive to chilling (Chinnusamy et al., 2007). Furthermore, 48 varieties from the same species can exhibit different levels of cold tolerance (Liu et al., 2012; 49 Zhang et al., 2012). Hence, determining the gene expression profile under cold stress could help 50 to elucidate the mechanism of cold acclimation in plants and can be effective for selecting 51 candidate genes. Moreover, candidate genes can be targeted to identify genetic variation and 52 develop molecular markers. 53 Hevea brasiliensis [(Willd. ex Adr. de Juss.) Muell-Arg], commonly known as the rubber tree, is 54 a perennial tree crop native to the Amazon rainforest. The species, belonging to the 55 Euphorbiaceae family, is monoecious, undergoes cross-pollination and has a chromosome 56 number of 2n = 2x = 36. Among the 2,500 species that produce natural rubber (cis-1,4-57 polyisoprene), H. brasiliensis is the only species that produces high-quality rubber in 58 commercially viable quantities, accounting for more than 98% of total worldwide production 59 (Pootakham et al., 2017). 60 Natural rubber is one of the most important raw materials for many industries and cannot be 61 replaced by synthetic alternatives due to its unique properties, such as flexibility, elasticity and 62 abrasion resistance (Sakdapipanich, 2007). Natural rubber is an essential commodity for the tire 63 industry and for the manufacture of more than 40,000 products. 64 Although the Amazon basin offers a suitable climate for this crop, Southeast Asia is the major 65 producer of rubber, being responsible for 92% of worldwide production. South America is 66 responsible for only 2% of worldwide rubber production, due to the occurrence of the fungus 67 Microcyclus ulei (P. Henn) v. Arx, which causes South American leaf blight (SALB). SALB was 68 responsible for devastating plantations in northern Brazil in the 1930s and remains a permanent 69 threat to the rubber industry (Pushparajah, 2001). To date, the rubber tree plantations in 70 Southeast Asia have not been affected by SALB, but other native pathogenic fungi are threats to 71 rubber production. The two major fungal pathogens in Southeast Asia (Phytophthora and 72 Corynespora) cause leaf fall and, consequently, significant losses of natural rubber yields 73 (Pootakham et al., 2017). Due to the occurrence of diseases, plantations have been expanded to 74 sub-optimal areas of some countries, such as northeastern India, the highlands and coastal areas 75 of Vietnam, southern China and the southern plateau of Brazil (Priyadarshan et al., 2005). These 76 areas are characterized by new stressful conditions, such as cold and dry periods. 77 The exposure of rubber trees to low temperatures can cause leaf necrosis, affecting tree 78 development and latex production (Priyadarshan et al., 2005;Mai et al., 2010). In addition, low 79 temperatures are responsible for halting latex production for 1-3 months per year (Rao et al., 80 1998;Jacob et al., 1999). 81 In recent years, there has been an exponential increase in genomic data acquisition for the rubber 82 tree, including transcriptome profiles ( expression mean-variance trend to facilitate Bayesian-moderated, weighted t-stastitics (Soneson 174 and Delorenzi, 2013), with at least 10 counts per million (CPM) in at least 3 samples. Three 175 biological replicates for each condition were provided for this analysis. We considered a gene to 176 be differentially expressed using a false discovery rate (FDR) cut-off ≤ 0.05. The pairwise 177 comparison was performed between RRIM600 and GT1 for each time-series of the cold 178 treatment: RRIM600 0h x GT1 0h, RRIM600 90min x GT1 90min, RRIM600 12h x GT1 12h, 179 RRIM600 24h x GT1 24h. 180

Gene Ontology Enrichment 181
The DEGs identified previously were subjected to GO enrichment analysis using GOseq with a 182 FDR cut-off ≤ 0.05 (Young et al., 2010

Alternative Splicing Identification 218
The filtered transcripts and AS events defined by PASA were processed using an in-house 219 pipeline. This pipeline identifies and re-classifies AS events simultaneously encompassing 220 alternative 5" and 3" splice sites (Chamala et al., 2015). Furthermore, a minimum of 10 reads 221 mapped at the splice junction was set as the threshold for considering an AS event. 222

Sequencing and Transcriptome Assembly 224
In the present study, we sequenced leaf tissue from the RRIM600 and GT1 genotypes. were considered a reference transcriptome and employed for further analysis. 247

Functional Annotation 248
The transcriptome was subjected to BLAST searches against the SwissProt/UniProt, Pfam, Gene 249 Ontology and KEGG databases with an e-value cut-off of 1e-5. In total, 63,983 (61%) transcripts 250 were annotated against the SwissProt/UniProt database using BLASTX. In addition, among the The GO annotation retrieved a total of 58,500 terms for the Biological Process (BP) category, 271 61,467 terms for Molecular Function (MF) and 62,795 terms for Celular Component (CC). In 272 addition, a total of 62,077 transcripts were annotated in the KEGG database. 273

Digital Gene Expression Analysis 274
To investigate the chilling stress response strategy between the rubber tree genotypes, we 275 performed a pairwise comparison for each time series between GT1 and RRIM600. Prior to 276 exposing the plants to cold stress (0 h), 624 genes were up-regulated in RRIM600 relative to 277 GT1 (RRIM600 0h x GT1 0h), while 732 genes were up-regulated in GT1 0h relative to 278 RRIM600 0h. After 90 minutes of cold stress exposure, we identified 514 genes that were up-279 regulated in RRIM600 90 min compared to GT1 90 min and 854 up-regulated genes in GT1 90 280 min relative to RRIM600 90min. Moreover, a total of 569 genes and 1034 genes were up-281 regulated in RRIM600 after 12 h and 24 h of cold stress exposure, respectively. Nevertheless, the 282 GT1 genotype exhibited 610 and 875 up-regulated genes relative to RRIM600 after 12 h and 24 283 h of cold treatment, respectively (Figures 2 and 3). 284 Since the DEG analysis was performed by comparing the RRIM600 and GT1 genotypes for each 285 time point of cold treatment, we compared the up-regulated genes previously identified for each 286 genotype across all time points in order to identify genes that were commonly and exclusively 287 up-regulated in the relative comparison between each treatment time. For the RRIM600 288 genotype, we detected a total of 229 genes that were exclusively up-regulated at 0 h ( Figure 3A). 289 After 90 minutes, 12 h and 24 h, we identified 100, 125 and 567 exclusively up-regulated genes, 290 respectively. Moreover, a total of 208 RRIM600 genes were commonly up-regulated within the 291 entire series ( Figure 3A and C). 292 Whereas, there were a total of 143 up-regulated genes identified in GT1 relative to RRIM600 293 that were exclusively up-regulated at time 0 h and 235, 110 and 390 genes that were exclusively 294 up-regulated after 90 minutes, 12 h and 24 h, respectively ( Figure 3B). Furthermore, a total of 295 255 genes were commonly up-regulated across the time series ( Figure 3B and D). 296 One of the genes up-regulated in RRIM600 relative to GT1 after 90minutes of cold stress was 297 identified as a putative stelar K + "outward-rectifying" channel (SKOR) based upon a high 298 BLAST similarity. Furthermore, two SKOR genes were up-regulated after 12 and 24 h of cold 299 stress, one of which was among the ten most highly expressed genes after 12 h of treatment 300 (  (Takahashi and Shimosaka, 1997). Notably, the log-fold-change in 313 RRIM600 24h relative to GT1 24h was greater than the log-fold-change at earlier periods of cold 314 treatment (Supplementary Table 2). 315 We also observed up-regulated genes in GT1 24h relative to RRIM600 24h , which were related 316 to cell growth, such as the gene with high similarity to the LONGIFOLIA 1 (LNG1) protein, 317 which in association with LONGIFOLIA 2 (LNG2) regulates leaf morphology by promoting 318 longitudinal  Table 2). 326

Protein Domain Homology among DEGs 327
Prior to cold stress, we detected four genes up-regulated in RRIM600 0h relative to GT1 0h with 328 the Apetala 2 (AP2) domain. The AP2 genes show high similarity to ERF119, which may be 329 involved in the regulation of gene expression by stress factors and by components of stress signal 330 transduction pathways. The IQ calmodulin-binding motif domain was detected in six up-331 regulated genes in GT1 0h relative to RRIM600 0h. At 90 minutes, the number of up-regulated 332 genes in GT1 containing this domain increased to eight (Table 4). The IQ calmodulin-binding is 333 a major calcium (Ca 2+ ) sensor and orchestrator of regulatory events through its interaction with a 334 diverse group of cellular proteins (Rhoads and Friedberg, 1997) (Supplementary Table 2). 335 We also identified five up-regulated genes in RRIM600 12h compared to GT1 12h containing 336 the VQ motif (Table 4). This domain is a plant-specific domain characteristic of a class of 337 proteins that regulates diverse developmental processes, including responses to biotic and abiotic 338 stresses, seed development, and photomorphogenesis (Jing and Lin, 2015). 339 After 24 h of cold treatment, the most abundant domain among the up-regulated genes in 340 RRIM600 24h was the NB-ARC domain, which is a signaling motif that is shared by plant 341 resistance products and is a regulator of cell death in animals (van der Biezen and Jones, 1998  (Table 5). 436 The up-regulated genes detected in the GT1 and RRIM600 genotypes were merged to identify 437 putative SSRs. A total of 1,034 dinucleotide SSRs were identified, followed by 629 tri-,18 tetra-438 ,19 penta-and 23 hexanucleotide SSRs (

qRT-PCR Validation 454
To validate the DEG analysis, a total of 20 genes were selected, and primer pairs were designed 455 to validate the analysis (Supplementary Table 1). All primer pairs were initially tested via PCR 456 using genomic DNA as a template to verify the amplification product. From the 20 primer pairs, 457 14 were successfully amplified and used for qRT-PCR. 458 The qRT-PCR assays were conducted using RH2b and YSL8 as housekeeping genes. were greater in GT1 than RRIM600, nevertheless this difference was not significant (Figure 4). 467 The only gene that is not in accordance with the in silico analysis was the protein ETHYLENE 468 INSENSITIVE 3 (EIN3) (PASA_cluster_52015). The in silico analysis showed higher 469 expression levels for EIN3 in RRIM600; however, qRT-PCR results showed that this gene is 470 significantly up-regulated in GT1 (Figure 4). 471

Alternative Splicing Detection 472
AS is an important mechanism involved in gene regulation that may regulate many physiological 473 processes in plants, including the response to abiotic stresses such as cold stress ( Due to the importance of AS, the reference transcriptome obtained in this study was also used to 478 detect AS events. A minimum depth of 10 reads was used as the threshold to identify isoforms. 479 A total of 20,279 AS events were identified, with intron retention (IR) representing the major AS 480 event, accounting for a total of 9,226 events (45.5%), followed by exon skipping (ES), 481 alternative acceptor (AltA) and alternative donor (AltD) events, at 4,806 (23.7%), 3,599 (17.7%) 482 and 2,648 (13%) events, respectively ( Figure 5). 483 Although the ES type accounts for the majority of AS events in humans, it has been reported that 484 IR events are the most abundant type in plants (Chamala et al., 2015). Before exposing the plants to cold stress in the present study, both genotypes presented up-511 regulated genes with similarity to peroxidase. However, GT1 exhibited two up-regulated gene, 512 while RRIM600 exhibited only one. After 90 minutes, only GT1 exhibited up-regulated 513 peroxidase genes. At this time point, we detected a total of four up-regulated genes with the best 514 BLAST hit to peroxidase. However, there were three up-regulated genes in RRIM600 with 515 probable peroxidase activity after 24h of treatment ( Figure 6, Supplementary Table 2). 516 SODs are the first line of defense of ROS (Alscher et al., 2002) and we identified two up-517 regulated putative SOD genes in both genotypes. The putative SOD gene that was up-regulated 518 in RRIM600 was detected after 12 h of cold stress, while the SOD gene in GT1 was up-regulated 519 after 24 h of treatment ( Figure 6, Supplementary Table 2). Furthermore, putative GST genes 520 were identified across the time series, where one GST gene was up-regulated in RRIM600 and 521 three were up-regulated in GT1 before cold stress treatment ( Figure 6). After 90 minutes, we 522 detected five genes with high similarity to GST genes. Among these five genes, four were up-523 regulated in GT1 and one was up-regulated in RRIM600. At the subsequent time point, we 524 observed an increase in up-regulated GST genes in RRIM600. After 24 hours, RRIM600 up-525 regulated nine putative GST genes, whereas GT1 had only four putative GST genes ( Figure 6, 526 Supplementary and are hypersensitive to cold and salt stress. 550 The KEGG annotation revealed up-regulated genes in the MAPK pathways for both genotypes 551 across the time series. GT1 and RRIM600 exhibited up-regulated MKK2 and MPK4/6 putative 552 genes before the initiation of the chilling treatment. However, pair-wise comparison between 553 RRIM600 and GT1 revealed that RRIM600 exhibited one MEKK1/ANP1 (EC 2.7.11.25) gene 554 that was exclusively uo-regulated after 90 minutes of cold stress. 555 Regarding MAPK signaling pathway, RRIM600 showed up-regulation of the basic endochitinase 556 B (CHI-B) (EC 3.2.1.14) gene at 12 h and 24 h. This gene product is involved in the 557 ethylene/jasmonic acid signaling pathway. Ethylene/jasmonic acid have also been shown to 558 neutralize chilling stress, activate ROS scavenging enzymes, and regulate the C-repeat binding 559 factor (CBF) pathway during cold stress (Sharma and Laxmi, 2016). In this study, the number of 560 genes exhibiting high similarity to CHI-B increased from 1 to 4 after 24 h. In this study, we detected a significant increase in up-regulated LecRLK genes in RRIM600 576 relative to GT1 after 24 h of cold treatment. Before initiating the cold treatment, the LecRLK 577 gene was not up-regulated in RRIM600; however, after 24 h, we detected eight up-regulated 578 genes. Interestingly, in GT1, the number of genes identified during cold treatment decreased 579 from seven to one gene ( Figure 7A, Supplementary CaM genes were detected in both genotypes ( Figure 7C). Relative to RRIM600, GT1 showed up-592 regulation of eight CaM genes prior to cold stress, while after 90 minutes, we identified a total of 593 nine CaM genes ( Figure 7C). We identified two up-regulated CaM putative genes in GT1 after 594 12 h and 24 h of cold stress. However, RRIM600 showed only one up-regulated CaM gene 595 before cold treatment relative to GT1 ( Figure 7C). This gene was up-regulated in RRIM600 596 during the entire treatment period. Interestingly, we also identified one up-regulated CIPK gene 597 in GT1 at 90 minutes (Supplementary Table 2). Calcium/calmodulin-regulated receptor-like 598 kinase 1 (CRLK1) is required for cold tolerance via the activation of MAP kinase activity (Yang 599 et al., 2010). This gene activates MEKK1 in response to cold in a calcium-dependent manner 600 (Furuya et al., 2013). The differential expression analysis between the two genotypes identified 601 one CRLK1 in GT1 after 12 h of cold stress. Additionally, we identified one up-regulated 602 CDPK2 gene in RRIM600 after 24 h of treatment. 603 Cysteine-rich receptor-like kinases (CRKs) can significantly affect plant development and stress 604 responses . It has been suggested that CRK transcript levels are elevated in 605 response to salicylic acid (SA), pathogens and drought. Additionally, CRKs are involved in 606 mediating the effects of ROS (Bourdais et al., 2015). It was observed that the number of up-607 regulated CRK genes increased in both genotypes due to cold treatment. Before chilling stress, 608 the comparison between GT1 and RRIM600 exhibited four and three up-regulated CRK genes, 609 respectively. However, after 24 h cold treatment, GT1 and RRIM600 exhibited five and four up-610 regulated CRK genes, respectively. Interestingly, all of the up-regulated genes identified in GT1 611 after 24 h are different from the previous up-regulated CRK genes ( Figure 7B, Supplementary 612 Table 2). 613 Proline-rich extensin-like receptor protein kinase (PERK) is structurally organized into a proline-614 rich N-terminal domain, followed by a transmembrane segment and a C-terminal 615 serine/threonine kinase domain (Florentino et al., 2006). The PERK genes are assigned to two 616 major groups. One group is expressed exclusively in pollen, and the other is expressed 617 throughout all tissues. PERKs may be involved in the early response to general perception and 618 the response to wounding and/or pathogen stimuli and turgor pressure (Osakabe et al., 2014). 619 GT1 0h and RRIM600 0h exhibited two and three up-regulated genes, respectively. During cold 620 treatment, the number of up-regulated genes increased in both genotypes. After 24 h of chilling 621 stress, the comparison between RRIM600 and GT1 revealed five and seven up-regulated genes 622 in each genotype, respectively (Supplementary Table 2). 623

Photosynthesis Activity and Stomata Closure 624
Cold stress can cause an imbalance between light utilization and energy dissipation through 625 metabolic activity. An excess of photosystem II (PSII) excitation pressure exists, which can be 626 reversible through the dissipation of excess absorbed energy or irreversible inactivation of PSII 627 and, consequently, inhibition of the photosynthetic capacity (Oquist and Huner, 2003). The 628 imbalance in the PSII caused by cold stress might generate ROS, which can then damage the 629 photosynthetic apparatus and the whole cell (Tyystjarvi, 2013). Therefore, tolerance to cold-630 induced photoinhibition can be considered a mechanism for cold tolerance. Additionally, 631 RubisCO plays a central role in CO 2 assimilation and photosynthesis efficiency. Crops that are 632 acclimated to cold, such as winter wheat and rye, adjust their RubisCO content and are able to 633 maintain a high CO2 assimilation rate (Yamori et al., 2010). 634 Low temperature can also affect the enzymes and ion channels responsible for maintaining the 635 guard cell osmotic potential (Ilan et al., 1995). Due to the reduction in the photosynthetic 636 capacity caused by cold stress, there is an increase in internal CO 2 in the substotamal cavity, 637 which reduces the stomatal aperture (Wilkinson et al., 2001). 638 In this study, GO enrichment analysis revealed enriched categories in up-regulated genes in GT1 639 related to photosynthesis compared to RRIM600 after 24 cold stress. Although we observed 640 enriched GO categories related to plant defense in GT1 after 24h, we also identified enriched 641 categories associated with photosystem II ( Table 3). Additionally, the qPCR validation 644 corroborated the in silico DEG analysis, in which the RubisCO gene was found to be up-645 regulated in GT1. 646 The pair-wise comparison between GT1 and RRIM600 revealed that the abscisic acid (ABA) 647 receptor PYL4 gene was up-regulated in RRIM600 at 0 h, 90 minutes and 24 h of chilling stress. 648 Additionally, after 24 h of cold stress, RRIM600 showed up-regulation of two abscisic acid 649 receptor PYL4 genes and one ABA receptor PYR1 gene (Supplementary Table 3 pressure on the PSII reaction diverting photon energy into heat via zeaxanthin (Sui et al., 2007). 661 Furthermore, we identified one up-regulated gene in RRIM600 after 12 and 24 h of chilling 662 stress with high similarity to the HT1 gene. In Arabidopsis, the HT1 serine/threonine kinase gene 663 is reported to be involved in the control of stomatal movement in response to CO 2 . Additionally, 664 genes with high similarity to hexokinase-1 and PtdIns3P 5-kinase (PI3P5K) were up-regulated in 665 RRIM600 at 24 h of treatment. These genes are also reported to be related to stomatal closure via 666 ABA (Bak et al., 2013) (Supplementary Table 2). 667 In Arabidopsis, SRK2E is a positive regulator in ABA-induced stomatal closure and is involved 668 in stress adaptation. Additionally, SRK2E acts as a transcriptional repressor involved in the 669 inhibition of plant growth under abiotic stress conditions (Yoshida et al., 2006). reported to be involved in stomatal closure were up-regulated in RRIM600. We also observed 735 that genes related to cell growth were up-regulated in GT1 after 24 h hours of cold stress. 736 Although we identified genes related to the defense response in GT1, the DEG and GO 737 enrichment analyses showed that GT1 remains active, displaying up-regulation of genes related 738 to photosynthesis during cold treatment. In addition, GT1 probably has a more efficient strategy 739 for strengthening the cell wall, as revealed through GO enrichment analysis. 740 Rubber tree breeding programs are interested in genotypes resistant to cold and that exhibit latex 741 production that does not cease during the cold season. The elucidation of different chilling 742 tolerance strategies linked to information about possible genes involved in such responses, 743 including the identification of molecular markers in these genes associated with information on 744 AS events, provides a powerful tool for the genetic and genomic analysis of the rubber tree for 745 breeding strategies and future studies involving GWAS and GS. 746

Conflict of Interest 747
The authors have no conflicts of interest to declare. 748