BLAST and Hidden Markov Model searches were conducted to identify H. brasiliensis CIPKs and CBLs using Arabidopsis, rice and poplar CIPK and CBL protein sequences as queries. A total of twelve CBLs (named HbCBL1 to 12) and thirty CIPKs (named HbCIPK1 to 30) were identified from the H. brasiliensis genome (Tang C. et al., 2016). Detailed information, including protein length, isoelectric point (pI), molecular weight (MW), number of introns and evolution group of identified HbCBLs and HbCIPKs, is listed in Table 1. The number of amino acid residues of the identified HbCIPKs ranged from 328 (HbCIPK22) to 541 (HbCIPK17), and HbCBLs ranged from 165 (HbCBL2) to 247 (HbCBL5). Their relative molecular mass ranged from 37.29 kDa (HbCIPK22) to 61.30 kDa (HbCIPK17) for HbCIPKs and 19.38 kDa (HbCBL2) to 28.44 kDa (HbCBL5) for HbCBLs. The isoelectric points of HbCBLs and HbCIPKs were between 4.64 (HbCBL6) and 9.21 (HbCBL2), and between 5.47 (HbCIPK28) and 9.51 (HbCIPK22), respectively. It is noteworthy here that the isoelectric point of HbCBL2, at 9.21, was an outlier that was significantly higher than those of the other HbCBLs (Table 1). Conserved domain analysis confirmed that most of the HbCIPKs identified harbored the Pkinase (PF00069) and NAF (PF03822) domains, while all the HbCBLs harbored the calcium-binding EF hand domains (PF13833 and PF13499) which are the hallmark of CBL family.

Phylogenetic analysis of the thirty HbCIPKs, together with 103 CIPKs from five other plants, classified the CIPK family into five clusters (Group A to E, Figure 1A and Supplementary Tables 1, 2). There were eleven HbCIPK members (HbCIPK1, 6–11, 16–17, 27 and 30) in Group A, thirteen (HbCIPK2, 4–5, 12, 15, 18, 22–26, 28–29) in Group C, and two each in the other three groups, i.e., HbCIPK13 and 21 in Group B, HbCIPK14 and 19 in Group D, and HbCIPK3 and 20 in Group E. As expected, CIPKs from Hevea generally had closer relationships with those from the two Euphorbiaceae family species, Manihot esculenta and Ricinus communis, than those from Arabidopsis, rice or poplar, in accordance with plant taxonomic status. Phylogenetic analysis based on predicted amino acid sequences identified some closely related pairs of HbCIPKs, such as HbCIPK9 and 11 in Group A, HbCIPK13 and 21 in Group B, HbCIPK2 and 23, HbCIPK4 and 5, HbCIPK18 and 24, and HbCIPK26, 28 and 29 in Group C, and HbCIPK3 and 20 in Group E (Figure 1A). Similar trends are also observed in other plant species (Hu et al., 2015), and closely related HbCIPK pairs often have similar cellular localization and functions. The exon-intron structures of the thirty HbCIPK genes were determined based on their predicted genomic sequences. As shown in Figure 2A, most HbCIPK members within the same groups shared similar gene structure in terms of intron number, domain localization and exon length. It is worth noting that the members in Group A had eleven to sixteen introns, while those in the four other groups had no intron or only one. A similar intron-rich or poor pattern is also observed in the CIPKs from other plant species, such as Arabidopsis, rice, poplar, cassava, pear, and soybean (Kolukisaoglu et al., 2004; Zhang et al., 2008; Hu et al., 2015; Tang J. et al., 2016; Zhu et al., 2016), reflecting a conserved feature of the CIPK family in gene structure. Considering the rate of intron loss is faster than that of intron gain after segmental duplication (Roy and Penny, 2007), the CIPKs of group A might represent the progenitor genes of the CIPK family.

