Two different approaches were used to identify MdTLPs in the apple genome. For the first approach, the published Arabidopsis TLP protein sequences were retrieved from the TIGR database¹ and used as queries in BLASTP searches against the Malus domestica genome database (Malus domestica Genome v1.0²). Previous studies failed to detect AtTLP4 in Arabidopsis using several experimental approaches, suggesting that it might be a pseudogene (Lai et al., 2004; Bao et al., 2014). Thus, we used 10 Arabidopsis TLP protein sequences as queries. To avoid excluding any possible candidates, the E-value cut-off was set to 0.001, as was done previously in a similar investigation (Li et al., 2011; Jia et al., 2013; Cui et al., 2015). Redundant sequences with the same accession numbers were removed from the data set. The retrieved sequences were then queried against the InterPro³ database to ensure the presence of tubby domains (Mitchell et al., 2015). For the second approach, the tubby-domain Hidden Markov Model Profile, which was downloaded from the Pfam⁴ database, was used as a query to search for all of the annotated proteins in the complete apple genome database² using HMMER 3.0 (Eddy, 1998). The candidate proteins’ sequences were extracted by a Perl script and then examined for the tubby domain using the InterPro database. All of the sequences identified through these two methods were submitted to SMART⁵ to ensure the integrity of the tubby domain sequence. A BLASTN search (E-value ≤E^-100) against the apple EST dataset was conducted to find the corresponding expression record for each putative family member⁶.

Multiple sequence alignments of the TLP amino acid sequences in apple (MdTLP1-9), Arabidopsis (AtTLP1-3, AtTLP5-11), rice (OsTLP1-14), and maize (ZmTLP1-15) were performed with ClustalX (Version 2.1⁷) (Saitou and Nei, 1987; Larkin et al., 2007). To make sure that all sequences were properly aligned, the sequences of the tubby domain were adjusted manually first, then let the other regions realigned, using software CLC Sequence Viewer 7⁸. An unrooted phylogenetic tree was constructed from the alignments of the full-length protein sequences according to the neighbor-joining method with 1,000 replications, and the phylogenetic tree was drawn with the MEGA5 program (Tamura et al., 2011).

Chromosomal location data were retrieved from apple genome annotations downloaded from the Genome Database for Rosaceae⁹. The chromosome map showing the physical location of all of the MdTLP genes was generated using a revised version of MapDraw (Liu and Meng, 2003). To explore the diverse exon–intron organizations of MdTLPs and AtTLPs (AtTLP1-3, AtTLP5-11), we compared the predicted coding sequences of MdTLPs and AtTLPs with their corresponding genomic sequences using GSDS software¹⁰ (Hu et al., 2015). The molecular weight and isoelectric point (pI) of the proteins were calculated with the ExPASY Compute pI/Mw Program¹¹ (Wilkins et al., 1999).

To identify potential protein motifs and detect any additional domains outside the apple tubby domains, sequences of full-length MdTLPs were queried against the InterPro database¹² (Mitchell et al., 2015). Protein subcellular locations were predicted using WoLF PSORT¹³, an extension of the PSORT II program (Horton et al., 2007).

All putative MdTLPs were analyzed by MEME (version 4.11.1¹⁴), a motif search tool for identifying conserved motifs of tubby domains (Bailey et al., 2009). To obtain the most significant conserved motifs in the nine MdTLPs, different numbers of motifs were tried with default parameters in normal mode. The identified motifs were annotated using SMART protein analysis software¹⁵.

The three-dimensional structure of the tubby domain was obtained by homology modeling using the website CPHmodels 3.2 Server¹⁶. Images were generated in the modeling package PyMOL v1.5¹⁷ (Nielsen et al., 2010).

To investigate cis-acting elements in the promoter sequences of MdTLPs, 1,500 bp of genomic DNA sequence upstream of the transcriptional start sites was obtained from the apple genome. The upstream sequences were subsequently scanned in the PlantCARE database¹⁸ for the presence of various cis-acting elements (Lescot et al., 2002).

To investigate the expression of MdTLPs under abiotic stress, 3-year-old apple (Malus sieversii) seedlings were treated with either chilling at 4°C, 20% PEG6000, 100 mM H₂O₂, or 100 μM ABA for 0, 1, 3, 6, and 9 h. The leaves and roots from five individual plants were collected, placed into liquid nitrogen and stored at -80°C until further use. With respect to the samples used for organ-specific expression, different organs were collected and also stored at -80°C. Total RNA was isolated from the leaves using the CTAB procedure (Gasic et al., 2004). The RNA concentrations and A260/A280 ratios were determined using a NanoDrop Spectrophotometer (ND-1000 Spectrophotometer, Peqlab). The integrity of the RNA samples was examined with an Agilent 2100 Bioanalyzer (RNA Nano Chip, Agilent, Santa Clara, CA, USA). Suitable RNA was used for cDNA synthesis and qRT-PCR.

