To identify MdCHK receptors, BLAST searches were performed against the Genome Database for Rosaceae (GDR; Jung et al., 2014) using A. thaliana cytokinin receptor sequences as queries (AHK2, AHK3, and AHK4). Five sequences were identified based on genome analysis and corresponding cDNA were amplified from various plant organs using specific primers (Supplementary Table S1). Sequences were registered in Genbank as MdCHK2 (KM114879), MdCHK3a (KM114880), MdCHK3b (KM114881), MdCHK4a (KM114883) and MdCHK4b (KM114882).

MdCHKs gene organization has been visualized using the FancyGene program (Rambaldi and Ciccarelli, 2009). Phylogeny analyses were performed on conserved domains and local similarities among proteins sequences. To this aim, multiple protein sequence alignments were done using the COBALT tool (Papadopoulos and Agarwala, 2007) and sequences were curated with Gblocks prior the construction of a bootstrap neighbor joining tree. Protein domain predictions were acquired using the SMART (Letunic et al., 2015) and PROSITE (Sigrist et al., 2002) programs, and transmembrane regions were identified with TMHMM (Krogh et al., 2001) and TMpred tools (Hofmann and Stoffel, 1993). Visualization of the transmembrane helixes has been performed with a helical wheel drawing program¹.

The S. cerevisiae strain YIL147C, deficient in SLN1 receptor (MATa/α, ura3, leu2, his3, can1Δ::LEU2-MFA1pro-HIS3/CAN1, sln1Δ::KanMX/SLN1) was used in complementation assays. Full-length coding sequences of MdCHKs were amplified and cloned into the yeast expression vector pYES2 under the control of the GAL1gene promoter using the NotI restriction site (for primers, see Supplementary Table S1). The S. Cerevisiae strain was transformed as follows. Cells were grown in 200 mL YPD liquid medium (150 rpm, 28°C) to 0.4–0.6 OD₆₀₀, harvested by centrifugation (3000 g, 10 min) and resuspended in 0.1 M lithium acetate, 10 mMTris–HCl, pH 7.5,1 mM EDTA, 10 mM DTT. After 1 h incubation at 30°C, cells were washed twice with ice-cold 1 M sorbitol and resuspended in 1–5 mL ice-cold 1 M sorbitol. Plasmid DNA (0.5–1 μg DNA) was added to a 180 μL cell suspension and transferred to a 0.2 cm gap width electroporation cuvette. Electroporation was performed using a Bio-Rad gene pulser with an electric pulse of 2.5 kV, 25 μF and 200 Ω. Cells were immediately washed out from the cuvettes after the electroporation with YPD, plated on selective CSM-URA medium and incubated at 28°C for 3–5 days. Fresh colonies were then grown 3 h in liquid YPD at 30°C and meiosis was induced by pouring the suspension cell on ACK medium (10 g/L potassium acetate, 2.5 g/L yeast extract, 20 g/L agar). These plates were incubated 1 week at 20°C and haploids were finally selected on MMAS medium plates [20 g/L galactose, 7 g/L YNB (WA), 0.6 g/L DOB-LEU-HIS-ARG-URA, 0.06 g/L L-canavanine, 0.2 g/L G418, 20 g/L agar] supplemented with 10 μM trans-zeatin at 28°C for 3–5 days. Suspensions of transformants were then spotted onto dropout media containing or not 10 μM of trans-zeatin with 2% galactose and grown for 48 h at 28°C. For specificity and sensitivity assays, complemented-yeast growth was carried out in liquid YCGal (7 g/L YNB, 0.8 g/L CSM-URA, 20 g/L galactose) supplemented with various types and concentrations of cytokinins for 48 h at 28°C. Cell growth was measured at 630 nm (BioHit Reader BP800).

Pure standards of isopentenyladenosine 5′-monophosphate, trans-zeatin riboside 5′-monophosphate, cis-zeatin riboside 5′- monophosphate, isopentenyladenine, trans-zeatin, cis-zeatin, dihydrozeatin, isopentenyladenosine, trans-zeatin riboside, cis-zeatin riboside, dihydrozeatin riboside, dihydrozeatin riboside 5′-monophosphate, trans-zeatin O-glucoside, trans-zeatin O-glucoside riboside, trans-zeatin N7-glucoside, 2- methylthio-isopentyladenine, 2-methylthio-trans-zeatin, 2-methylthio-cis-zeatin, 2-methylthio-isopentyladenosine, 2-methylthio-trans-zeatin riboside, 2-methylthio-cis-zeatin riboside were purchased from Olchemim (Olomouc, Czechia).