Phylogenetic analysis put the CBL genes of Hevea and the other five plants into four groups (Group I to IV, Figure 1B), with three, two, two and five members in Group I, II, III, and IV, respectively. This analysis also identified six closely related orthologous pairs between Hevea and cassava (HbCBL9 and MeCBL8, HbCBL7 and MeCBL7, HbCBL2 and MeCBL5, HbCBL3 and MeCBL9, HbCBL8 and MeCBL1, HbCBL10 and MeCBL6), suggesting the existence of an ancestral set of CBL genes prior to the divergence of Hevea and cassava. The results of gene expression analysis (Figure 3) showed that HbCBL1, 8 and 10 might be the main direction of HbCBLs evolution. Our analysis of genomic organization of the HbCBLs revealed that eight of the twelve members contained seven introns, two contained eight introns, one contained four introns and one contained six introns (Figure 2B). Except for HbCBL2 and HbCBL8, most HbCBL members within the same groups shared very similar gene structures in terms of intron number, domain localization, and exon length.

The locations of HbCIPK and HbCBL genes in Hevea chromosomes were as shown in Figure 4A. The HbCIPK genes were mapped onto fourteen of the eighteen Hevea chromosomes, while the HbCBL genes were located on seven chromosomes. However, distribution of these two sets of genes on the chromosomes was uneven. There was no HbCIPK or HbCBL members mapped onto chromosomes 6, 11, 13, and 16, whereas the other chromosomes contained one to five members. Chromosome location analysis also found some tandem repeat genes, such as HbCIPK19 and 26, HbCIPK28 and 29, HbCIPK14 and 15, HbCIPK3, 4 and 5, which were located on chromosomes 2, 7, 12, and 17, respectively. The same was true of HbCBL1 and 11, HbCBL5 and 6, HbCBL3 and 4 located on chromosome 5, 8, and 10, respectively. Collinearity analysis found some HbCIPK members might have been amplified by chromosome duplication, such as in the case of HbCIPK26 and 28, HbCIPK24 and 18, HbCIPK11 and 9, HbCIPK13 and 21, HbCIPK6 and 7, HbCIPK20 and 3, HbCBL5 and 8, HbCBL10 and 12, HbCBL4 and 9, HbCBL2 and 7. This was consistent with their clustering patterns shown in the phylogenetic tree (Figures 1A, 2). Chromosomal location and colinearity analysis of CBL and CIPK genes within Hevea (Figure 4A) and across Hevea and three other plant species (Figure 4B) were investigated in order to explore the potential evolutionary relationships. The results revealed a higher homology between H. brasiliensis, M. esculenta and P. trichocarpa than that between P. trichocarpa and A. thaliana. Some HbCIPK and HbCBL genes were collinear with CIPK and CBL genes in P. trichocarpa, A. thaliana and M. esculenta, suggesting their important roles in plant evolution. These results can be useful for subsequent comparative studies of CIPK and CBL genes with known functions.

To study the expression profiles and deduce the gene function of HbCBL-HbCIPKs, transcriptome analyses were conducted in seven Hevea tissues: latex, bark, leaf, root, seed, and female and male flowers using Solexa sequencing and quantitative PCR as described previously (Xiao et al., 2017). As shown in Figure 3, transcripts of thirteen HbCIPK genes (HbCIPK1, 3–6, 8, 12, 17, 19–20, 22, and 28–29) and six HbCBL genes (HbCBL2, 4–6 and 12) were barely detectable in almost all the tissues examined, indicating the loss of their functions during evolution. In contrast, HbCIPK14, 23, 30 and HbCBL1 were expressed in most of the tissues examined, indicating their functional conservation in the course of evolution. Some of the HbCIPK and HbCBL genes displayed preference expression in specific tissues, such as HbCIPK15, 18 and 26 and HbCBL1 and 3 in latex, HbCIPK13-14 and 23 and HbCBL1 in the root, HbCIPK14 and 30 and HbCBL10 in the bark, HbCIPK9, 14 and 30 and HbCBL8 in the leaf, and HbCIPK13-15 and 21 and HbCBL1 in the seed. It is interesting to note that some of the HbCIPK and HbCBL members showed distinct expression patterns in female and male flowers. For example, the expression levels of HbCIPK2, 9, and 30 and HbCBL10 were much higher in the female flowers than in the male flowers. On the other hand, HbCIPK10-11, 16, 23 and 27 and HbCBL1 and 8 were much higher expressed in the male flowers than the female flowers. The expressions of randomly selected HbCIPK and HbCBL genes including HbCBL3 and 8 and HbCIPK14, 16, 23 and 27 were further investigated by qPCR (Supplementary Figure 1), revealing patterns basically consistent with transcriptome analysis.

