Phenotype segregation and analysis of inheritance
By using leafy cutting propagated plants of F1 hybrids from Malus robusta Rehd. ‘Baleng Crab (BC)’ × M. pumila Mill. ‘M9’, salt injury index (SID), alkali injury index (AID), and salt-alkali injury index (SAID) were phenotyped in 2015–2017. The phenotype exhibited broad segregation ranging with averages of 0.46, 0.48, and 0.57 for SID, AID, and SAID, respectively (Fig. 1A; Fig. S1; Table 1; Table S1-S3). The overall mean phenotype value of the population for SID, AID, and SAID was approximately equal to that of to their mid-parent value, respectively (Table 1). The segregation spectrums were nearly Gaussian distributed (Fig. S1). The broad sense heritability of SID, AID, and SAID was 63.28%, 62.26% and 60.69%, respectively (Table 1).
The Pearson’s correlation coefficient was 0.1474 (p < 0.05) between SID and AID, indicating a weak positive correlation. The Pearson’s correlation coefficient between SID and SAID was -0.0172 (p > 0.1). For AID and SAID, the Pearson’s correlation coefficient was -0.1959 (p < 0.01) (Fig. 1C), indicating a negative correlation between AID and SAID. We found that not all of the salt-tolerant plants were alkali-tolerant, while salt-sensitive plants were not necessarily all alkali-sensitive (Fig. 1B).
Identification of QTLs via BSA-seq
Based on the phenotyping data from 2015–2017, six segregant bulks were developed, each including 11–26 hybrids with extreme phenotypes for tolerance or sensitivity to salt, alkali, or salt–alkali stress (Fig. 2A; Table S4). Sequencing of the six bulk DNA samples generated a total of 611,037,786 clean reads, of these reads, 97.77% were mapped to the GDDH reference genome and 74.54% were uniquely mapped (Table S5). Seventeen QTLs for SID were mapped on chromosomes 2, 7, 11, 14, 15, and 16 using BSATOS software [43]. This included three QTLs with large G’ values, S-M02.1, S-BC15.1, and S-H16.1 (Fig. S2). For AID, 13 QTLs were identified, and the G’ value of A-BC09.1, A-M06.1, and A-M16.1 were greater than the other 10 QTLs (Fig. S3). Only two QTLs were detected to be associated with SAID, and the G’ value of SA-BC16.2 was relatively large (Fig. S4). Among these QTLs, S-H16.1, A-M16.1, and SA-BC16.2 on chromosome 16 overlapped each other (Fig. 2B); QTLs A-M02.1, S-M02.3, and S-H02.1 on chromosome 2 were partially overlapping (Table S6)
Candidate gene prediction via parental resequencing data
DNA sequences of 2,657 genes within the 32 QTL intervals were downloaded from the GDDH apple genome database, and then the gene expression was analyzed using RNA-seq data. Of the total number of genes, 199 were excluded because they did not show expression. Another 410 genes were culled because the genetic variations were inconsistent with the parent from which the QTL was mapped. We also excluded 883 genes for which functional annotation was not related to the targeted trait. In addition, 475 genes with variations on the promoter region but absent from the list of differentially expressed genes (DEGs) were excluded. Finally a total of 425, 248, and 17 candidate genes were predicted for SID, AID, and SAID, respectively (Table S7-8).
RNA-seq analysis for salt, alkali, and salt–alkali tolerances
To show the differences in gene expression between hybrids that were tolerant or sensitive to salt, alkali, and salt-alkali stress, RNA-seq was performed. RNA-seq generated 289.94 Gb clean data with a Q30 value of 92.85. For each sample, over 3.04 Gb clean reads were generated, with an average GC content of 46.93%. The average proportion of total reads mapped to the reference genome was 86.05% (Table S9-S11). A total of 875, 658, and 5,359 DEGs were identified between hybrids which were tolerant/sensitive to salt, alkali, and salt-alkali, respectively (Table S12- S14).
The expression profiles differed among salt, alkali, and salt–alkali stress tolerant hybrids compared to sensitive hybrids. The expression of most abscisic acid (ABA) signaling related genes were lower throughout the experiment in saline or alkaline sensitive hybrids; some ABA related genes were considerably higher in salt–alkali sensitive hybrids compared to the other hybrids (Fig. 3A-3C). Conversely, the transcription of genes related to secondary metabolites were upregulated in saline and alkaline sensitive hybrids (Fig. 3A and C). In salt or salt–alkali tolerant hybrids, the expression of nearly all the abiotic stress response genes were higher, but most ethylene biosynthesis or signaling genes were expressed at lower levels (Fig. 3A and C). A large proportion of cell wall associated genes had a greater level of expression in alkali and salt–alkali sensitive hybrids (Fig. 3B and C). Most jasmonate (JA) related genes were highly expressed in alkali sensitive hybrids, while exhibited upregulation in salt–alkali sensitive hybrids (Fig. 3B and C).
To test the interaction among DEGs, the co-expression network was analyzed. The majority of DEGs responded at early stage (day 0–4) to salt stress than to alkali (at day 4) or salt–alkali stress (at day 3) (Fig. 4). A GATA transcription factor (TF) gene, MD16G1234900, and several ubiquitin related genes co-expressed with DEGs in tolerant to salt, alkali, and salt–alkali hybrids compared to sensitive hybrids (Fig. 4; Table S15). MD16G1234900 was located in the confident interval of QTL SA-BC16.1 (Table S8). Similarly, 13 transcription factors co-expressed with a group of DEGs between hybrids that were tolerant to saline or alkaline stress (Fig. 4; Table S15). In response to salt stress, a specific co-expression network was found to include a GATA12 gene (MD15G1092500) and several heat shock protein genes. Another co-expression network of 29 DEGs responded to alkaline stress (Fig. 4; Table S15). MD16G1235900, a voltage-gated potassium channel subunit beta gene MdKCAB, also located within QTL SA-BC16.1 interval, co-expressed with DEGs of saline–alkaline tolerant hybrids (Table S8).
Development of GAP models
To develop the GAP models, a SNP was chosen from every candidate gene to design markers for kompetitive allele-specific PCR (KASP). A total of 17, 13, and two KASP markers were designed for the traits SID, AID, and SAID, respectively (Table S16).
The marker SH14275 from QTL S-H14.1 exhibited the largest marker effect (0.4235); the marker effect (0.0354) of SM1111 from QTL S-M11.1 had the least effect on SID (Table S17). For alkali tolerance, the largest marker effect (0.2134) on AID was estimated in AB1614 derived from QTL A-BC16.1, whereas marker nn174 from QTL A-M06.1 (0.0119) had the lowest marker effect (Table S17). For SAID, the effects of both the markers SAB161 and SAB162 designed on QTLs SA-BC16.1 and SA-BC16.2, respectively, were as large as 0.2391 and 0.2947 (Table S17).
To calculate genomics predicted value (GPV), the genotype effects of all markers for SID, AID, and SAID were added to the overall value of the mean phenotype for an individual hybrid (Table S18–S20). The prediction accuracy, represented by Pearson’s correlation coefficients between GPV and the observed phenotype value (OPV) of individual hybrids, was 0.6569, 0.6695, and 0.5834 for SID, AID, and SAID, respectively (Fig. 5). Five-fold cross validation confirmed the accuracy of the predictability of GPV by 0.6587, 0.6573, and 0.6534 for SID, AID, and SAID, respectively (Table S21).
In the simulative selection, when GPV criteria were set as 0.2, the selection rates were 16.84% and 18.06% for SID and AID, respectively (Table 2). The selection efficiencies, represented by the percentage of individuals for which OPV is consistent with the corresponding GPV, were 58.33% and 40.74% for SID and AID, respectively (Table 2). For SAID, if the GPV criterion was 0.32, the selection rate was as high as 42.78%, and the selection efficiency was as relatively low as 31.33% (Table 2).
Validation of functional variations in the candidate gene MdRGLG3
From the overlapping QTLs, S-H16.1, A-M16.1, and SA-B16.2 on chromosome 16, MdRGLG3 (MD16G1282700) was predicted as a candidate gene for SAID. The full length coding sequence (CDS) of MdRGLG3 was 1,107 bp, encoding an E3 ubiquitin-protein ligase with 369 amino acid residues. In the CDS of MdRGLG3, two nonsynonymous heterozygous SNPs were found at +182 bp (SNP182) and +932 bp (SNP932) from the ATG codon in the maternal parent ‘BC’. SNP182 G/T may change a leucine to an arginine at the vWFA-domain; SNP932 T/C leads to an amino acid substitution of valine with alanine (Fig. 6A). Two nonsynonymous SNPs, SNP168 C/G at +168 bp and SNP936 G/C at +936 bp from the ATG codon, were detected in the pollen parent ‘M9’, causing a change from phenylalanine to leucine and from lysine to asparagine, respectively (Fig. 6A). These SNPs were then confirmed by Sanger sequencing (Fig. 6A; Fig. S5). There were 16 SNPs and insertion/deletions (InDels) located in the 1,600 bp DNA fragment upstream of the ATG codon of MdRGLG3 in ‘M9’ but not ‘BC’ (Supplementary File 1). In the root tissue of ‘BC’ and ‘M9,’ no differences in MdRGLG3 expression were detected in fragments per kilobase per million (FPKM) by RNA-seq or in the relative expression determined by qPCR, indicating that the promoter variations did not alter the gene expression (Supplementary File 1; Fig. S6A).
To further confirm the function of the variations in MdRGLG3 CDS, the variants were transformed into apple callus. The in vitro growth of both MdRGLG3 over-expression and MdRGLG3-RNAi apple calli was similar with the untransformed wild type (WT) under normal conditions. However, the MdRGLG3-RNAi line exhibited severely less fresh weight and higher malondialdehyde (MDA) content 14 days after treatments with either salt, alkali or salt–alkali (Fig. 6B and D). After 14 days under salt, alkali, or salt–alkali stress conditions, the fresh weight of 35S:MdRGLG3-182G apple calli was significantly higher, whereas the MDA content was significantly lower, than the other transformants and the WT (Fig. 6C and D). The phenotype values of SID, AID, and SAID in F1 hybrids with the MdRGLG3 SNP182G allele were significantly lower than that in hybrids without the SNP182G allele (Fig. 6E). SNP182G was closely linked to the tolerance genotype of the KASP marker A-BC16.1 (Table S22). These data indicated that MdRGLG3 SNP182G allele was a functional variation leading to increased tolerance to salt, alkali, and salt–alkali conditions. One hundred Malus accessions were then genotyped for the MdRGLG3 SNP182 allele; the genotype frequency of the SNP182G:G homozygote was as low as 1%, compared with 27% and 72% for SNP182G:T and SNP182T:T, respectively, indicating a possibly lethal effect of the SNP182G:G genotype (Table S24).
Validation of functional variations in the candidate gene MdKCAB
From QTL SA-B16.1, MdKCAB (MD16G1235900) on chromosome16 was predicted as a candidate gene for salt–alkali tolerance. The full length CDS of MdKCAB was 1,002 bp, encoding 334 amino acid peptide. The expression of MdKCAB did not differ between the parents, ‘BC’ and ‘M9’ (Fig. S6B). No variations altering the cis-element in the 1.5 kb upstream region were detected in either ‘BC’ or ‘M9’ (Supplementary File 2). However, two nonsynonymous SNPs were detected at +11 bp and +761 bp from the ATG codon in ‘BC’ but not in ‘M9’ (Fig. 7A). These SNPs were confirmed by Sanger sequencing (Fig. S7). SNP11 A/G and SNP761 G/A led to a substitution from aspartic acid to serine and from arginine to lysine, respectively (Fig. 7A; Fig. S7).
Under normal conditions, transformed apple calli overexpressing MdKCAB or MdKCAB-RNAi showed no difference in in vitro growth compared to the WT, irrespective of genotypes. Under salt, alkali, or salt–alkali treatments, the growth of MdKCAB-RNAi calli was significantly reduced after 14 days of culture (Fig. 7B and C). When grown under salt, alkali, or salt–alkali stress conditions, MdKCAB overexpressed calli showed better proliferation and lower levels of MDA than the untransformed WT and MdKCAB-RNAi line (Fig. 7C and D). After 14 days of salt, alkali, or salt–alkali stress conditions, the callus transformed with 35S:MdKCAB-11G761A had a significantly higher fresh weight and lower MDA content, while other transformants had no significant difference (Fig. 7C and D). The phenotype values of SID, AID, and SAID in F1 hybrids with SNP11G/SNP761A alleles of MdKCAB were significantly lower than that of hybrids without these alleles. SNP11G761A was closely linked to the tolerance genotype of the SA-BC16.1 KASP marker (Fig. 7E, Table S23). These data indicated that SNP11G and SNP761A of MdKCAB both contributed to the improved tolerance of apple rootstock to salt, alkali, and salt–alkali stress conditions.
Changes in osmolality of Hoagland’s solution caused by the addition of NaCl and/or Na2CO3
The addition of NaCl and/or Na2CO3 in the nutrient solutions in alkaline, saline, and alkaline-saline stress treatments led to significant increases in osmolality of the solutions (Fig. S8). The osmolality values in the solutions for alkali (1.75 MPa), salt (2.73 MPa), and alkali-salt (5.12 MPa) treatments were 2, 3, and 5 times that of non-stress control (0.82 MPa) (Fig. S8).