Extraction of total RNA from M. domestica organs was performed using the NucleoSpin RNA extraction kit (Marcherey-Nagel), with improved lysis step (McKenzie et al., 1997). First-strand cDNA were synthesized from 1 μg of total RNA using the iScript cDNA Synthesis Kit (Bio-Rad). Quantitative real-time PCR measurements were carried out in triplicate using SsoAdvanced Universal SYBR Green (Bio-Rad) in a 15 μL final volume containing 6 μL diluted template cDNA and specific primers (0.5 μM) (Supplementary Table S1). Amplification was performed on a CFX96 Touch real-time PCR system (Bio-Rad) with the following conditions: 95°C for 7 min and 40 cycles at 95°C for 10 s and 60°C for 40 s. Amplification was followed by a melt curve analysis. Absolute quantification of transcript copy number was assessed with calibration curves. Transcript levels were then normalized with EF1α.

Subcellular localization of MdCHK receptors were studied in Catharanthus roseus C20D cells transiently transformed using plasmid-coated particles bombardment as described in Guirimand et al. (2009). The full length MdCHK sequences were amplified and cloned into the SpeI restriction site of pSCA-YFP plasmid (for primers, see Supplementary Table S1), in frame with the 5′ extremity of the YFP coding sequence. The endoplasmic reticulum (ER) cyan fluorescent protein (CFP) marker (Guirimand et al., 2010) was used in co-transformation assays.

Dynamic localization of MdCHK receptors was also studied in yeast S. cerevisiae strain WT303 (MATa/α, leu2, trp1, ura3, ade2, his3) transformed by pESC-LEU plasmids (Foureau et al., 2016) containing MdCHK sequences, except for MdCHK2, cloned in pYES2 and transformed in sln1Δ yeast strain to bypass the sequence toxicity in microorganisms. The CYP450 T16H2 sequence was clone in the pESC-TRP plasmid in fusion with the 5′end of the CFP sequence and used as an ER marker (Besseau et al., 2013). MdCHK sequences were cloned under galactose inducible promoter and fused at the 5′ end with the YFP sequence. Transformed colonies were cultivated on selective plates (CSM-LEU or CSM-URA, supplemented by 2% glucose, respectively, for pESC-LEU and pYES2) at 30°C for 48 h and then transferred in inducing liquid media (CSM-LEU or CSM-URA, supplemented by 2% galactose, respectively, for pESC-LEU and pYES2) with or without iP (5 μM) for additional overnight culture at 28°C.

An Olympus BX51 epifluorescence microscope equipped with the Olympus DP71 digital camera and Cell^∗D imaging software (Soft Imaging System Olympus) was used for image capture and merging false-colored images of both C20D cells and S. cerevisiae colonies expressing YFP.

Bimolecular Fluorescent Complementation (BiFC) experiments were conducted using SPYNE and SPYCE plasmids (Waadt and Kudla, 2008). MdCHK sequences were amplified and cloned into the SpeI restriction site (for primers, see Supplementary Table S1), in frame with the 5′ extremity of a truncated YFP coding sequence. Transient transformation of C. roseus cells by particle bombardment and YFP imaging were performed according to Guirimand et al. (2009) with adaptation for BiFC assays (Guirimand et al., 2010). Interactions were tested in triplicates using three independent plasmid clones.

Available RNAseq data for M. domestica was downloaded from NCBI via SRA toolkit2.6.2. The recovered SRA files were transformed in fastq format with the “fastq-dump” command from SRA toolkit. The files were cleaned with Trimmomatic 0.36 with default parameters and using provided adapter sequences for TruSeq2 and TruSeq3. The transcription quantification was performed with Salmon 0.6.1 using the Variational Bayesian EM algorithm and biase correction. TPM (transcripts per million) values from the resulting quant.sf files were combined under R 3.3.0 in an expression matrix containing 95,232 predicted genes (Malus domestica v3.0) × 250 experimental conditions. Using the expression matrix, Pearson Correlation Coefficients (PCC) and further Highest Reciprocal Ranks (HRR) computation were performed using a homemade program written in C [HRR (gene A, gene B)] = max [rank (gene A, gene B), rank (gene B, gene A)] to establish the co-expressed genes lists for MdCHK and GO enrichment tests. For each MdCHK highly co-expressed genes, i.e., genes with a HRR ≤ 500 were selected. The procedure was repeated with publicly available A. thaliana RNAseq data. We similarly prepared an expression matrix containing 33,604 transcripts (Arabidopsis TAIR v10 genome annotation) and 1,676 samples. Co-expressed genes lists were obtained after calculating PCC and ranking them with HRR. For each AtCHK (AT5G35750.1, AT1G27320.1, AT2G01830.2), gene pairs having an HRR < 500 were considered to be significantly co-expressed. Orthology between Arabidopsis and apple tree was obtained from Plaza 3.0 (Van Bel et al., 2012). The functions represented by coexpressed genes of each MdCHK were analyzed with the Gene Ontology classification. A BlastX was performed on the M. domestica genome v3.0 to recover correspondent protein sequences and Pfam domains were identified using Hmmer. The functional annotations of the M. domestica genome v3.0 were generated using Trinotate on the previous data. In order to determine potential functional enrichment for every target gene, enrichment of GO terms was tested by comparing effectives to a hypergeometric distribution (p-value cut-off = 0.001) using the R “phyper” function. To compare redundancies in the five co-expressed genes lists, a Venn diagram was drawn using the venneuler package 1.1 for R.