To obtain information on the functions of HbCIPK and HbCBL genes in the course of leaf development, we examined their expression levels by RNA-seq at four progressive leaf stages (bronze, color change, pale-green, and mature). As shown in Figure 5, some of the HbCIPK genes, i.e., HbCIPK13-14, 16, 20–21, 25–26, and 30, were obviously up-regulated during leaf development. Of the HbCBL genes, only HbCBL7 was similarly up-regulated with leaf development. Some HbCIPK and HbCBL genes were up-regulated during the first three leaf stages (bronze, color change, and pale-green) but down-regulated toward leaf maturity, as seen for HbCIPK2, 10, 19, 23, 27 and HbCBL8. A small number of HbCIPK and HbCBL genes were down-regulated during leaf development, such as HbCIPK11 and HbCBL1 and 10. Some closely related gene pairs exhibited similar expression profiles. For example, HbCIPK2 and 23, HbCIPK28 and 29, HbCIPK13 and 21, HbCIPK9 and 11 had similar expression patterns during leaf development from soft young leaves to fully hardened mature leaves. Nevertheless, there are more pairs of closely related genes showed dissimilar expression patterns. Examples of these were HbCIPK1 and 10, HbCIPK3 and 20, HbCIPK8 and 16, HbCIPK14 and 19, HbCIPK17 and 27, HbCBL1 and 11, HbCBL2 and 7, HbCBL10 and 12. The expressions of HbCIPK14-16, 23, 26 and 27 were further investigated by qPCR (Supplementary Figure 2), all of which were obviously up-regulated during leaf development, consistent with the results of transcriptome analysis.

Ethephon is an ethylene generator which is widely used in Hevea to stimulate latex flow, and hence rubber yield. The mechanisms of ethylene signal transduction and yield stimulation are, as yet, poorly understood. Since the CBL-CIPK complex plays an important role in calcium signaling in relation to various plant responses (Luan, 2009), we examined the expression levels of HbCBL and HbCIPK genes in latex upon ethephon treatment. As shown in Figure 6, HbCIPK14-16 and 27 and HbCBL10 were obviously up-regulated following the ethephon treatment, while HbCIPK7 and HbCBL3 were down-regulated. Expressions of HbCIPK11 and 30 displayed a transient high expression at 12 h after ethephon treatment but weakly expressed at 24 h. The expressions of HbCBL1, 3 and 8 and HbCIPK10, 14–16, 18 and 27 were further investigated by qPCR (Supplementary Figure 3); all showing patterns consistent with those of the transcriptome analysis. The above results implicated the HbCIPK-HbCBL complexes in ethylene-simulated latex production.

To assess the respective effect of latex tapping on the expression of HbCIPK and HbCBL genes, experiments were performed on previously untapped cultivated trees of the clones PR107 and ReYan8-79, and a number of HbCIPK and HbCBL genes with substantial expression in latex were examined. The first tapping of an untapped rubber tree normally produces very little latex. The latex yield increases progressively with successive tappings at regular intervals to reach a steady output after seven to ten tappings (Tang C. et al., 2010). As shown in Figures 7A,B, the expressions of HbCIPK27 and HbCBL3 were clearly up-regulated when untapped PR107 and ReYan8-79 trees were first brought into tapping. On the other hand, HbCIPK16, 18, 23, 30 and HbCBL1 were significantly down-regulated from the second tapping onward. The above results suggested that the CIPK-CBL complexes might play a role in tapping-stimulated latex production. The expression patterns of HbCIPK and HbCBL genes in the two Hevea clones were consistent in most cases, indicating the ways of their response to the tapping treatment being similar under different Hevea genetic backgrounds.

