Comparative genomic analysis of the IDD genes in five Rosaceae species and expression analysis in Chinese white pear (Pyrus bretschneideri)

Identification of IDD genes in five rosaceae plants

In this study, the Chinese white pear genome database was obtained from the Pear Genome Project (http://gigadb.org/dataset/10008). The genome databases of three Rosaceae plants (F. vesca, R. occidentalis, Prunus avium) were downloaded from genome database for Rosaceae (GDR) (https://www.rosaceae.org/). The genome sequence data of Prunus mume were downloaded from the Prunus mume genome project (http://gigadb.org/dataset/10008). DNATOOLS software was used to create a local database containing the amino acid sequences of the IDD genes in the five Rosaceae plants (Cao et al., 2016). To identify the IDD genes of these five Rosaceae species, we used the following methods. First, the A. thaliana (16) and Malus domestica (20) IDD gene sequences were collected for use as query sequences in a TBlastN search with the default E-value (Table S1). We compared these sequences with the local database sequences of the above-mentioned five Rosaceae species. Using the SMART online software program, the IDD candidate gene sequences initially screened by BLAST were then tested to determine whether they contained a zinc finger domain (http://smart.embl-heidelberg.de/) (Letunic, Doerks & Bork, 2012). Each candidate sequence containing a zinc finger domain was subsequently manually checked for an IDD domain, and the protein sequences lacking a full IDD domain and redundant sequences were discarded. We then isolated the candidate IDD gene sequences and used the ExPASy online site to investigate the molecular weights of the IDD proteins (http://web.expasy.org/protparam/) (Artimo et al., 2012). The subcellular localization of all IDD proteins was predicted using the WoLF PSORT website (http://www.genscript.com/wolf-psort.html) (Horton et al., 2007). Phtre2 website (http://www.sbg.bio.ic.ac.uk/phyre2) for protein three-dimensional structures prediction (Kelley et al., 2015). Gene ontology (GO) annotations analysis we use BLAST2GO software (Conesa et al., 2005).

Phylogenetic analysis

Sequence alignment of all IDD proteins was done using the ClustalW tool in MEGA 7.0 software. The phylogenetic tree was constructed with MEGA 7.0 software using the maximum likelihood (ML) (bootstrap = 1,000) (Kumar, Stecher & Tamura, 2016). The Malus domestica IDD gene sequence used in phylogenetic tree was derived from the results of Fan et al. (2017). The sequences of A. thaliana, O. sativa and Z. mays were obtained from the article by Colasanti et al. (2006).

IDD gene structures and conserved motif prediction

INDETERMINATE DOMAIN gene structures were compared using Gene Structure Server (http://gsds.cbi.pku.edu.cn) (Guo et al., 2007). The motifs of the IDD genes in these five Rosaceae species were analyzed using MEME online analysis software (http://meme.sdsc.edu/meme4_3_0/intro.html) (Bailey et al., 2015). Parameters for the conserved motif prediction were motif width greater than 6 and less than 200. The number of identified motifs was 20.

Chromosomal location, gene duplication and Ka/Ks ratio analysis

Five Rosaceae species chromosome start positions and other relevant information about the IDD genes were obtained from the public genome database. The chromosomal physical locations of the IDD genes of all five Rosaceae species were mapped using MapInspect software (http://mapinspect.software.informer.com) (Hu & Liu, 2011; Lin et al., 2011; Zhu et al., 2015). To define gene duplication events, we mainly relied on the following principles: (1) The similarity of the two gene-coding sequences was more than 80%. (2) If the two genes were located on the same chromosome and separated by at least 200 kb, we considered these two genes tandem-duplicated genes. (3) If the two genes were located on different chromosomes, they were called segmentally duplicated genes (Long & Thornton, 2001). DnaSP v5.0 software was used to calculate non-synonymous (Ka) and synonymous substitution (Ks) values and perform sliding window analysis. Parameters for sliding window analysis were window size 150 bp and step size nine bp (Librado & Rozas, 2009).

Microsynteny analysis and Chinese white pear IDD gene promoter Cis-acting element analysis

We obtained the promoter sequence of each IDD gene from the Chinese white pear genome database, including the DNA sequence of the initiation codon (ATG) located 1,500 bp upstream of each IDD gene. The online software Plantcare was used to analyze the cis-acting elements of the promoter region (http://bioinformatics.psb.ugent.be/webtools/plantcarere/html/) (Rombauts et al., 1999). The microsynteny of six Rosaceae species (we added the IDD genes of Malus domestica) was identified using the Multiple Collinearity Scan toolkit (MCSscanX) (Abdullah et al., 2018).

RNA extraction and qRT-PCR

The plant material used in this experiment was the “Dangshan Su” pear, which grows in the Yuanyichang agricultural park (Dangshan County, Anhui Province, China). Samples of current-year flowers, flower buds, stems, mature leaves and fruits were obtained from 40-year-old Pyrus bretschneideri cv. Dangshan Su trees. The fruits were picked 15 days after pollination (DAP), 39 DAP, 47 DAP, 55 DAP, 63 DAP, 79 DAP and 145 DAP, and the fruits obtained 39 DAP were used for the expression analysis in different organs. “Dangshan Su” pear buds were treated with gibberellin (GA) and sucrose. Specifically, on a clear morning, terminal buds on current-year spurs were sprayed (<5 cm) with GA (700 mg·L⁻¹) and sucrose (20,000 mg·L⁻¹) using a low-pressure hand wand sprayer. Some plant materials were collected immediately after spraying with GA and sucrose and stored at −80 °C for use as controls. Materials collected 2 h (H), 4 H, 6 H, 8 H, 12 H and 24 H After GA and sucrose treatment were stored at −80 °C for gene expression analyses and samples collected 0 h post-treatment (HPT), 2 HPT, 4 HPT, 6 HPT, 8 HPT, 12 HPT and 24 HPT were also collected. RNA was extracted from the plant materials (including fruits and various organs) using a plant RNA extraction kit from Tiangen (Beijing, China). Reverse transcription was performed using a PrimeScript^™ RT reagent kit with gDNA Eraser (Takara, Kusatsu, Japan), and each reaction consisted of one μg of RNA. The qRT-PCR primers for the pear IDD genes were designed using Beacon Designer 7 software (Table S8). Each 20-μL qRT-PCR system consisted of 10 μL of SYBR® Premix Ex Taq^™ II (2×) (Takara, Kusatsu, Japan), two μL of cDNA, 6.4 μL of water and 0.8 μL of PbIDD-F and PbIDD-R. The pear Tubulin gene (No. AB239680.1) (Wu et al., 2012) was used as an internal reference. The reaction procedure was performed following the instruction manual, and three repetitions were run for each sample. The relative expression levels were calculated using the 2^−ΔΔCt method (Livak & Schmittgen, 2001).

Identification, characterization and phylogenetic analysis of IDD genes

A total of 68 IDD proteins were identified and used for further analysis (Table 1; Table S1), which included 16 IDD proteins (PbIDD1-PbIDD16) in Chinese white pear. We also identified a total of 52 IDD proteins in the other four species, including 14 in F. vesca (FvIDD1-FvIDD14), 13 in Prunus mume (PmIDD1-PmIDD13), 14 in R. occidentalis (RoIDD1-RoIDD14) and 11 in Prunus avium (PaIDD1-PaIDD11). Although only 11 IDD proteins were identified in Prunus avium, the numbers of IDD proteins in the other species were very similar but fewer than the number of IDD proteins in Malus domestica (20). We used ExPASy software to calculate the physicochemical parameters of the IDD genes (Table S2). We identified 68 IDD genes with isoelectric points (pIs) greater than 7. In these Rosaceae species, the lowest pI value was 7.76 (PaIDD8), whereas the highest pI value was 9.45 (FvIDD4). The molecular weights of these IDD genes were quite different, ranging from 42.0 kDa for PaIDD10 to 79.6 kDa for PmIDD8. The numbers of amino acids in the IDD proteins were also quite different. PbIDD14 contained the fewest (only 381) amino acids, whereas, PmIDD8 had the largest number (728). The grand average hydropathicity (GRAVY) of all 68 IDD genes was less than 1, and the smallest GRAVY value (−0.914) was obtained for PbIDD15.

A phylogenetic tree of 98 IDD proteins was constructed using the ML method (Fig. 1). Based on the phylogenetic tree, we identified four phylogenetic groups (groups 1–4). Group 1 had the most members and contained all six Rosaceae species and 44 IDD genes: Pyrus bretschneideri (8), F. vesca (7), Prunus mume (6), R. occidentalis (7), Prunus avium (6), Malus domestica (10). Groups 3, 2 included 21 and 14 IDD genes, respectively. Group 4 had the fewest members: two IDD genes in Pyrus bretschneideri, two IDD genes in F. vesca, two IDD genes in Prunus mume, two IDD genes in R. occidentalis, one IDD gene found in Prunus avium and no IDD gene in Malus domestica. ZmID1 was grouped into a class by itself, and most PbIDD genes were more tightly grouped with MdIDDs. In addition, the IDD genes of some Rosaceae species were closely related to those of A. thaliana and O. sativa: PbIDD3, PbIDD5, FvIDD3, RoIDD4, PmIDD3 and PaIDD4 were included in group 1 and were closely related to OsIDD2; FvIDD12, RoIDD9, PmIDD11, PaIDD5 and AtIDD3 and 8 are closely related, and AtIDD1 was included in group 3 and clustered with seven Rosaceae IDD genes. The results of GO annotations (Table S9) showed that all 68 IDD genes had the function of nucleic acid binding (GO:0003676).

Structural and conserved motif analysis of IDD proteins

To more comprehensively analyse the structural diversity of the IDD genes in the five Rosaceae species, we created exon-intron pattern maps of all 68 IDD genes (Fig. 2). A total of 25 of the 34 IDD genes in group 1 had three exons and the other six genes had two exons. Additionally, FvIDD7 (4), RoIDD6 (5) and PmIDD8 (6) had a large number of exons. The analysis of group 2 revealed that FvIDD11 had two exons, PbIDD7 had four exons and all the other members had three exons. Group 3 comprised a subfamily with a complex gene structure: 11 genes had four exons and four genes had three exons. All members of group 4 had a gene structure comprising three exons. The lack of high similarity among the 68 IDD genes allowed a better understanding of the conserved motifs in these IDD genes. We used MEME software to identify 20 conserved motifs (Fig. 2; Table S3). Motifs 2, 3, 1 and 4 encoded the ID domain, whereas the remaining motifs did not have functional annotations. As shown in Fig. 3, motifs 2, 3, 1, 4 represented the zinc fingers ZF1, ZF2, ZF3 and ZF4, respectively. The four conserved motifs constituted a conserved ID domain, whereas these four conserved motifs represented two C2H2 and two C2HC structures (ZF1 and ZF2 belonged to C2H2, and ZF3 and ZF4 belonged to C2HC). Two C2H2 and two C2HC structures were identified through a sequence alignment of the conserved ID domains of five species and the results showed that the ID domain was highly conserved. The cysteine (C) and histidine (H) residues in each zinc finger domain were conserved in each species, consistent with the results of previous studies (Figs. S1–S5). All 68 IDD genes had motifs 1–4, which were the most conserved. In addition, some members of the same group and the more closely related members had highly similar motif compositions (e.g., FvIDD7 and RoIDD6 in group 1; PmIDD10 and PaIDD9 in group 2). These genes were closely related and had exactly the same motif compositions. We also identified certain group-specific motifs, such as motifs 12, 17, 20 in members of group 4 and motif 19 in some members of group 1.

We also identified a NLS in the N-terminal border of these IDD genes (Fig. S6). According to previous studies, this NLS is usually KKKR or KRKR (Colasanti et al., 2006). Not all members in the five Rosaceae species, PbIDD13, 14, FvIDD10, 11, 12 and PmIDD11, have a NLS. In addition to the conserved IDD domain, we also found two C-terminal domains (MSATALLQKAA and TRDFLG), which are encoded by motifs 6, 7, respectively, in some members (Fig. 2; Figs. S1–S5). All nine members of group 4 lacked both motifs 6 and 7. Among the members of group 1, PaIDD10 did not contain motif 6, and PbIDD14 and PmIDD8 lacked motif 7. All members of groups 2, 3 contained these two motifs. The lack of these two conserved sequences at the C terminus of some IDD genes might explain why the members of group 4 were clustered into a single category and might also underlie the differences among the functions of these proteins (Fig. 1).

Chromosomal location and duplication events of IDD family genes in rosaceae

Based on the complete genome sequences of the five Rosaceae species, the exact chromosomal locations of all 68 IDD genes were identified (Fig. 4). In Pyrus bretschneideri, the 16 IDD genes were distributed on chromosomes 3, 4, 8, 9, 11, 12, 14, 15, 16, 17. In F. vesca, the IDD genes were distributed on chromosomes 1, 2, 3, 4, 5, 6. In R. occidentalis, the IDD genes were mainly distributed on chromosomes 2, 4 and 6. Several other IDD genes were distributed on chromosomes 1, 3 and 5. In Prunus mume, except PaIDD5, which could not be mapped to any chromosome, the IDD genes were distributed on chromosomes 1, 2, 4, 6, 7, 8. However, in Prunus avium, four IDD genes were distributed on chromosome 1 and the remaining 7 IDD genes were distributed on chromosomes 3, 5, 6, 7, 8. Only four pairs of duplicated genes were identified in pear (Fig. 4). All duplicated genes identified in Pyrus bretschneideri were segmental duplications.

To clarify the driving forces of gene duplication and explore the impact of these genes on evolutionary processes, we calculated the Ka, Ks and Ka/Ks ratio for the four duplicated gene pairs in Pyrus bretschneideri. Ka/Ks = 1 is the cut-off value that indicates neutral selection, Ka/Ks < 1 represents negative selection and Ka/Ks > 1 represents positive selection (Bitocchi et al., 2017). For all four duplicated gene pairs identified in Pyrus bretschneideri, the Ka/Ks values were less than 0.2726 (Table S4). In the evolutionary processes of genes, positive selection may be overshadowed by strong negative selection (Han et al., 2016). Therefore, to comprehensively explain the selection pressure of IDD genes, we performed a sliding window analysis of the four pairs of IDD paralogues in pears (Fig. S7). Most coding site Ka/Ks ratios were less than one, with exceptions for one or several distinct peaks (Ka/Ks > 1). The Ka/Ks ratios of the conserved IDD domains were less than 1.

Microsynteny and evolutionary analysis of the IDD gene family

For the identification of interspecific microsynteny, we visualized the location of homologous and orthologous genes. A total of 30 orthologous gene pairs were identified in five Rosaceae species (Fig. 5; Table S5) and these included one, eight, three, six and 12 collinear gene pairs between Pyrus bretschneideri and F. vesca, Prunus mume, R. occidentalis, Prunus avium, Malus domestica, respectively. Notably, PbIDD2 also showed a collinear relationship with the members of the MdIDD, PmIDD and PaIDD families. In contrast, two IDD genes from pear matched one pair of Prunus mume or R. occidentalis genes. For example, PbIDD9 and PbIDD10 were orthologous to PmIDD9, whereas PbIDD9 and PbIDD10 were orthologous to RoIDD7.

However, we did not identify any orthologous gene pairs between Pyrus bretschneideri and O. sativa, Z. mays and A. thaliana. We also compared the other four Rosaceae species with these three species, and the results did not reveal any orthologous gene pairs. To further study the evolution of IDD genes, we collected the IDD genes from four monocotyledons (M. acuminata, Sorghum bicolor, Z. mays and O. sativa) and two dicotyledons (Vitis vinifera and A. thaliana). A phylogenetic tree of the genomes of these 12 species was obtained to perform an evolutionary analysis of the IDD gene family (Fig. S8A). Several Rosaceae species exhibit a closer relationship and are more distant from monocots. This result is basically consistent with the results of the phylogenetic tree constructed from all IDD genes of 12 species (Fig. S8B). The results identified 12 classes (classes I–XII). The Pyrus bretschneideri and Malus domestica, Prunus avium and Prunus mume IDD genes often cluster into one class, and the IDD genes of monocotyledons are closely related (i.e., classes X and XI are only composed of IDD genes of monocotyledons). ZmID1, OsID1 and SbID1 are closely related, but we did not identify an orthologous gene for the ID1 gene in several other species. Class III contained all species except Malus domestica, and the IDD genes in this subfamily did not have complete MSATALLQKAA and TRDFLG domains.

Analysis of Cis-acting elements in the promoters of the IDD genes in Chinese white pear

Plant IDD genes may be involved in multiple biological processes, and hormones affect their expression. We predicted the cis-acting elements of 16 Chinese white pear IDD gene promoters (Fig. 6; Table S7). A total of 13 IDD genes (except PbIDD7, 9, 10) contained at least one G-Box, which is a light-responsive element, indicating that the expression of these genes may be regulated by light. Eight Pyrus bretschneideri IDD genes had HSE, and 11 members contained MBS components. TC-rich repeats, which are cis-acting elements involved in defence and stress responsiveness, were identified in PbIDD1, 6, 7, 9, 10, 12, 14 and 16. In addition, there were many cis-acting elements related to the responses to hormones, including responses to Abscisic acid response element (ABRE), methyl jasmonate (MeJA) (CGTCA-motif, TGACG-motif), GA (P-box, GARE-motif), auxin (TGA-element), and salicylic acid (TCA-element). The auxin-related TGA-element was identified in only five members (PbIDD1, 2, 4, 8, 13), which was the least number of elements components found. A total of 41 Cis-acting elements associated with MeJA were found in nine members, and these were the most abundant Cis-acting elements. Cis-acting elements responding to GA were the most widespread, and twelve of the Pyrus bretschneideri 16 IDD genes had this element. The related Cis-acting elements of abscisic acid and salicylic acid were identified in 10 and nine members, respectively. We found that PbIDD1 has Cis-acting elements related to the above five hormones, whereas PbIDD9 only has cis-acting elements that are responsive to MeJA. Of the 16 pear IDD genes, 12 contained 19 CGTCA-motifs and TCA-elements appeared in the promoter regions of 10 members, for a total of 22 genes. In addition, the CAT-box and the CCGTCC-box were found in 11 members (PbIDD1, 2, 3, 4, 5, 6, 7, 8, 9, 12, 14); these are Cis-acting elements related to meristem expression.

Expression characteristics of Chinese white pear IDD genes

To obtain a more in-depth understanding of the Chinese white pear IDD gene family, we investigated the expression patterns of PbIDDs in different organs (Fig. 7). No PbIDD14 expression was detected in any of the organs of Chinese white pear. Among the remaining 15 members, PbIDD4, 9, 11 were highly expressed in buds but showed extremely low expression in other organs. PbIDD1 and 10 were mainly expressed in flowers. PbIDD2, 6, 8 and 15 were mainly expressed in fruits, flowers and buds. Among these four genes, PbIDD2 and 8 showed the highest expression levels in buds, whereas PbIDD6 and 15 exhibited the highest expression levels in fruit. PbIDD5, 12, and 16 were detected in multiple organs, but the main difference in the expression of these genes was found in the leaves: (1) PbIDD12 was highly expressed in leaves; (2) although PbIDD5 was expressed in all five organs, its lowest expression was detected in the leaves; and (3) PbIDD16 showed almost no expression in leaves. PbIDD3, 7 and 13 were mainly expressed in the fruit of Chinese white pear, and PbIDD3 and 7 showed almost no expression in other organs.

We investigated the expression patterns of 16 IDD genes in Chinese white pear fruit at seven developmental stages (Fig. 7). The expression level of PbIDD1 peaked at 39 DAP and was low during the other six fruit developmental stages. The expression peaks of PbIDD2, 7, 8, 11, 13, 16 also appeared at 39 DAP, but high expression of these genes was also detected at the early developmental stages (15 and 47 DAP). At the late developmental stages of fruit (79 and 145 DAP), these genes were expressed at a moderate level. PbIDD6, 15 exhibited a similar expression pattern as the above-mentioned genes, with the exception that PbIDD6, 15 were detected at the early developmental stages of fruit but were expressed at extremely low levels at the late developmental stages of fruit. The expression level of PbIDD4 was higher only at the early stage of fruit development (15–55 DAP) and at 79 DAP. PbIDD9 had a high transcription level at 63 and 145 DAP, and PbIDD12 was mainly expressed at 79 DAP. PbIDD10 had two expression peaks, one at 39 DAP and another at 55 DAP. PbIDD3 and 5 had the same expression pattern: these two genes were mainly expressed at the early stage of fruit development, and their expression level then gradually decreased.

Gibberellin and sugar response pattern analysis of PbIDDs

Chinese white pear PbIDD2, 5, 6, 8, 9, 12, 16 are mainly expressed in bud and flower. We treated the bud of Chinese white pear with exogenous GA and sucrose to investigate the expression pattern of PbIDD2, 5, 6, 8, 9, 12, 16 under GA and sucrose treatment (Fig. S9). Under exogenous GA treatment, PbIDD2, 5, 6, 8, 9, 12, 16 showed three different expression patterns (activation, inhibition and invariance). There was no significant change in the transcriptional level of PbIDD5, 9 under exogenous GA treatment. The expression levels of PbIDD2, 8 continued to decrease after exogenous GA treatment, indicating that the expression of PbIDD2, 8 was inhibited by GA. PbIDD6, 12 showed a significant increase in expression after GA treatment (reaching a peak at 3HPT); thereafter, the expression level gradually decreased to untreated levels. PbIDD16 reached the maximal expression levels at 2, 4HPT. Overall, the expression of PbIDD6, 12, 16 showed activation by GA. The transcription levels of PbIDD2, 8 increased gradually after sucrose treatment in Chinese white pear bud. The expressions of PbIDD5, 12, 16 were not always induced by sucrose. The lowest expression levels of PbIDD5, 12 were found at 4HTP, while the lowest expression level of PbIDD16 was observed at 1HTP. The expressions of PbIDD6, 9 were inhibited by sucrose.