Based on the CHASE domain of the three A. thaliana CHKs, apple tree genome was browsed to identify putative CHK sequences in M. domestica. Five candidates were identified and named MdCHK2 (locus tag MDP0000258078), MdCHK3a (locus tag MDP0000310800), MdCHK3b (locus tag MDP0000155347), MdCHK4a (locus tag MDP0000151825) and MdCHK4b (locus tag MDP0000242242) according to their distribution within the three classical groups homologous toAHK2, AHK3, and AHK4, as shown by phylogenetic analysis (Figure 1). Genomic sequences revealed that MdCHK2 possesses 13 exons and is located on chromosome 9 (Supplementary Figure S1). Both MdCHK3a and MdCHK3b, respectively, located on chromosomes 16 and 13, display 10 exons with a similar organization regarding intron positions and intron/exon sizes. The same observation was made for MdCHK4a and MdCHK4b respectively located on chromosomes 13 and 10 and containing 11 exons. Such similarity may reflect the gene duplication events leading to the couples of CHKs homologous MdCHK3a/MdCHK3b and MdCHK4a/MdCHK4b which share respectively 94.67 and 95.67% nucleotide identity. The full-length cDNAs of the five MdCHKs were cloned and deposited at NCBI under the GenBank accession numbers KM114879 to KM114883. They contained large open-reading frames ranging from 3027 to 3612 bp encoding proteins of 1008 to 1203 amino acids.

The computational analysis of the protein sequences of MdCHKs revealed the presence of the four basic conserved domains called Cyclases/Histidine kinases Associated SEnsory (CHASE), Histidine Kinase (HK), Receiver (REC), and Receiver-like (REC-like) (Figure 2 and Supplementary Figure S2). While these modular MdCHKs share a similar multidomain architecture in their cytoplasmic C-terminal part (HK, REC-like and REC domains), they differ in the TM domain topology of the N-terminus part. Indeed, the MdCHK4a and MdCHK4b have two predicted α-helices transmembrane domains TM1 and TM2 (Supplementary Figure S3) bordering the CHASE domain (Figure 2). MdCHK3a and MdCHK3b possess equivalent TM1/CHASE/TM2 organization with an additional predicted transmembrane helice (TM3) in N-terminus. Finally, the N-terminus of MdCHK2 includes a fourth predicted transmembrane domain (TM4). Consequently, in addition of the extracytoplasmic loop containing the sensing CHASE domain, the TM3 and TM4 domains border a supplemental extracytoplasmic loop (140 aa) absent in MdCHK3a/b and MdCHK4a/b structures (Figure 2). Regarding the cytoplasmic part, each MdCHK possesses an HK domain containing the conserved phosphorylatable histidine residue as well as the C-terminal REC domain including the conserved phospho-accepting aspartate residue (Figure 2 and Supplementary Figure S2). In addition, the five receptors contain a second receiver domain located between the HK domain and the REC domain called pseudo-receiver domain or REC-like (Figure 2). Interestingly, the putative phospho-accepting aspartate residue in the REC-like domain of MdCHK2 is conserved, suggesting that this domain may be functional in terms of phosphorelay reaction, whereas the corresponding residues in MdCHK3a/MdCHK3b and MdCHK4a/MdCHK4b are substituted with glutamate (Supplementary Figure S2).