Calcineurin B-like proteins and calcineurin B-like interacting protein kinases are functionally linked through physical interactions, and participate in plant development various environmental stresses (Luan, 2009). The interaction relationships between the HbCBL and HbCIPK proteins with substantial gene expression in latex or other Hevea tissues were detected by yeast two-hybrid (Y2H), i.e., spotting co-transformed yeast cells on various types of selective media. As shown in Figure 8A and Supplementary Figure 4, multiple combinations of HbCBL and HbCIPK fusion proteins showed obvious interactions although the strength of interaction varied. HbCBL1 and its closely related homolog, HbCBL10, displayed overlapping strong or moderate interactions with more than half of the HbCIPKs examined. In contrast, HbCBL3 and its close homolog HbCBL8 did not interact with most of the HbCIPKs investigated. Similar overlapping patterns of interaction were observed for the closed related HbCIPK homologs, e.g., HbCIPK11 and 18, and HbCIPK14 and 15. Expression profiles of the genes encoding HbCBL1 and its four strong interaction HbCIPK partners revealed by the Y2H assay (Figure 8A), HbCIPK14, HbCIPK15, HbCIPK16, and HbCIPK26, were examined in different tissues and treatments based on Solexa sequencing (Figure 8B and Supplementary Table 4). The expressions of HbCIPK26 and HbCBL1 displayed a strong correlation in different tissues and treatments, particularly in latex. To validate the results of the Y2H assay, the interactions of HbCBL1-HbCIPK15 (strong interaction) and HbCBL1-HbCIPK30 (weak interaction) were further investigated by the BiFC experiments conducted in N. benthamiana leaves. As shown in Figure 9, the results were consistent with those of the Y2H analysis, with the co-transformation of the HbCBL1-HbCIPK15 pair yielding a stronger green fluorescence than that of the HbCBL1-HbCIPK30 pair.

Sequences of Arabidopsis thaliana, Populus trichocarpa, and Oryza sativa CBL and CIPK genes were downloaded from GenBank¹. The genome and protein sequences of Arabidopsis thaliana, Oryza sativa, Populus trichocarpa, Manihot esculenta, and Ricinus communis were downloaded from Phytozome². The Hevea Genome and transcriptome were obtained from GenBank³. Local BLAST and Hidden Markov Model searches were conducted to identify Hevea CBLs and CIPKs using Arabidopsis, rice and poplar CBL and CIPK protein sequences as queries to search against the proteome of each species for the candidate CIPKs and CBLs from H. brasiliensis, Arabidopsis thaliana, Oryza sativa, Populus trichocarpa, Manihot esculenta, and Ricinus communis. All putative candidates were manually verified with the InterProScan server⁴ to confirm the presence of relevant protein domains. The website ProtParam⁵ was used to predict protein isoelectric point and molecular weight.

Multiple alignments of the amino acid sequences of CIPK and CBL from Hevea brasiliensis and five other species were set up and phylogenetic trees were constructed with MEGA6.0 by employing the Neighbor-Joining (NJ) method with a bootstrap test for 1,000 replicates. Exon/intron structures of HbCBL and HbCIPK genes were analyzed by comparing the cDNA and their genomic DNA sequences, and chromosomal locations analyzed through mapping the genes to chromosomes, both using the TBtools software (Chen et al., 2020). The collinearity analysis was completed using MCscanX software (Wang et al., 2012).