cDNA fragments were synthesized from total RNA using TransScript^TM One-step gDNA Removal and cDNA Synthesis SuperMix (TransGen Biotech, Beijing, China). To ensure the cDNA samples obtained were qualified, two stress-specific genes (DREB genes: MDP0000147009 and MDP0000218344) in apple genome were selected as marker genes, which expressions have been demonstrated to be up-regulated under stress treatments (Zhao et al., 2012). Gene-specific primers for amplification from cDNA were designed based on target gene sequences using the Beacon Designers 8.10 software. The primer sequences used in this investigation are listed in Supplementary Table S1. qRT-PCR was performed with a CFX96 real-time system (Bio-Rad, USA) in a final volume of 20 μl containing 0.8 μl of cDNA, 10 μl of 2 × SYBR Premix Ex Taq (SYBR Green; Tiangen, China) and 0.8 μl of (10 μM) primers. The thermal cycling conditions were as follows: 44 cycles of 95°C denaturation for 15 s, 55°C annealing for 30 s and 72°C extension for 15 s. The apple actin gene was used as an internal control. The real-time PCR experiment was carried out at least three times under identical conditions. The relative expression levels were calculated as 2^{-(ΔCtoftreatment-ΔCtofcontrol)}. The relative expression levels of MdTLPs in stressed samples (1, 3, 6, and 9 h) were compared to the controls (0 h) with parametric one-way ANOVA at significance levels of P ≤ 0.05 and P ≤ 0.01.

To identify the TLP protein-coding genes in apple, two approaches were used. For the first strategy, the peptide sequences of the TLPs of Arabidopsis were used as BLAST queries against the apple genome (Malus domestica Genome v1.0). To ensure that potential TLPs were not excluded and to obtain credible results, the E-value was set to 0.001, as done in a similar investigation (Li et al., 2011; Jia et al., 2013; Cui et al., 2015). Using this approach, 25 potential TLPs were identified in the apple genome. To determine whether these proteins contained tubby domains, the sequences were compared against the InterPro Database. Using this approach, eight potential TLPs were identified from the apple genome. For the second strategy, the tubby domain Hidden Markov Model Profile (Pfam01167) from the Pfam database¹⁹ was used to search the apple genome. A total of 10 putative TLPs were obtained. These sequences were also analyzed using the InterPro Database, and eight potential TLPs contained a tubby domain. All of the potential TLPs obtained by the above two strategies were submitted to the SMART database²⁰ to verify the integrity of the tubby domain sequence and confirm the presence of apple TLP genes. Finally, eight genes were identified.

In our previous study, we cloned a TLP gene (named MdTLP7) from apple, whose expression increased significantly under cold stress (Du et al., 2014). However, the MdTLP7 gene sequence could not be identified in the published apple genome. To determine whether MdTLP7 was present in the apple genome and to check the expression of this and other identified apple TLP genes, the gene sequences were searched against the apple EST database at NCBI²¹. In BLASTN analysis, each of the nine gene sequences hit several apple ESTs sequences (Supplementary Table S2), which indicated that these nine genes were genuinely expressed in apple. Therefore, we believe that the apple genome contains nine MdTLP genes (Table 1). To remain consistent with our previous studies, the MdTLP7 gene name was retained, and the other eight MdTLP genes were named MdTLP1-6 and MdTLP8-9 based on their distribution in the phylogenetic tree (Figure 1). The MdTLP peptides in length ranged from 269 to 693 amino acids, with predicted molecular weights between 29.9 and 78 kDa (Table 1).