To assess the function of the five MdCHKs as cytokinin receptors, we exploited the S. cerevisiae sln1Δ deletion mutant strain which carries a lethal mutation in the SLN1 gene encoding its unique osmosensing histidine kinase (Inoue et al., 2001). The sln1Δ yeast mutants carrying pYES2-MdCHK2, pYES2-MdCHK4a or pYES2-MdCHK4b were lethal (Figure 3A). However, the addition of trans-zeatin (tZ) in culture medium allowed recovering a normal growth (Figure 3A). By depending on the presence and perception of cytokinin to complement sln1 mutation, the three recombinant yeast strains clearly demonstrated that MdCHK2, MdCHK4a, and MdCHK4b act as cytokinin receptors in this heterologous system. Concerning the recombinant strains carrying pYES2-MdCHK3a and pYES2-MdCHK3b, they displayed an original phenotype since they exhibited a basal growth in absence of cytokinin, especially for MdCHK3b (Figure 3A and Supplementary Figure S4). However, the addition of tZ clearly induced the yeast growth pointing out that both MdCHK3a and MdCHK3b sense cytokinins (Figure 3A and Supplementary Figure S4). The basal constitutive activity of these two cytokinin receptors raised the question of their putative additional sensing function.

We further evaluated the substrate specificity and sensitivity of M. domestica cytokinin receptors. In this way, we measured the growth of the yeast cells in presence of various cytokinin-types at different concentrations including free-bases, ribosides, glucosides, and methylthio cytokinins. The specificity as well as the sensitivity (the minimal cytokinin concentration that induced yeast growth) were reported on Figure 3B and Supplementary Figure S4. Three distinct profiles corresponding to MdCHK2, MdCHK3, and MdCHK4 groups were observed. First, MdCHK2 clearly perceived a wide range of cytokinin forms since each cytokinin type activated the receptor except tZ7G. Furthermore, this receptor presented a remarkable higher sensitivity than other MdCHKs (until 0,1 nM for iP and 2MeSiP forms) and was the only one to be activated by some ribosides and glucosides cytokinin-types. Secondly, the strictly cytokinin-dependent MdCHK4a and MdCHK4b were activated by the free-bases iP and tZ and the methylthio-forms 2MeSiP, 2MeStZ, 2MeScZ, and 2MeSiPA (Supplementary Figure S4). High concentrations (10 μM) of cZ and DHZ were also effective on MdCHK4b. Nevertheless, the cytokinin-sensitivity of both MdCHK4a and MdCHK4b was obviously lower than MdCHK2. Finally, even if MdCHK3a and MdCHK3b showed a basal growth in absence of cytokinin, we were able to detect a significant difference of growth in presence of iP, tZ, DHZ, 2MeSiP, and 2MeStZ (Supplementary Figure S4).

To examine the gene expression of MdCHKs, RT-qPCR was carried out using distinct plant organs including roots, stems, leaves, flower buds and flowers. Transcripts of the five MdCHKs were detected in all the tested organs, but with distinct expression pattern. Thus, MdCHK2 reached higher expression level in leaves and stems (Figure 4). MdCHK4a and MdCHK4b displayed similar expression profiles with high expression level in stems whereas gene expression was hardly detected in flowers (Figure 4). MdCHK3a and MdCHK3b disclosed differential pattern of expression. While MdCHK3a was mainly expressed in roots, MdCHK3b showed its highest expression in flowers (Figure 4). Finally, substantial differences in the overall expression level of each MdCHKs were also observed. Interestingly, MdCHK2 displayed the higher expression level whilst MdCHK3b/MdCHK4a and MdCHK3a/MdCHK4b retained a 10- and a 100-fold lower expression, respectively.

To complete qPCR analysis, we used the available RNAseq data to generate an expression matrix of apple tree genes (Supplementary Table S2). The best co-expressed genes with each MdCHK through our expression matrix were investigated and compared (Supplementary Table S3). A very weak degree of overlap among lists of genes co-expressed with each MdCHK was found. MdCHK2 shared up to 18 genes with other cytokinin receptors. MdCHK3 and MdCHK4 groups did not have common co-expressed genes. MdCHK3a/MdCHK3b homologs as well as MdCHK4a/MdCHK4b homologs shared respectively 88 and 26 genes (Figure 5). This low overlapping of co-associated genes might support the limited functional redundancy of MdCHKs. Each MdCHK co-expressed genes list was compared to the list established for their Arabidopsis ortholog in order to highlight potential shared genes. We found a relatively weak overlap between functions associated to either CHK (Supplementary Figure S5A and Table S4). For example, only 19 genes were similarly co-expressed between MdCHK2 (1.8%) and AtCHK2 (8.7%). While such a weak overlap could be due to the initial datasets used to calculate correlations which differ in size and experiments, conserved co-expressed genes may be good candidates for a further investigation of the cytokinin pathway. In addition, we found very small overlaps between co-expressed gene lists of AtCHKs, as observed for MdCHKs, reinforcing a potential specificity in CHK functions (Supplementary Figure S5B). Enrichment tests of Gene Ontology (GO) terms performed on each list of co-expressed genes also gave an overview of possible specific physiological processes associated with each receptor. For example, “Embryo development ending in seed dormancy” and “Response to cadmium ion” were exclusively enriched for MdCHK2 whereas “Plant-type secondary cell wall biogenesis” and “Regulation of growth” were specifically enriched for MdCHK4a and MdCHK3b, respectively (Figure 6 and Supplementary Figure S6).