The Solexa sequencing data for various tissues/organs, developmental stages and ethephon stimulation available at the NCBI Sequence Read Archive (SRA) database were used for the expression analysis of Hevea brasiliensis as described previously (Project Accession: PRJNA310171, Supplementary Table 3) (Xiao et al., 2017). The Hevea tissues included (a) latex, bark, leaf, root, seed, female flower and male flower, (b) leaves of four developmental stages (bronze, color-change, pale-green and mature), and (c) latex samples collected at 0, 3, 12, and 24 h after ethephon stimulation. Raw RNA-seq reads were processed to trim terminal low-quality bases and adapter sequences via an in-house custom pipeline. The clean reads were then mapped to the Hevea brasiliensis genome using Bowtie2, and RSEM software was used for quantifying transcript abundance with default parameters (Li and Dewey, 2011). The normalized value expression profiles of HbCBL1 and its interaction partner genes were visualized in bar-plots by R⁶.

To verify the results obtained by Solexa sequencing and to examine genes expression, qPCR was performed on a number of the Hevea brasiliensis CIPK and CBL genes as described previously (Li et al., 2011). The primer pairs used for the CIPK and CBL genes were listed in Supplementary Table 1. For internal control, the genes RH2b and YLS8 were used as described previously (Li et al., 2011).

To analyze the relationship between CBL and CIPK family members, the expression correlations among family members were analyzed. The input data of correlation analysis was the FPKM values of RNA-seq, included the Hevea (a) tissues latex, bark, leaf, root, seed, female flower and male flower, (b) leaves of four developmental stages (bronze, color-change, pale-green and mature), and (c) latex samples collected at 0, 3, 12, and 24 h after ethephon stimulation. R package Corrplot⁷ was used for analysis and drawing.

The full length cDNA of eight HbCIPK and four HbCBL genes with substantial expressions in latex were amplified by PCR with their respective primers (Supplementary Table 1) and inserted into the Y2H vectors of pGADT7 and pGBKT7, respectively. The pGBKT7-HbCBLs and pGADT7- HbCIPKs vectors were co-transformed into the Y2HGold yeast (Saccharomyces cerevisiae) strain and examined for their interactions using the MatchMaker yeast two-hybrid system (Clontech, United States). The detailed protocols were as described (Kleist et al., 2014). To test the functionality of the latex expressed HbCBL-HbCIPK complexes in salt tolerance SOS pathway, complementation test was conducted using the yeast mutant strain AXT3K (4ena1:HIS3:4ena4,4nha1:LEU2, and 4nhx1:KanMX4) that lacks the main plasma membrane Na⁺ transporters (Quintero et al., 2011; Zhou et al., 2015). The full length coding regions of HbCBL1/8 and HbCIPK14/15/27/30 were amplified by PCR (Supplementary Table 1) and cloned into the yeast expression vector p414 (Mo et al., 2018), whereas the coding regions of HbSOS1 and AtSOS1 were cloned into the yeast expression vector p416 (Yin et al., 2020). The three plasmids (p414-HbCBLs, p414-HbCIPKs, and p416-HbSOS1 or AtSOS1) were co-transformed into the yeast strain, and the salt tolerance tests were performed according to the detailed process as described previously (Yin et al., 2020).

Open reading frames (ORFs) of HbCIPK15, HbCIPK30 and HbCBL1 were amplified by RT-PCR and cloned into the pEG100-YFP vectors to construct fusion proteins at the N- and C-termini. Primers are listed in Supplementary Table 1. Each cDNA was under the control of the 35S promoter. For transient expression in Nicotiana benthamiana leaves, plasmids were introduced into Agrobacterium tumefaciens strain GV3101 competent cells for infiltration of 5-week-old N. benthamiana leaves. For microscopic observation, the reconstructed green fluorescence protein (GFP) signals of the lower epidermal cells of leaves cut 4 day after infiltration were examined using a Nikon A1RHD25 confocal microscope (Nikon, Japan).