Multiple sequence alignment showed that all plant TLPs had a highly conservative tubby domain at C-terminal almost with a conservative proline residue at the beginning (Supplementary Figure S1). To gain insight into the evolutionary relationships among all plant TLP proteins, a phylogenetic tree was constructed based on the full-length amino acid sequences from apple, Arabidopsis, rice, and maize (Lai et al., 2004; Yang et al., 2008; Yulong et al., 2016). All of the TLP encoding genes from the above four plant species were divided into three distinct groups: A, B, and C (Figure 1). Of 48 members, 41 plant TLPs belonged to group A and were further divided into the four subgroups A1–A4. Group B contained five members, and group C only contained two maize TLPs. Group A and group B contained TLP members from both dicot and monocot plants. The results showed that in some cases, the evolutionary relationship of plant TLPs between dicot and monocot was closer than that among dicot or monocot plants (Figure 1). For example, MdTLP9 was closer with ZmTLP13 than with MdTLP8 in group B (Figure 1). To check whether this result was valid, a detailed sequence alignment was checked (Supplementary Figure S1). Indeed, MdTLP9 had higher similarity with ZmTLP13 than MdTLP8 in sequences, especially in the segment of core α helix in tubby domain (indicated in the box of Supplementary Figure S1). These results suggested that the main characteristics of plant TLPs in group A and group B had been established before the dicot–monocot plants split. The TLPs from apple (MdTLPs) were distributed into two groups, group A, including A1, A2, and A4, and group B (Figure 1). The members in subgroup A3, AtTLP3 and AtTLP9 from Arabidopsis, have been studied in detail, and it has been indicated that they are involved in the stress response (Lai et al., 2004; Reitz et al., 2012, 2013; Bao et al., 2014). Whether MdTLPs in the different subgroups or different groups have similar functions requires further investigation.

Based on the chromosomal distribution map of MdTLP genes generated in this study, all of the MdTLP genes were distributed across six of the seventeen apple chromosomes, including 8, 9, 10, 11, 15, and 17 (Supplementary Figure S2). Only chromosomes 8 and 11 contained two MdTLPs, while other MdTLPs were mapped to different chromosomes. MdTLP1 and MdTLP2 were both found on chromosomes 8, located near each other. They were also clustered together in the phylogenetic tree. To understand the possible structural evolution of MdTLPs, the intron-exon structure of MdTLPs and AtTLPs was analyzed in this study (Supplementary Figure S3). Intron/exon organizations for all MdTLPs were determined based on their exon position and gene length. The number of introns in the MdTLPs varied, ranging from 2 to 9 (Supplementary Figure S3). Unlike other MdTLPs, MdTLP5 had much longer introns, although the number of introns (8) was similar to that of other MdTLPs. The MdTLPs in group B shared similar intron-exon structure distribution characteristics.

Most MdTLPs in group A contained an additional functional domain, an F-box domain at the N-terminus (Figure 2). F-box-containing proteins constitute a large family in eukaryotes, which are characterized by a conserved F-box domain consisting of approximately 50 amino acids at their N-terminus. The appearance of both domains in one protein suggested that an interplay between tubby and F-box domains may play roles in physiological processes. Of nine MdTLPs, three members had no F-box domain (MdTLP1, MdTLP8, and MdTLP9) (Figure 2). Although these members were quite similar in patterns of protein structure, multiple alignments, and phylogenetic analysis showed that MdTLP1 was far from MdTLP8 and MdTLP9 (Figure 1; Supplementary Figure S1).

The subcellular location of MdTLPs was predicted, and the results are shown in Table 1. Most MdTLPs were predicted to be located in the nucleus. Only MdTLP2 and MdTLP3 in subgroup A1 were predicted to be located in chloroplasts. In mice, TLPs were mainly localized within the nucleus as transcription factors (Boggon et al., 1999). Here, we speculate that the nuclear-localized MdTLPs may also function as transcription factors. Future investigations are needed to experimentally confirm their location and transcription factor activity.

All of the MdTLP peptide sequences were submitted to MEME²², and three types of motifs were identified in the tubby domain of MdTLPs (Figure 3A). The consensus sequences of motifs are shown in Figure 3B, with the height of each stacked letter representing the probability of that amino acid appearing at each position. Motif 1 and 2 were more highly conserved and usually located near the C-terminus of the polypeptide, harboring significant tubby domain characteristics similar to other plant and animal TLPs. Motif 3 often existed in the middle of the protein with a highly conserved amino acid sequence (RGPRRM), suggesting that this sequence may have an important biological function.

Three-dimensional structures of MdTLP tubby domains were established by homology modeling of a central α helix surrounded by a β barrel (Figure 4). Some MdTLP members had a typical tubby architecture with a closed β barrel formed by 12 anti-parallel strands and a central α helix, for example, MdTLP1, MdTLP2, MdTLP7, and MdTLP8. Other members contained an incomplete β barrel (less than 12 anti-parallel strands) and a central α helix, for example, MdTLP3, MdTLP4, and MdTLP6, which had 10, 11, and 6 anti-parallel sheets, respectively. MdTLP9 consisted of a complete β barrel without a central α helix, while MdTLP5 had an incomplete β barrel without a central α helix. Differences in the three-dimensional structures may lead to the functional diversification of different MdTLPs.