We investigated the subcellular distribution of the MdCHKs using C-terminal YFP tagging to ensure the correct anchoring of the transmembrane domains. MdCHK-YFP constructs were transiently expressed in C. roseus cells that constitute a reliable model for studying protein subcellular localization (Foureau et al., 2016). In transiently transformed cells, the fusion proteins displayed a fluorescence signal located in the endoplasmic reticulum (ER) network throughout the cell as well as in the perinuclear space (Figures 7A,E,I,M,Q). The signal perfectly co-localized (Figures 7C,G,K,O,S) with the specific ER-CFP marker (Figures 7B,F,J,N,R), confirming that the five MdCHKs are located in the ER of plant cells. Noteworthy, cytokinin addition in plant cell medium did not alter the ER localization of the MdCHK-YFP fusions (data not shown).

Since plant cells may produce their own pool of cytokinins preventing the study influence of exogenous cytokinins on the localization of the MdCHKs, we therefore investigated the subcellular distribution of MdCHKs in the yeast S. cerevisiae, by using YFP fusion proteins. In absence of cytokinins, a punctate fluorescence pattern was observed for the five MdCHKs (Figures 8G1,I1,K1,M1,O1), which accumulated in the ER forming structures comparable as organized smooth ER (Snapp et al., 2003) that is described to result from protein interactions. Upon cytokinin treatment, the ER localization of the five MdCHKs did not change. But interestingly, a reorganization of the fluorescent pattern was observed reflecting the decrease or disappearance of aggregate structures (Figures 8H1,J1,L1,N1,P1). Indeed, MdCHK2 and MdCHK4b displayed a strong perinuclear localization (Figures 8H1,P1). Concerning MdCHK3a, MdCHK3b, and MdCHK4a, the fluorescence signal appeared in a discontinuous pattern as well as in the perinuclear space (Figures 8J1,L1,N1). In order to ensure that cytokinins themselves had no impact on the architecture of the ER, we used T16H2-CFP construct as a specific ER marker (Besseau et al., 2013). Upon cytokinin treatment, no redistribution of fluorescence signal was observed, reinforcing the plausibility of a specific reorganization of MdCHKs in response to cytokinin signal (Figures 8E1,F1).

Cytokinin receptors were previously proposed to interact each other to enable the trans-phosphorylation of the HK domain after cytokinin perception (Dortay et al., 2006; Caesar et al., 2011; Hothorn et al., 2011; Wulfetange et al., 2011). Considering the multiple possibilities of interactions between the five MdCHKs, homo- and hetero-dimerization were investigated by BiFC assays in planta. The full coding sequences of MdCHKs were cloned upstream of the coding sequence of the two split-YFP fragments (YFP^N and YFP^C) to generate the MdCHK-YFP^N and MdCHK-YFP^C fusion proteins. BiFC analysis revealed that MdCHK2, MdCHK4a, and MdCHK4b were able to form homodimers within the ER network (Figures 9A,S,Y) whereas no BiFC complex reconstitution was observed when testing the MdCHK3a and MdCHK3b homodimers (Figures 9G,M). Moreover, no signal was detected with the MdCHK3a and MdCHK3b heterodimer combination (Figures 9H,L). Additionally, a cytokinin application did not result in the formation of a fluorescent signal within the three configurations (data not shown). Nevertheless, MdCHK3a and MdCHK3b were able to heterodimerize with MdCHK2, MdCHK4a and MdCHK4b within the ER (Figures 9B,C,F,I–K,N,O,Q,R,V,W). In addition, MdCHK2 and both MdCHK4 homologs shape heterodimers with each other (Figures 9D,E,P,U,T,X).