Inheritance of dwarfing traits
Phenotype of dwarfing traits, tree height, trunk diameter, and canopy width segregated significantly among ‘Red Fuji’ scions grafted on the 1123 hybrids of ‘M. baccata’ × ‘M9’ (Fig.1 A and B; Table 1). None of the frequency distributions of these dwarfing traits showed a Gaussian distribution (Fig. 2 A-C). Broad sense heritability was estimated to be 73.0%, 79.7%, and 82.3% for tree height, trunk diameter, and canopy width, respectively (Table 1).
Significant linear regressions were shown between tree height, diameter and canopy width. The Pearson’s correlation coefficients were 0.9421, 0.8663, and 0.9388 (P<0.01) for tree height/trunk diameter, tree height/canopy width, and trunk diameter/canopy width, respectively (Fig. 2 D-F).
Marker validation and candidate gene prediction
Three SSR markers, Hi01c04, Hi04a08, and CH03a09 near Dw1 on chromosome 5, and the other two markers, Hi07d11 and CH02d08 close to Dw2 on chromosome11, were genotyped in 30 hybrids where the scion exhibited the dwarfing phenotype [2]. All five of these markers were heterozygous in ‘M9’ rootstock, homologous in ‘M. baccata’ rootstock, and displayed a distorted segregation ratio in dwarfing F1 hybrids, indicating the association of these markers with dwarfing effects (Fig. 3).
To further narrow down the intervals of Dw1 and Dw2, four and nine InDel markers flanked by the above-mentioned SSRs near Dw1 and Dw2, respectively, were selected from the parental resequencing data set (Table S1). All markers were heterozygous in ‘M9’, and were homologous in ‘M. baccata’. These InDel markers were also genotyped in the 60 hybrids, the scion grafted on each of the 30 hybrids showed extremely dwarfing or extremely vigorous tree architecture. Whenever a marker recombined with the phenotypes, it’s considered not to be a diagnostic marker at all and the corresponding gene should not be the candidate gene. Of these markers, the segregation ratio of C10 and B06 near Dw1 and M90 close to Dw2 was 30:0 in dwarfing hybrids and 0:30 in vigorous hybrids (Fig. 4). According to the apple genome, C10, B06, and M90 were located in the upstream regions of MdLBD3, MdARF6, and MdG3OX3, respectively. Therefore, these three genes were selected as potential candidate genes for further validation.
Validation of allelic variation of the candidate genes
The relative expression of MdLBD3 and MdG3OX3 was significantly higher in the fine roots of ‘M9’ than in ‘M. baccata’, but MdARF6 was more extensively expressed in the roots of ‘M. baccata’ than in ‘M9’ (Figure 5A). The 0~1104 bp, 0~1289 bp, and 0~1176 bp segments upstream of the transcription start sites (TSS) of MdLBD3, MdARF6, and MdG3OX3, respectively, were cloned from the DNA of ‘M9’ and ‘M. baccata’ and Sanger sequenced (Fig. S1; Fig. 5C). Besides some SNPs that were heterozygous on ‘M9’, an 11bp and a 10 bp deletion were found -339 bp and -278 bp upstream of the TSS of MdLBD3 and MdG3OX3, respectively, while a 14 bp insertion was located -954 bp prior to the TSS of MdARF6 (Fig. 5C). Non-synonymous variations were not found on the coding sequences of these three genes (Fig. S1). To test whether the gene expression differs between genotypes of these InDels in MdLBD3, MdARF6, and MdG3OX3, fine roots of 108 F1 randomly chosen F1 hybrids of ‘M. baccata’ × ‘M9’ were sampled for a gene relative expression assay (Table 2). The relative gene expression of hybrids with MdLBD3 Del:Ins, MdG3OX3 Del:Ins, and MdARF6 Del:Del genotypes was significantly greater than that of hybrids with MdLBD3 Ins:Ins, MdG3OX3 Ins:Ins, and MdARF6 Del:Ins genotypes, respectively (Fig. 5B). These promoter variants of each gene were then constructed fusion with a GUS reporter and transiently expressed in Nicotiana benthamiana leaf epidermal cells. The promoter activity of MdLBD3-pro-Del and MdG3OX3-pro-Del from ‘M9’ was significantly stronger than that of MdLBD3-pro-Ins and MdG3OX3-pro-Ins from ‘M9’ or ‘M. baccata’, whilst the promoter activity of MdARF6-pro-Ins from ‘M9’ was lower than that of MdARF6-pro-Del from ‘M. baccata’ or ‘M9’ (Fig. 5D).
To further examine the functional variations on the promoter sequences of these candidate genes, the promoter activity was re-examined by promoter truncation. Truncating off segments of promoter sequences, including the InDels of MdLBD3, MdARF6, and MdG3OX3 in ‘M9’, destroyed the promoter activity completely. The SNPs did not interfere with promoter activity, but the InDels on the each gene’s promoter altered the transcription activity significantly (Fig.5 A, B and D). Consequently, InDels on the promoter of MdLBD3, MdARF6, and MdG3OX3 were utilized as functional markers for dwarfing effects.
Functional analysis of allelic variation of the candidate genes
Prediction revealed that the 11bp deletion ‘TCAACAACATA’ located -339 bp from the TSS of MdLBD3 may have caused the generation of a new cis-element with the critical motif ‘TGGAATGGGTATAG’, which could be bound by MdWRKY2 transcription factor (Fig. 6A). Yeast one hybrid assay (Y1H) confirmed the binding of MdWRKY2 protein to the promoter of MdLBD3-Del in ‘M9’ but not to MdLBD3-Ins in ‘M. baccata’ (Fig. 6B). No differences in MdWRKY2 relative expression were identified between either ‘M9’ and ‘M. baccata’ nor their F1 hybrids with dwarfing and non-dwarfing effects (Fig. 6A).
The 14 bp insertion ‘GAGAGAGCGAGAGA’, on the promoter of MdARF6 in ‘M9’, destroyed the MdGAMYB binding motif ‘GGCGGTTGTG’ based on the prediction from PlantCare software (Fig. 6C). Y1H showed that the binding ability was completely absent between AD-MdGAMYB and BD-MdARF6-Ins (Fig. 6C). The MdGAMYB expression levels were quite analogous among ‘M9’, ‘M. baccata’, and their F1 hybrids, irrespective of the dwarfing effects (Fig. 6A).
The 10 bp deletion ‘ATCATATTTT’, on the MdG3OX3 promoter, produced a new MdABI5 binding motif ‘ACGTTTTTGT’, which was confirmed by Y1H assay (Fig. 6D). Similar MdABI5 expression levels were observed among ‘M9’, ‘M. baccata’, and their F1 hybrids, irrespective of the dwarfing effects (Fig. 6A). These data suggested the use of the InDel variations MdLBD3, MdARF6, and MdG3OX3 as potential markers for dwarfing effects.
Prediction of dwarfing ability with marker genotype estimated genetic values (GEGV)
To estimate the genotype effects of the functional markers, MdLBD3(Del/Ins) /MdARF6(Del/Ins) /MdG3OX3 (Del/Ins) (abbreviated as Ld/Li, Ad/Ai, and Gd/Gi), 108 F1 hybrids of ‘M. baccata’ × ‘M9’ were randomly chosen and genotyped. The marker genotype effects and effects of marker genotypic combinations on tree height, trunk diameter, and canopy width of the grafted trees were estimated (Table 2). AdAi had the largest marker genotype effect, and LdLi had the least (Table 2). The marker GEGV of the F1 hybrids were calculated by adding up the genotype effects of the three functional markers and the overall mean phenotype value. The Pearson’s correlation coefficients for tree height, trunk diameter, and canopy width between GEGV and observed phenotype value (OPV) were all significant (Fig. 7A, B and C). Five-fold cross validation produced an average correlation coefficient of r = 0.93 (P<0.01), indicating a very high predictability of GEGV in the current hybrid population. To test the applicability of this GEBV in Malus germplasm, 64 Malus rootstock accessions were genotyped for Ld/Li, Ad/Ai, and Gd/Gi markers. The GEGVs were consistent with the accession dwarfing ability reported by previous literatures, implying that the GEGV model should be applicable for a broad spectrum of apple rootstocks (Table 3).