Transcriptional control of stress-responsive genes is a crucial part of the plant response to a range of abiotic and biotic stresses. Transcription factors have the potential to activate or repress genes through cis-acting sequences in promoter regions that respond to specific stresses (Singh et al., 2002). In plants, some TLPs play a role in responses to abiotic stress (Wardhan et al., 2012; Bao et al., 2014; Du et al., 2014). A search for putative cis-acting elements within the 1,500 bp of the genomic sequence upstream of the MdTLP 5′-UTRs was performed. Many stress-response related cis-acting elements were found in the promoter regions of MdTLPs, including ABREs, DRE, LTRE and MYB and MYC transcription factor elements (Table 2). MYB elements exist in all of the MdTLP promoters. MYB has been demonstrated to be involved in stress-induced drought, low temperature, salt, and ABA responses (Xiong et al., 2002; Dai et al., 2007). ABRE responds to drought and ABA via ABRE binding proteins in Arabidopsis (Hobo et al., 1999; Li et al., 2012). DRE binding proteins participate in drought, salt, low temperature, and ABA responses in rice (Zhang et al., 2009). LTRE primarily contributes to the regulation of low temperature responses in poplar (Jiang et al., 1996; Maestrini et al., 2009). The presence of abiotic stress-responsive elements suggests that MdTLPs may be regulated by various stresses. In addition to the cis-acting elements mentioned above, other types of cis-acting elements were also detected and are listed in Supplementary Table S3. Two promoters of identified MdTLPs were not analyzed because their promoter sequences could not be found in the apple genome database.

Based on the promoter analysis results, MdTLPs may be associated with the abiotic stress response. Thus, to further investigate the potential functions of MdTLPs under abiotic stress conditions, cDNA samples were obtained using apple seedlings exposed to either PEG, H₂O₂, exogenous ABA, or cold stress for 0, 1, 3, 6, and 9 h. To detect the quality of cDNA samples, two stress-sensitive genes (DREB genes: MDP0000147009 and MDP0000218344) were selected as marker genes. Similar with the previous report (Zhao et al., 2012), their expressions were significantly up-regulated under different stresses both in leaves and roots (Supplementary Figure S4). As to the expression of MdTLPs, almost all of the MdTLPs were also up-regulated in both leaves and roots under different stresses and exogenous ABA. During PEG treatment, six MdTLPs (MdTLP1-5, MdTLP9) were up-regulated significantly in leaves (Figure 5A). For example, the relative transcript levels of MdTLP1 and MdTLP9 were both up-regulated by approximately 35-fold compared with the control condition. Only three MdTLPs (MdTLP6, MdTLP7, and MdTLP8) were down-regulated under PEG stress. In the roots, five MdTLPs (MdTLP3, MdTLP4, and MdTLP7-9) were up-regulated significantly in response to PEG treatment (Figure 5B). Among them, the expression of MdTLP4 changed maximum, reaching nearly 30-fold. Under cold stress, almost all MdTLPs in leaves and in roots showed significantly up-regulated transcript levels, suggesting that all MdTLPs are cold-responsive genes. Among these genes, MdTLP4 expression increased nearly 70-fold in leaves in 6 h and about 150-fold in roots in 3 h. In response to H₂O₂ or ABA treatment, five members (MdTLP1-5) were significantly up-regulated in leaves, while in roots, the expression of MdTLP4, MdTLP6, MdTLP8, and MdTLP9 showed a significant increase. Notably, MdTLP4 was up-regulated in all the treatments both in leaves and roots.

To provide further clues to the putative roles of the MdTLP genes in apple development, the expression patterns of all 9 MdTLP genes were investigated by qRT-PCR in six different organs (leaves, roots, stems, flowers, seeds, and buds), and the results are indicated in Figure 6. Relatively high levels of MdTLPs were found in roots, stems, and leaves; in flowers, the transcript levels of MdTLPs were rather low. Compared with other MdTLPs, MdTLP8, and MdTLP9 had low expression levels in most organs. MdTLP1, MdTLP2, MdTLP3, and MdTLP4 were highly expressed in roots, stems, and leaves, with lower levels of expression in flowers, seeds and buds. Conversely, the MdTLP6, MdTLP7, MdTLP8, and MdTLP9 transcript levels were relatively high in buds, which suggest that these genes play an important role in the sprouting process.