Combining Genetic and Transcriptomic Approaches to Identify Transporter-Coding Genes as Likely Responsible for a Repeatable Salt Tolerance QTL in Citrus

The excessive accumulation of chloride (Cl−) in leaves due to salinity is frequently related to decreased yield in citrus. Two salt tolerance experiments to detect quantitative trait loci (QTLs) for leaf concentrations of Cl−, Na+, and other traits using the same reference progeny derived from the salt-tolerant Cleopatra mandarin (Citrus reshni) and the disease-resistant donor Poncirus trifoliata were performed with the aim to identify repeatable QTLs that regulate leaf Cl− (and/or Na+) exclusion across independent experiments in citrus, as well as potential candidate genes involved. A repeatable QTL controlling leaf Cl− was detected in chromosome 6 (LCl-6), where 23 potential candidate genes coding for transporters were identified using the C. clementina genome as reference. Transcriptomic analysis revealed two important candidate genes coding for a member of the nitrate transporter 1/peptide transporter family (NPF5.9) and a major facilitator superfamily (MFS) protein. Cell wall biosynthesis- and secondary metabolism-related processes appeared to play a significant role in differential gene expression in LCl-6. Six likely gene candidates were mapped in LCl-6, showing conserved synteny in C. reshni. In conclusion, markers to select beneficial Cleopatra mandarin alleles of likely candidate genes in LCl-6 to improve salt tolerance in citrus rootstock breeding programs are provided.


Introduction
Salinity in soil and irrigation water constitutes a serious global threat to food security, affecting 2.1% of dry agricultural land and up to 19.5% of irrigated land, which accounts for one-third of world food production [1].Losses in agricultural production globally due to salinity are estimated to be USD 27 billion per year [2].Citrus is ranked among the most salt-sensitive tree crops [3].Tree growth and fruit yield of citrus species are impaired at a soil salinity of approximately 2 dS/m soil saturation, without any concomitant expression of leaf symptoms [4,5].Additionally, salinity predisposes the citrus tree root to attacks by Phytophthora [6], nematodes [7], and bacterial pathogens [8].As citrus varieties of sweet oranges, mandarins, grapefruits, pummelos, and lemons are always propagated by budgrafting onto a seedling rootstock in order to obtain a uniform orchard with tolerance to soil pathogens and well-adapted to local edaphoclimatic conditions, salt tolerance is often a target trait in citrus rootstock breeding programs [9].In addition, the root (contributed by the rootstock) is the first plant organ to come into contact with salts and plays a major role in water and nutrient uptake.
However, citrus breeding for salt tolerance is not an easy task due to the following factors: (1) salt tolerance is a complex, quantitative trait that depends on the level and duration of salinity, its definition (vegetative growth, fruit yield), and the progeny (parents) under study; and (2) the long juvenility period of citrus (a mode value of 8 years; [10]) and segregation for nucellar embryony in citrus rootstock progenies [10].Fortunately, the rich citrus germplasm has great halotolerant genetic potential to improve citrus rootstocks [11][12][13][14][15][16][17].
Citrus species used as rootstocks were originally classified into three groups [18]: Citrus reshni (Cleopatra mandarin), with a good level of salt tolerance; Citrus volkameriana (Volkamer lemon) and Citrus aurantium (sour orange), with a medium level of salt tolerance; and Poncirus trifoliata (trifoliate orange), with low salt tolerance levels.The literature on the ranking of species since that time has been abundant, and the differences observed depend on the specific species and the set of accessions per species under study.However, Cleopatra mandarin has always been ranked at the top of these classifications [11,13,14], which is useful for breeding improved rootstocks [11,17].The identification of beneficial gene alleles in the halotolerant germplasm using QTL analysis of salt tolerance would facilitate their introduction into new citrus rootstocks.This would enable marker-assisted selection, make breeding programs more efficient-particularly in relation to fruit trees with long juvenility-and make the genetic modification of this complex trait feasible.
Under saline conditions, in several fruit crops, such as citrus, grapevine, avocado and persimmon, the excessive accumulation of Cl − in leaves (but not of Na + ) has been found to be related to decreased transpiration, photosynthesis, yield and quality of crops, and eventually the death of the plant [19][20][21][22][23][24].All the data indicate that the accumulation of Cl − ions in shoots is more often associated with a reduction in the growth and yield of citrus under salinity conditions than with the accumulation of Na + in shoots.However, the genetic regulation of the movement of Cl − from the root to the shoot as compared to that of Na + has been rarely studied [25].Cl − , which is an essential micronutrient for plants and acts as a cofactor in photosynthesis and in numerous enzymatic activities.However, given its mM accumulation levels in plants, Cl − is also a beneficial macronutrient and is involved in osmotic functions, affecting key processes such as stomatal movement, charge balance, and cell expansion [26,27].In fact, low concentrations (1-5 mM) of Cl − in irrigation water increase water (WUE), nitrogen (NUE), and CO 2 use efficiency in well-watered tobacco plants and improve their ability to withstand drought [26][27][28].However, under saline conditions, a reduction in Cl − concentrations in the xylem is the key step in reducing Cl − toxicity in the aerial part, which can be achieved by limiting the entry of Cl − into the xylem and/or by enhancing its removal [23,25].Genomic-wide transcriptomic information on putative Cl − transporters involved in citrus salt tolerance is limited [20,29,30], particularly regarding the root [29,30].The experimental design of these transcriptomic studies based on just one genotype [29,30] or two genotypes from different species [20] prevents distinguishing heritable tolerance responses from global responses or species-specific responses to salinity.
The detection of genomic regions containing QTLs controlling net Cl − entry into the xylem could play a useful role in the search for candidate genes coding for the transporter(s) involved.Early attempts at QTL analysis of this trait have been reported by Tozlu et al. [31] and Raga et al. [17].However, genetic studies of citrus are usually based on a small number of genotypes, which greatly affects the power of detection of QTLs, even after increasing the number of replicates per genotype by using nucellar seedlings.

Results
This study covers the results obtained from two long-lasting salt tolerance experiments using the same mapping population, named R×Pr.In one experiment, a citrus commercial variety was grafted onto the nucellar seedlings from apomictic R×Pr hybrids (grafted population, GP experiment, Figure S1A,B).In the other experiment, new nucellar nongrafted seedlings from apomictic R×Pr hybrids were used (NG experiment, Figure S1C).

Salt Tolerance in the GP Experiment
Under salinity conditions, the estimated heritabilities for some traits, such as leaf [Cl − ], fruit soluble solids content, and total fruit yield (Cl_L_S, SSC_S, and TFW_S, respectively), were substantial (higher than 0.3) in this experiment (Table S1).Salt tolerance, measured as TFW_S, was highly significantly correlated with Cl_L_S (r = 0.59; Table S2).In addition, TFW_S was positively correlated with TChl_S and SSC1_S and negatively correlated with K_L_S and Na_L_S.The distributions of relevant traits evaluated in this experiment are shown in Figure S2, in which the means of the parents (Cleopatra, Flying Dragon) and the reference rootstock Carrizo are also included.
A few significant QTLs are shown in Table 1, most of which correspond to Cleopatra segregation.A clustering of QTLs for fruit yield traits, leaf concentrations of Cl, Fe, and K under salinity conditions, and soluble solid content (SSC) of fruit juice, also under salinity conditions, was observed in linkage group 4c (C.clementina scaffold 6) near marker CR23,750 (in Ciclev10013177m.g),where a major QTL (percentage variance explained (PVE) = 52%) controlling Cl_L_S was detected (Figure S3A).The directions of correlations mostly fitted genotypic means at CR23,750, with Cleopatra allele a (in genotypes ac and ad) being associated with low Cl_L_S and high TFW, NFp, and Fe_L_S, while allele b was associated with high Cl_L_S and low fruit yield (Figure S3B).Therefore, salt tolerance conferred by the R×Pr rootstock on the satsuma variety in terms of fruit yield is genetically related to the leaf Cl − concentration, as QTLs for both traits were found to be located in the same genomic region around marker CR23,750 in this experiment.

Salt Tolerance in the Non-Grafted Population (NG Experiment)
A very wide distribution of leaf Cl concentrations under salinity conditions was observed as compared to that for the root (Figure S4A,B).In the plant, under control conditions, Cl − concentrations are generally higher in the root than in the leaf with respect to most of the hybrids; however, under salinity conditions, around 50% of the population presented higher Cl − levels in the leaf than in the root (Figure S4C), meaning an excess.
Cl − was transported to the leaves, where it then accumulated.Eight hybrids, including hybrid 107, showed no significant differences in the distribution of plant Cl − between salinity levels and showed similar leaf Cl − levels to those of the salt-tolerant parent, Cleopatra (Figure 1C).These hybrids are considered salt tolerant.Moreover, similar differences in leaf Cl − levels between salt-sensitive and tolerant genotypes (trifoliate orange versus Cleopatra, and 90 versus 107) under salinity (Figure 1A) were also found in the previous (GP) experiment (Table 2) where leaves belong to a grafted mandarin.
Table 1.List of positions in centiMorgan (cM), log of odd ratio (LODs,) as well as nearest markers or marker intervals of quantitative trait loci (QTLs) detected by IM in the GP experiment in the integrated Citrus reshni-Poncirus trifoliata genetic linkage map (LG) using the cross-pollinated model.The QTLs that were also detected in the individual parental linkage maps are indicated by putting the parental linkage groups in parentheses (R and Pr for C. reshni and P. trifoliata, respectively).The four genotypic means (ac, ad, bc, and bd, being C. reshni ab and P. trifoliata cd), as well as the percentage of explained variance (PEV), are also indicated.LOD values higher than genome-wide significant LOD scores are indicated in bold.Results from the mixed model analyses and estimated heritabilities are shown in Table S3.As expected, G×E interactions were highly significant (p < 0.0001) for Cl − -related traits, thus indicating important differences between the behaviors of hybrids when salinity levels change.Differences in phenotypic plasticity with respect to Cl − traits can be easily observed as differences in the slopes of their reaction norms, particularly for leaf [Cl − ] and the leaf-to-root distribution of Cl − concentrations, where the salt-tolerant hybrids showed almost no change (Figure S5A,C).Thus, salt tolerance is associated with stability (robustness) in leaf and leaf-to-root Cl − concentrations.

Trait
Under salinity conditions, leaf [Cl − ] is positively related to leaf [Na + ] and [Ca 2+ ] and negatively related to N_L_S (Table S4).In general, root mass, total root dry weight (TRDW), and fine root dry weight (FRDW) are positively related to leaf [Cl − ], meaning that the larger the root, the more Cl − is accumulated in the leaf.It is noteworthy that, under salinity conditions, TRDW and FRDW are additionally related (negatively) to root [Cl − ]; this means that the larger the root, the more Cl − is accumulated in the leaf and the less Cl − is accumulated in the root.FRDW increased under salinity conditions, an increase which was particularly high in some hybrids (Figures S5E and S6), making the G×E interaction significant for this trait (Table S3).S1 and S2 for GP and NG experiments, respectively.0.8 ± 0.0 1.0 ± 0.1 0.9 ± 0.1 0.8 ± 0.0 0.9 ± 0.0 Ca_L A list of significant QTLs is presented in Table 3.Most QTLs corresponded to the segregation of Cleopatra alleles, with, again, a clustering of QTLs being observed in linkage group 4c (Figure S7).Again, the QTL for leaf [Cl − ], henceforth named LCl-6, which had a particularly strong effect (39.4%) under salinity conditions, was located in the genomic region between markers 15R,750 at Ciclev10011720m and CR28,270 at Ciclev10011175m in C. clementina scaffold 6.The QTLs influencing root dry mass (TDRW_C and FRDW_C) and leaf Ca, C, and Fe (Ca_L_S, C_L_S, and Fe_L_S) were also present in this genomic region.Salt tolerance in terms of leaf Cl content was again associated with the Cleopatra allele a at marker CR23,750 (genotypes ac and ad in Figure S8).
The root transcriptomics of two full-sibs from the R×Pr population, which differed considerably both in terms of salt tolerance (Figures 1 and S6) and of the Cleopatra allele at LCl-6 (numbered 107 and 90, salinity-tolerant and -sensitive, respectively), were comparatively studied by bioinformatically analyzing the RNA sequencing in the NG experiment after 383 days of salt treatment (Table S6).Taking into account the whole genome of C. clementina as reference, around 80% of genes were expressed at the root under our experimental conditions; in addition, the number of differentially expressed genes (DEGs) depended on the method used for their assessment (iDEP versus RSeqFlow with the use of an eBayes or treat procedure) and the comparison under consideration as shown in Table 5.Most DEGs corresponded to the comparison between hybrids and their responses to salinity ((107_15 + 107_0) vs. (90_15 + 90_0)), followed by the comparison between hybrids under control conditions.When considering only the genes included in this QTL (452 or 465, depending on the reference genome), 322 (71%) or 303 (65%) genes were expressed at the root.From the 322 QTL-expressed genes (C. clementina as reference and RSeqFlow together with Treat as the procedure used), the comparison between 107_15 and 90_15 resulted in eight DEGs, the comparison between 107_0 and 90_0 resulted in 21 DEGs, while the comparison between (107_15 + 107_0) and (90_15 + 90_0) resulted in 92 DEGs.It is worth noting that no differential transcriptomic root response to salinity was hybrid-independent ((107_15 − 90_15) vs. (107_0 − 90_0)) at either the whole genome or LCl-6 QTL level.
Table 5. Number of differentially expressed genes (DEGs) per comparison, depending on the method used for their assessment, iDEP versus RSeqFlow, using the eBayes or Treat procedure, for the whole root transcriptome or for genes in the LCl-6 QTL only.The C. clementina genome was used as reference for RNAseq readings here, and the Treat procedure was chosen for further analyses of DEGs.A total of 78 genes (68 annotated genes) in the QTL were involved in these DEGs as a whole.All DEGs and the corresponding genes in the salt tolerance QTL are listed in Table S7.Thus, from our initial list of 23 functional candidate genes in Table 4, 11 became relevant when their differential mRNA expression at the root was tested; of these 11 genes, 4 are NPFcoding genes: Ciclev10013337m (NPF5.9),Ciclev10011381m (NPF8.1),Ciclev10013821m (NPF8.2),and Ciclev10011341m (NPF5.12).It is also worth noting that the only candidate gene whose expression in the tolerant hybrid decreased was Ciclev10011341m (NPF5.12).

Comparison
Considering that P. trifoliata is the reference genome, some candidate genes in the QTL also showed differential mRNA expression (Table S8).Most Citrus and Poncirus DEGs correspond to the same candidate genes except those coding for PIP2.1.Thus, gene Ptrif.0006s0418.2 (ortholog of Ciclev10012379m) is differentially expressed only when Poncirus is considered to be the reference genome, and gene Ciclev10012375m is differentially expressed when C. clementina is the reference genome.The C. clementina ortholog of gene Ptrif.0006s0399(coding for a reverse transcriptase), which corresponds to a DEG in two comparisons, was absent; this suggests that differences in Poncirus-specific retrotransposon activity exist between the hybrids in response to salinity at the root.
Once gene alleles a/b from Cleopatra and c/d from Flying Dragon could be assigned by whole genome sequencing of trees 107, 90, and 22-7, whose genotypes at CR23,750 were ac, bc, and bb, respectively, it was possible to infer the translated sequence of each allele of the differentially expressed candidate genes.The alignments of translated sequences for the alleles of Ciclev10013337m/Ptrif.0006s0823.1 (NPF5.9),Ciclev10012375m/Ptrif.0006s0419.1, coding for plasma membrane aquaporin PIP2.1, and Ciclev10011745m/Ptrif.0006s0610,coding for a major facilitator superfamily (MFS) protein, are shown in Figures S9, S10, and S11, respectively.
The agglomerative hierarchical clustering of DEGs in LCl-6 based on their expression patterns, represented as the normalized mean of counts per million reads mapped, resulted in three clusters (Figure S13).Genes in cluster AHC-1 exhibited a slight increase in expression under salt stress conditions in genotype 107, while decreasing in genotype 90 under similar conditions.The co-expression network of DEGs in AHC-1 is shown in Figure 2. Thus, one can clearly observe several processes usually related to plant development and en-vironmental stress, such as cell wall biosynthesis mechanisms, specifically glucuronoxylan and xylan biosynthesis, as well as lignin metabolism, in addition to secondary metabolism (alkaloid/phytoalexin biosynthesis).
The agglomerative hierarchical clustering of DEGs in LCl-6 based on their expression patterns, represented as the normalized mean of counts per million reads mapped, resulted in three clusters (Figure S13).Genes in cluster AHC-1 exhibited a slight increase in expression under salt stress conditions in genotype 107, while decreasing in genotype 90 under similar conditions.The co-expression network of DEGs in AHC-1 is shown in Fig- ure 2. Thus, one can clearly observe several processes usually related to plant development and environmental stress, such as cell wall biosynthesis mechanisms, specifically glucuronoxylan and xylan biosynthesis, as well as lignin metabolism, in addition to secondary metabolism (alkaloid/phytoalexin biosynthesis).S13) in the salt tolerance quantitative trait locus (QTL) LCl-6.Names of transcripts and genes are provided, together with their putative function, showing that processes related to cell wall biosynthesis (specifically, glucuronoxylan and xylan biosynthesis, as well as lignin metabolism) and secondary metabolism (alkaloid/phytoalexin biosynthesis) at the root could be involved in salt tolerance.

Incorporating LCl-6 Candidate Genes into the Linkage Maps to Improve the Mapping Resolution of QTL Analysis
The analysis of DEGs was used as a tool to rank the different functional candidate genes within LCl-6.Thus, the whole reference mapping population R×Pr was genotyped for the most relevant candidates, NPF5.9, MFS, NPF8.1, PIP2.1, and ABCG.It is worth noting that Flying Dragon and Rich were genotypically different regarding candidate gene NPF5.9, with just one allele in common.All the candidate genes mapped in linkage group 4c were observed to maintain C. clementina synteny (Figure S14).All these candidate genes, except NPF8.1, which only segregated for the Poncirus alleles, were mapped in C. reshnii linkage group R4c (Figure S14B).QTL analyses of chloride-related traits evaluated in GP and NG experiments, and that previously reported [17], showed that the closest gene to the LOD peak was usually MFS by using both the cross-pollination genotyping configuration and the double haploid configuration for the female (Cleopatra) segregation (Figure 3A,B, respectively).In the case of the non-grafted experiment, leaf [Cl − ], and the difference between root and leaf chloride concentrations relative to that of the root under salinity, (ClR-ClL/ClR)_S, LOD peaks were flatter and fell between MFS and ABCG.It is worth noting that LOD scores at QTL peaks of these Cl-related traits increased when the candidate genes were incorporated into the linkage map.When Kruskal-Wallis statistic K was used to study marker-trait associations, similar results were obtained (Figure 3C).In the case of (ClR-ClL/ClR)_S, the K values for MFS, ABCG, and PIP2.1 were similar (16.486, 16.327, and 16.707, respectively).Log of odd ratio (LOD) profiles of quantitative trait loci (QTLs) for Cl-related traits in linkage group 4c (A) and R4c (B) after incorporating candidate genes into the linkage maps and using Kruskal-Wallis statistic K for the association between phenotype and Cleopatra marker/gene segregation (C).The red horizontal line indicates significant levels.Cl*_S corresponds to the salt tolerance QTL reported by [17].To distinguish between experiments, the prefix GP or NG has been added to the trait name (Cl_L_S).
Figure 3. Log of odd ratio (LOD) profiles of quantitative trait loci (QTLs) for Cl-related traits in linkage group 4c (A) and R4c (B) after incorporating candidate genes into the linkage maps and using Kruskal-Wallis statistic K for the association between phenotype and Cleopatra marker/gene segregation (C).The red horizontal line indicates significant levels.Cl*_S corresponds to the salt tolerance QTL reported by [17].To distinguish between experiments, the prefix GP or NG has been added to the trait name (Cl_L_S).
The resolution of QTL analysis of fruit yield (TFW, NFp), root mass (TRDW, FRDW) under control conditions, as well as leaf Fe concentration under salinity conditions (GP_L_Fe_S), was improved by incorporating the candidate genes into the linkage map (Figure S15).Thus, the LOD profile for GP_L_Fe_S reached a maximum exactly at NPF5.9, where TFW_S and NFp_S also displayed one of two peaks.Similarly, root growth under control conditions clearly peaked around ABCG.On the other hand, other QTLs reported in linkage group 4c, such as L_Ca_S and L_K_S in the GP experiment (Table 1), became non-significant.
From a practical perspective, markers at MFS and NPF5.9 (Figure 4) constitute new tools to help select beneficial Cleopatra alleles in rootstock breeding programs that use Cleopatra progenies to improve the salt tolerance of citrus plants.

The Salt Tolerance QTL LCl-6 Has Been Consistently Detected across Three Experiments
There have been two previously reported QTL analyses of salt tolerance in citrus plants [17,31].The latter study used the same R×Pr reference population of hybrids as the GP and NG experiments in this study.It is worth noting that [17] also detected a QTL for leaf [Cl − ] in the same region as LCl-6, around the SSR marker CR23,750, where Cleopatra allele segregation was also involved (Cl*_S in Figure S3).Thus, the salt tolerance QTL LCl-6 has been consistently detected across three independent experiments (one grafted and two non-grafted experiments) using the same R×Pr reference population derived from the intergeneric cross between Cleopatra mandarin (R) and Poncirus trifoliata (Pr).Therefore, LCl-6 can be considered a confirmed QTL [32], whose responsible genes are acting at the root, which is in line with previous reciprocal grafting studies of citrus plants [33].
Although other QTLs for traits correlated with leaf [Cl − ] (Tables S2 and S4) were also detected in the same genomic region as LCl-6, their detection generally depended on the experiment carried out (Figures S3 and S7).For example, a QTL for leaf K + concentration under salinity conditions (K_L_S; Figure S3) was detected in the GP experiment and was also reported by [17]; QTLs for Fe_L_S and Ca_L_S were detected in both the GP (Figure S3) and NG (Figure S7) experiments.In addition to Cl-related traits, the LCL-6 region contained QTLs for mandarin fruit yield (TFW, NFp) in the GP experiment (Figure S3) and for root growth (TRDW_C and FRDW_C) in the NG experiment (Figure S7).Despite the significant correlation between Cl_L_S and N_L_S (r = −0.5,Table S4), suggesting a certain degree of competition between Cl − and NO3 − under excess soil Cl − conditions, no QTL for N_L_S was detected at or close to LCl-6.
Regarding other QTL analyses of stress tolerance traits, LCl-6 overlaps with CD/FS-2016-t6, two QTLs for Huanglongbing tolerance, in terms of canopy damage and foliar symptoms, mapped to the Poncirus trifoliata genome [34].In addition, these tolerance

The Salt Tolerance QTL LCl-6 Has Been Consistently Detected across Three Experiments
There have been two previously reported QTL analyses of salt tolerance in citrus plants [17,31].The latter study used the same R×Pr reference population of hybrids as the GP and NG experiments in this study.It is worth noting that [17] also detected a QTL for leaf [Cl − ] in the same region as LCl-6, around the SSR marker CR23,750, where Cleopatra allele segregation was also involved (Cl*_S in Figure S3).Thus, the salt tolerance QTL LCl-6 has been consistently detected across three independent experiments (one grafted and two non-grafted experiments) using the same R×Pr reference population derived from the intergeneric cross between Cleopatra mandarin (R) and Poncirus trifoliata (Pr).Therefore, LCl-6 can be considered a confirmed QTL [32], whose responsible genes are acting at the root, which is in line with previous reciprocal grafting studies of citrus plants [33].
Although other QTLs for traits correlated with leaf [Cl − ] (Tables S2 and S4) were also detected in the same genomic region as LCl-6, their detection generally depended on the experiment carried out (Figures S3 and S7).For example, a QTL for leaf K + concentration under salinity conditions (K_L_S; Figure S3) was detected in the GP experiment and was also reported by [17]; QTLs for Fe_L_S and Ca_L_S were detected in both the GP (Figure S3) and NG (Figure S7) experiments.In addition to Cl-related traits, the LCL-6 region contained QTLs for mandarin fruit yield (TFW, NFp) in the GP experiment (Figure S3) and for root growth (TRDW_C and FRDW_C) in the NG experiment (Figure S7).Despite the significant correlation between Cl_L_S and N_L_S (r = −0.5,Table S4), suggesting a certain degree of competition between Cl − and NO 3 − under excess soil Cl − conditions, no QTL for N_L_S was detected at or close to LCl-6.
Regarding other QTL analyses of stress tolerance traits, LCl-6 overlaps with CD/FS-2016-t6, two QTLs for Huanglongbing tolerance, in terms of canopy damage and foliar symptoms, mapped to the Poncirus trifoliata genome [34].In addition, these tolerance QTLs included candidate gene Ptrif0006s0773.1,coding for a laccase 7-related protein and responsive to infection by Candidatus Liberibacter, the Huanglongbing pathogen [35], which was also expressed in the root (ESR42493) in LCl-6 in our salinity experiment.The overlapping of the salt tolerance QTL LCl_6 and these Huanglongbing tolerance QTLs genetically supports the relationship between salinity and bacterial pathogen effects on citrus plants [8].
Therefore, the genomic region where LCl-6 is located appears to contain genes that control relevant agronomic traits, ranging from root development to leaf nutrient content and fruit yield-related traits mediated by the rootstock, which merit further molecular research.

Salt Tolerance Mechanism(s) behind LCL-6
Both leaf Na + and Cl − concentrations were studied in the two salt tolerance experiments described here.Under salinity conditions, both parameters were positively correlated and negatively related to fruit yield in the GP experiment.Only a leaf Cl − QTL (LCl-6) was consistently detected under salinity conditions in the same genomic region in independent experiments.Moreover, this region also contained a fruit yield QTL that genetically supports the correlation between leaf Cl − and fruit yield in the GP experiment under salinity conditions.These results in the R×Pr population are in line with those from numerous scientific studies, which conclude that Cl − toxicity is more important than Na + toxicity in citrus plants and other woody species under salinity conditions [12,[36][37][38][39].The differences in Cl − tolerance exhibited by plants are usually related to the ability to restrict Cl − transport.This Cl − transport in long-term stressed plants involves two main steps: (1) a net influx (from the influx/efflux balance) of Cl − from the soil to the root and (2) a net loading (from the loading/unloading balance) of Cl − onto the xylem to follow the transpiration stream up to the leaves.Regarding the comparison of the behavior of the R×Pr population with respect to the distribution of Cl − between leaf and root (Figures 1C and S5C), although salinity increases Cl − in both organs, this increment is proportionally smaller in the leaves of salt tolerance hybrids like 107 (Figure 1C).Thus, while hybrids 107 and 90 have similar root Cl − levels (Figure 1B), they greatly differ with respect to leaf Cl − concentrations under salinity conditions (Figure 1A), resulting in practically no change in Cl − organ distribution between treatments in salt-tolerant hybrids (Figures 1C and S5C).Salt-sensitive hybrid 90 accumulates Cl − in the leaves, while salt-tolerant hybrid 107 shows higher Cl − levels in the root than in the leaves under salinity conditions, and both hybrids show similar root Cl − levels.It, therefore, seems reasonable to hypothesize that the main mechanism involved here to prevent leaf Cl − accumulation takes place in the first step (net influx of Cl − from soil to the root).A greater influx of Cl − into hybrid 90 would facilitate higher leaf Cl − levels in 90, and a greater efflux of Cl − out of hybrid 107 (such as reported in poplar species [38]) would explain the salt tolerance behavior of hybrid 107.However, taking into account the differences in root growth between hybrids 107 and 90 under salinity conditions (Figure S6), a simpler hypothesis, presented below, could explain the higher leaf Cl − levels in hybrid 90 than in 107, as citrus leaf Cl − accumulation has been reported to be linked to water use [40].
The root, which is the first plant organ to sense salinity in the soil, plays an important role in water uptake; therefore, its growth response to salinity has to be included in our hypothesis regarding leaf Cl − exclusion.In general, although salinity reduces root mass [41], root growth responses to salinity reported in citrus-related literature are contradictory [29,[42][43][44]; this is in line with our results regarding G×E interactions (Table S3) and reaction norms for fine root mass (Figure S5E).In our experiment, the root mass of some R×Pr hybrids, particularly that of salt-sensitive hybrid 90, increased considerably under salinity conditions (Figure S6).Since root dry mass (TRDW and FRDW) and leaf [Cl − ] (Cl_L) are positively correlated (Table S4), it is reasonable to hypothesize that, under salinity conditions, the larger root of hybrid 90 could simply facilitate a higher net influx of Cl − from the soil to the root, and finally to the shoot, than that for hybrid 107.
Co-expression approaches have previously been used to assign function to genes involved in root elongation and other related traits because genes functioning in the same pathway or required for the same process tend to express in a transcriptionally coordinated manner [45].Thus, following the agglomerative hierarchical clustering (AHC) analysis of DEGs in LCl-6, the co-expression network of DEGs in AHC-1 led to a slight increase in gene expression in hybrid 107 but to a decrease in hybrid 90 under salinity conditions (Figure S13).This shows that several processes related to plant development and in response to environmental stress, such as cell wall biosynthesis mechanisms, specifically glucuronoxylan and xylan biosynthesis, as well as lignin metabolism and secondary metabolism (alkaloid/phytoalexin biosynthesis), could be involved in the salt tolerance determined by LCl-6 (Figure 2).The effect of these underlying biosynthesis processes in this QTL region on the differences in root growth under salinity conditions between both hybrids needs to be further investigated.It would not be the first time that DEGs involved in cell wall loosening are associated with reduced root growth in a salt-tolerant citrus cultivar [29].Additionally, increased cell wall lignification in the root of the salt-tolerant hybrid could render the apoplastic movement of Cl − more difficult.

Salt Tolerance Candidate Genes Underlying LCl-6 in the R×Pr Population
Numerous genes (23) coding for aquaporins and transporters, in general, were detected in LCl-6 (Table 4).Some of these genes could be involved in the movement of Cl − from root to shoot, as has previously been reported regarding PIP1 and PIP2 coding genes in citrus plants [46,47].Thus, differences between the genotypes Cleopatra (C.reshni), Carrizo (a C. sinensis × P. trifoliata hybrid), and P. trifoliata in PIP1 expression were found to be related to Cl − exclusion from leaves, probably due to the effects on water movement, although salinity did not affect its gene expression [46].Taking into account the primer sequences provided [46], an isomorph could only be assigned to PIP1 (PIP1.4),corresponding to Ciclev10012384m, which is included in LCl-6 (Table 4), but with no differential expression between hybrids 107 and 90.In Arabidopsis, PIP2.1 substrate transport activity can be switched between ion and water channel models using phosphorylation [48].Some NPF proteins have been reported in previous studies to have a major impact on Cl − homeostasis in citrus leaves (NRT1-2; [20] and Arabidopsis roots (NPF2.5;[23]).None of the citrus candidate genes identified by Brumós et al. [21] for Cl − homeostasis are located in LCl-6 or even in chromosome 6.The rootstock root is the organ of the citrus plant that contributes to the uptake of water and elements from the soil, with these genes coding for transporters belonging to large gene families.The next step, to rank the salt-tolerant candidate genes in the Cleopatra mandarin genome, involved studying root differential gene expression in LCl-6 using hybrids 90 and 107, two salt tolerance-contrasting full-sibs of the R×Pr population, which share the Poncirus allele but differ in relation to the Cleopatra allele in the LCl-6 region.Thus, 11 genes coding for transporters, putatively involved in ion homeostasis, showed differential expression in LCl-6, with C. clementina and P. trifoliata used as reference genomes (Table 4 and Table S5, respectively).Four of the seven NPF coding genes located in LCl-6 exhibited differential expression in at least one of the comparisons carried out: coding genes for NPF5.12 and NPF8.1 in two comparisons and coding for NPF5.9 in three comparisons (Table 4).Noteworthily, Zhao et al. [30] found that salinity induced the differential expression of these NPF coding genes in the roots of Poncirus trifoliata, here named Ptrif.0006s0823(Pt6g013450), Ptrif.0006s2462(Pt6g014250), Ptrif.0006s0814(Pt6g013550), and Ptrif.0006s0815(Pt6g013550) (Table S5).One member of them, coding for NPF5.12, is the only gene whose expression level in the root is higher in the saltsensitive genotype than in the salt-tolerant genotype (Tables S7 and S8), which could be related to a higher influx of Cl − from soil to root in the salt-sensitive genotype.In Arabidopsis, He et al. [49] found that NPF5.12, which is a vacuolar nitrate efflux transporter, is preferentially expressed in root pericycle cells and xylem parenchyma cells and plays an important role in modulating the allocation of nitrate among roots and shoots.In Brassica rapa, NPF5.12 was found to be upregulated in roots under low nitrate conditions, suggesting that it plays a positive role in nitrate absorption [50].In addition to NPF5.9, another gene coding for a transporter and showing differential expression in most of the comparisons is Ciclev10011745/Ptrif.006s0610, coding for a major facilitator superfamily (MFS) protein (Tables 4 and S5).In Arabidopsis, its ortholog regulates salt, drought, heat stress, and turgor-dependent growth through the ABA-dependent signal transduction pathway [51].The differential expression of these genes has not been validated by qPCR because it is documented [52,53] that results obtained with RNA-seq correlate very consistently with qPCR, and the drastic differences in the expression of genes coding for NPF5.9 and MSF (Figure S12) make them robust enough to not require validation by qPCR or other approaches [54].
The incorporation of relevant candidate genes into the genetic maps improved the resolution of the QTL analyses, particularly for leaf Fe accumulation under salinity conditions in the GP experiment (comparing Figures S15 and S3) and root growth under control conditions in the NG experiment (comparing Figures S15 and S7); this points to NPF5.9 and the ABCG transporter coding genes as the most appropriate candidate genes, respectively.Thus, the fine mapping of leaf [Fe] in the GP experiment under salinity conditions at exactly NPF5.9 (Figure S15) concurs with the results reported by Chen et al. [55], who found that NPF5.9 and NPF5.8 mediate iron (and nitrate) long-distance transport and homeostasis in Arabidopsis.
In the case of leaf Cl − accumulation across the three long-lasting salinity experiments (Figure 3), interval mapping and Kruskal-Wallis analyses point to MFS as the most significant candidate gene regarding Cl*_S (leaf [Cl − ] in the experiment reported by Raga et al. [17]), leaf Cl − accumulation in the GP experiment (GP_L_Cl_S), and the phenotypic plasticity of leaf [Cl − ] (dCl_L) in the NG experiment.However, given the shape of the LOD and K profiles, the presence of more than one closely linked candidate gene controlling these leaf Cl-related traits cannot be ruled out, and the ability of MFS and NPF candidate genes to transport Cl − needs to be tested in future experiments.In addition, it must also be noted that phenotypic differences might not be related to differences at the mRNA level but at the protein level (protein abundance, amino acid sequence, etc.).
An important achievement of the present study is the biotechnology developed to speed up and increase the efficiency of rootstock breeding programs in order to confer salt tolerance on citrus varieties as a strategy to mitigate the effects of climate change on the availability of good-quality water.Thus, the markers developed for candidate genes MFS and NPF5.9 facilitate the identification of the Cleopatra allele that increases salt tolerance in rootstock breeding programs, in which Cleopatra has been used as a salt-tolerance donor (Figure 4), thus making marker-assisted selection for this quantitative trait feasible.

Plant Materials
The R×Pr mapping population, consisting of 151 hybrids obtained at IVIA in Spain, is composed of controlled crosses between Citrus reshni Hort.ex.Tan.(Cleopatra mandarin) as female parent (salt and iron chlorosis-tolerant and apomictic) and two apomictic and disease-resistant varieties of Poncirus trifoliata (L.) Raf.(trifoliate orange, Pr): 83 Flying Dragon hybrids and 68 Rich hybrids as pollinators [10].Seedlings from mature fruit yielding R×Pr trees were analyzed using molecular markers in order to discard the zygotic seedlings [56] to obtain nucellar plants for the GP (grafted population) and NG (non-grafted) experiments.
The procedures and steps involved in obtaining both mapping and grafted populations were previously described [57].In the GP experiment (2013-2014), six nucellar seedlings, obtained from each of the 62 R×Pr hybrids, showing apomictic reproduction and parents (Cleopatra and Flying Dragon), were grafted with Clausellina mandarin (Citrus unshiu (Mak.)Marc.) and maintained for several years until full production prior to the salinity experiment at IVIA (Valencia, Spain).In the NG experiment (2020), 12 nucellar seedlings, obtained from each of the 42 apomictic R×Pr hybrids and parents (Cleopatra and Rich) in 2018, were grown at the Estación Experimental La Mayora (IHSM-CSIC), Malaga, Spain.
In the GP experiment, two or three trees from the six repetitions (nucellar grafted plants) from each R×Pr hybrid were randomly assigned to the control and salinity treatments from 15 June 2013 to 20 October 2014 (16 months) in a greenhouse (Figures S1B and S1A, respectively).Plants were grown in pots (17 L) using cocofiber as a substrate.The greenhouse had an automatic roof ventilation and heating system to maintain the interior air temperature above 8 • C. A high-frequency fertirrigation system, together with 4 L h −1 drippers, was used and regulated to ensure the homogeneity of the nutrient solution in the roots of all plants cultivated simultaneously.The nutrient solution (pH: 6.4) contained the following concentrations of macronutrients (in mM) NO 3 − 8.1, H 2 PO 4 − 4, SO 4 2− 1, NH 4 + 0.9, K + 4.2, Ca 2+ 3.5, and Mg 2+ 1, in addition to the following concentrations of micronutrients (in µM): Mn 2+ 8, Zn 2+ 2.3, B 20, Cu 2+ 7, Mo 4+ 0.5, and Fe 2+ 15.3.The water for the nutrient solution was previously subjected to reverse osmosis treatment.Although the salinity-treated plants were similarly irrigated, the nutrient solution was supplemented with 30 mM NaCl.
In the NG experiment (Figure S1C), six plants from the twelve repetitions (eightmonth-old nucellar seedlings) of each apomictic R×Pr hybrid and parent were randomly selected to be subjected to control and salinity treatments over a period of a year from 1 October 2019 to 1 October 2020 in a greenhouse under natural light conditions with no temperature control.Plants were grown under natural greenhouse conditions in pots (4 L) using a vermiculite substrate.The nutrient solution (electrical conductivity:1.68,pH: 7.13) contained the following concentrations of macronutrients (in mM) NO 3 − 8.45, H 2 PO 4 − 0.74, SO 4 2− 1.84, K + 6.08, Ca 2+ 4.25, and Mg 2+ 1.33, in addition to the following concentrations of micronutrients (in µM) Fe 2+ 66, Mn 2+ 1.33, Zn 2+ 2.3, B 17, and Cu 2+ 2.3, supplemented with either 0 mM NaCl for control (1.67 dS m −1 ) or 15 mM NaCl for the saline treatment (4.67 dS m −1 ) until completion of the experiment.Plants were watered automatically using 2.5 L/h drippers three times per week, on alternate days, receiving 200 mL at each irrigation event.

Trait Evaluation
For the GP experiment, several physiological (Na, Ca, K, Fe, and Cl concentrations in mature leaves) and agronomic (related to fruit yield and quality) traits of the grafted satsuma variety were evaluated (see Table S1 for abbreviations) under both control and salinity conditions, denoted by the suffixes_C and _S, respectively.Three fully developed leaves per plant were sampled from vegetative spring shoots after 15 months of treatment in order to measure total chlorophyll leaf concentration (TChl, µmol m −2 ) using an SPAD-502 Plus chlorophyll meter (Konica Minolta Inc., Tokyo, Japan) and the function TChl = [(0.2861×SPADvalue) − 6.9501].Dry tissue samples of these leaves, maintained at 80 • C for 3 days, were prepared for mineral analysis by digestion in a HNO 3 :HClO 4 (2:1, v/v) solution.Inorganic solutes were determined in parts per million (ppm) (mg Kg −1 ) using the ICAP 6500 DUO/IRIS Intrepid II XDL spectrometer at the Segura Center for Soil Science and Applied Biology ionomics facility (CEBAS-CSIC; Murcia, Spain).Foliar concentrations of Cl − (milligrams per liter) were also evaluated using a chloride analyzer (Model 926, Sherwood Scientific, Cambridge, UK) and the methodology described by Gilliam [58].
Fruit yield was evaluated in terms of the number of normal, ripe fruit (FN) and total fruit weight (TFW, in g).Original trait FN was analyzed using a Poisson distribution (FNp) for QTL analysis.A minimum of five randomly sampled fruits per tree were also evaluated for the following internal fruit-quality traits at two harvest times (1: mid-September 2014; 2: mid-October 2014): fruit weight (FW, in grams), juice titratable acidity (A, %), and soluble-solids content (SSC, as • Brix) using a Pallete PR-101 digital refractometer (Atago, Tokyo, Japan).
At the end of the NG experiment, after one year of treatment, 8-10 leaves from each plant were harvested and dried to determine Cl, Na, K, Ca, and Fe, as described for the GP experiment (see Table S3 for abbreviations).Total C and N leaf concentrations (g/100 g of dry weight) were also determined using a Leco TruSpec CN628 elemental analyzer at the ionomics facility (CEBAS-CSIC).The roots of each plant were exhaustively washed with tap water and dried.Total root dry weight (TRDW, in grams) was estimated for each plant.Fine-root mass per plant was also evaluated by weight (FRDW, in grams).Cl, Na, K, Ca, Fe, C, and N were determined in the fine roots of each plant.To characterize the distribution of Cl and, similarly, of Ca in the plant, the function ([Cl] root − [Cl] leaf )/[Cl] root was calculated and used as an additional trait (ClR-ClL/ClR) at both salinity levels.The phenotypic plasticity of leaf and root Cl concentrations (dCl_L and dCl_R, respectively) and other traits were estimated as the difference between salinity and control means relative to the mean under control conditions [59].The responses of genotypes to salinity, also viewed as phenotypic plasticity, regarding chloride-related traits and root mass, were depicted as reaction norms [60].

Statistical Analysis
The GP experiment used a split-plot design with four blocks, using NaCl treatments as the main plot and rootstocks as the subplots.With regard to the statistical analysis of the experiment, the blocks were random.On the other hand, to study the genotype (G) and salinity (E) effects, as well as the G×E interactions of the evaluated traits, the effects of genotype and treatment were classed as fixed.In the NG experiment, genotypes were distributed at random within each NaCl treatment.Considering R×Pr hybrid genotypes as a random effect factor, broad-sense heritability (H 2 ) was estimated for all traits with respect to nucellar rootstocks (repetitions) derived from apomictic R×Pr hybrids under control or salinity conditions based on the genotypic (V G ) and environmental (V E ) variance estimators calculated using the minimum variance quadratic unbiased estimator (MIVQUE), as described elsewhere [61].
Spearman's correlation coefficients were used to study the relationships between the different traits.
4.4.Genetic Analyses 4.4.1.Linkage Map and QTL Analyses QTL analyses were carried out using the genotypic and map data from [10] based on SSR, IRAP, and SCAR markers, as well as the adjusted means of traits.The interval mapping (IM) procedure and the Kruskal-Wallis rank sum test (nonparametric QTL approach), together with the MapQTL ® 6 program [62], were used to identify QTLs.QTL analyses were carried out in two different ways.Firstly, we analyzed the data as a cross-pollinated (CP) population type in order to examine intralocus interactions.Secondly, using a two-way pseudo-test cross approach, we analyzed data for each parental meiosis separately [63].Although this second approach takes advantage of the computational benefits of the twogenotype QTL model, intralocus interactions were ignored, rendering this approach less powerful and realistic [62].JoinMap 4.1 software [64] was used to translate and split the marker data in order to separate the two meiosis.Being parent-specific, some linkage groups or linkage group parts (R9a, R6, R4a, Pr1, Pr4a, and Pr9b) were ignored when using the CP data for QTL analysis.The Cleopatra map contains 86 markers distributed among 10 linkage groups, covering 1127.127cM of the C. reshni genome.Similarly, the Poncirus map contains 73 markers distributed among 11 linkage groups, covering 1416.759cM of the Poncirus trifoliata genome.No genotypic differences between the male parents, Flying Dragon and Rich, were detected with respect to these mapped markers, which is in line with the high degree of relatedness found among the P. trifoliata accessions [34,65].The CP map contains 93 markers, spread among nine linkage groups, covering 1406.761cM of the integrated genome.
With respect to IM, experiment-wise significance was assessed to be 5% using 1000 permutation tests.These LOD critical values ranged from 1.1 to 2.0, depending on the specific trait and linkage group in the two-way pseudo-testcross analysis (population type DH).On the other hand, the critical LOD values ranged from 2.1 to 3.4, depending on the specific trait and linkage group in the CP analysis.Only significant QTLs with LOD ≥ 2.36 for heritable traits (H 2 > 0) are reported here.With regard to the Kruskal-Wallis procedure, the significance level for individual tests was fixed at 0.005, as recommended in the manual [62].

Candidate Genes and Linkage Analysis
A genomic region around the CR23,750 marker on the integrated linkage group 4c was found to be particularly rich in QTLs and to have QTLs for Cl-related traits.With respect to this region, markers from the CP map were anchored to the physical map of C. clementina using primer and/or EST sequences, as well as the BLASTN tool (https: //phytozome.jgi.doe.gov/pz/portal.html#!search?show=BLAST accessed on 18 February 2019).All genes from marker 15R,750 (at Ciclev10011720m.g)down to marker CR28,270 (at Ciclev10011175m.g) on integrated linkage group 4c were downloaded from both the C. clementina and Poncirus trifoliata genomes at https://phytozome.jgi.doe.gov.C. clementina was chosen instead of C. sinensis because C. reshni (Cleopatra mandarin) is genetically closer to C. clementina than to C. sinensis [66,67].The annotation of some genes downloaded from https://phytozome.jgi.doe.gov(accessed on 18 February 2019) was tested by blasting their peptide sequence at https://ncbi.nlm.nih.gov(accessed on 18 February 2019).Gene ontology (GO) enrichment analysis of genes was carried out using the Singular Enrichment Analysis tool [68] on the AgriGo platform (http://systemsbiology.cau.edu.cn/agriGOv2/accessed on 20 February 2019).
Aligned DNA sequences of relevant candidate genes were used to find divergent regions in order to design primers that reveal insertion/deletion (InDel) polymorphisms in the R×Pr population.Forward and reverse primers covering a divergent region for each candidate gene were designed using https://bioinfo.ut.ee/primer3/ (accessed on 20 February 2019) in order to develop sequence-characterized amplified region (SCAR) markers in or near these genes to examine their genetic location and quantitative effects on traits (Table S9).Genomic DNA extractions from leaf tissue were carried out for each hybrid.Polymerase chain reaction (PCR) conditions were specific to each marker, and the resulting product was analyzed by electrophoresis in 10% DNA sequencing-type polyacrylamide gels and was revealed by silver staining.All procedures are described elsewhere [56].PCR products from homozygous plants were sequenced at the López-Neyra Institute of Parasitology and Medicine DNA sequencing and genomics facility (Granada, Spain) to test their gene identification/physical location (Table S10).
JoinMap 4.1 mapping software [64] was used for segregation and linkage analysis of candidate genes in the R×Pr population (151 hybrids) using a linkage criterion LOD ≥ 10, a recombination fraction of 0.5, and the Kosambi mapping function.The CP population type, resulting from a cross between two heterogeneously heterozygous and homozygous diploid parents, was selected to analyze R×Pr progeny with no previous knowledge of the marker linkage phase.This strategy uses nn × np or lm × ll codifications for dominant markers that segregated into either P. trifoliata or C. reshni, respectively, in addition to ab × cd, ef × eg, and hk × hk codifications for codominant marker loci, at which 4, 3, and 2 different alleles are segregating, respectively.

DNA Sequencing of Selected Materials
For DNA and RNA sequencing, we selected two full-sibs of the R×Pr population, which differed in the Cleopatra allele at marker CR23,750, and the degree of salt tolerance measured in terms of both fruit yield and leaf [Cl − ] under salinity conditions in the GP experiment.These hybrids, which share the male parent Flying Dragon at CR23,750, were: hybrid 107 (salt tolerant with an ac genotype at QTL marker CR23,750) and hybrid 90 (salt sensitive with a bc genotype at QTL marker CR23,750).Additionally, using marker analysis, a homozygote (bb) at CR23,750 (22-7) was selected from a progeny obtained by self-pollination from monoembryonic hybrid 22, which belongs to the R×Pr mapping population and was also derived from Flying Dragon.DNA extractions, sequencing, and bioinformatic sequence analysis were carried out by Sistemas Genomicos S.L. (Paterna, Valencia, Spain).Leaves (around 5 g) from the 107, 90, and 22-7 trees were used for DNA extraction with the aid of a DNeasy Plant Mini Kit (Qiagen, Barcelona, Spain).Three DNA libraries were generated using the NEBNext DNA Library Prep kit for sequencing on the Illumina HiSeq platform, with an expected output per sample of 40 Gb (2 × 150 pb, single index, 100× coverage).The quality control of initial reads was analyzed using the FastQC method.The reads were then mapped against the Citrus clementina genome (GCF_000493195.1).The low-quality mapping reads were filtered using the samtools method [69], and the PCR duplicates were eliminated with the aid of Picard tools.
The reads of candidate genes in the salt tolerance QTL were used to perform de novo assembly of alleles (a, b, and c/d) at each locus based on three steps: (1) mapping of initial reads to the whole genome using the bwa algorithm [70] and discarding low-quality reads and duplicated PCR products through post-processing alignment; (2) selection of reads mapped to the candidate genes in the QTL, sequence extraction to perform de novo assembly; and (3) haplotype identification using the WhatsHap algorithm [71] with default parameters.The sequences of candidate alleles were obtained using the haplotypes identified in the previous step and through multiple alignments with the aid of the Muscle program [72].The same two full sibs of the R×Pr population, hybrids 107 and 90, were selected for the root transcriptomic study at the end of the NG experiment.After 383 days of salt treatment, the whole root of five plants per genotype and treatment was taken and maintained with liquid nitrogen for further RNA extraction and to estimate Cl, Na, K, Ca, Fe, C, and N concentrations.
Total RNA from the whole root was purified as described elsewhere [73].Briefly, three independent biological samples, with one plant per sample and treatment, were used for analysis in an initial purification step with CTAB extraction buffer, followed by separation in chloroform:isoamyl alcohol (24:1, v:v) and LiCl precipitation [74].The resuspended RNA was further purified using the Aurum TM total RNA mini kit (Bio-Rad Laboratories S.A., Madrid, Spain), together with RNAse-free DNase in-column treatment (Promega Biotech Ibérica SL) according to the manufacturer's instructions.RNA quantity, quality, and integrity were determined based on absorbance ratios at 260 nm/280 and 260 nm/230 nm, using a mySPEC NanoDrop spectrophotometer (VWR International bvba, Leuven, Belgium) and agarose gel electrophoresis and further verified with the aid of a 2100 Bioanalyzer system (Agilent Technology, Madrid, Spain) at the López-Neyra Institute of Parasitology and Biomedicine (IPBLN-CSIC, Granada, Spain) genomic facility.Sample RIN values ranged from 5.9 to 8.2.

Library Preparation and Sequencing
The cDNA libraries were constructed and sequenced using pair-end sequencing (2 × 74 bp) with the aid of the Illumina NextSeq 500 System MidOutput 150-cycle at the IPBLN-CSIC genomic facility (Granada, Spain).A total of 298 million reads, ranging from 21.8 million in NGS083-20-4 to 26.1 million in NGS047-20-1, were generated and filtered for high-quality reads (Q30, 91.62%).Raw reads were deposited under the name Bioproject PRJEB61142.

Differential Expression, Correlation, and Clustering Analyses
Transcript counts were submitted to our pipeline RNAseqFlow [77] for (1) filtering (removal of transcripts with <1 count per million [CPM] reads in >9 samples); (2) normalization using the trimmed mean of M-values (TMM) method implemented in the edgeR v3.36.0 library [78]; (3) heteroskedasticity removal using the limma::voom() function [79] to facilitate reliable linear adjustment; (4) under a negative binomial distribution, fitting a generalized linear model with the limma v3.50.3 library; as in the present study, with its small number of replicates, its power advantage may outweigh the possibly exaggerated number of false positives [80]; (5) differential expression using the limma::treat() function [81] to enhance the biological implications of the results [79] and to prevent the misuse of statistical significance [82,83]; and (6) correlation-based clustering and further definition of communities using the igraph v1.3.5 library in order to detect gene co-expression networks [84][85][86].A differentially expressed gene (DEG) has at least an adjusted p-value < 0.1 and an absolute fold change |FC| > 1.2, as indicated by [81].Minimal Pearson correlation was set to 0.75 [87].Those genes being DEG in all contrasts or having a Kleinberg's centrality score > 0.9 [88] and a degree >10 were considered for further analyses.

Functional Analyses
Functional interpretations based on gene ontology (GO) and KEGG pathways were conducted using the closest Arabidopsis ortholog to the Citrus or Poncirus transcript.Orthologs were obtained using the Arabidopsis thaliana TAIR10 proteome from ENSEMBL and a reciprocal comparison based on the fast and sensitive protein aligner DIAMOND v.2.0.14 [89] using default parameters and -max-target-seqs 1 in -ultra-sensitive mode to obtain the best hit.The same protocol was used to obtain the orthology between the C. clementina v1.0 and Poncirus trifoliata v1.3.1 transcripts.REVIGO [90] was used for a list of GOs and KEGGs in which genes are involved.AgriGO v2.0 [68] and ShinyGO 0.76 [91] were used for functional enrichment with the aid of the background for A. thaliana, running on default parameters (in all cases, FDR < 0.05).

Conclusions
In this original approach to determine genes conferring salt tolerance in citrus, we describe the integration of QTL and RNA-seq analyses to identify transporter coding genes as likely involved in the leaf chloride exclusion mechanism of the salt-tolerant Cleopatra mandarin.Thus, we showed that the salt tolerance of Cleopatra mandarin rootstock is inherited by its R×Pr progeny, which is consistently controlled by at least one QTL on chromosome 6 (LCl-6) across experiments using both grafted and non-grafted apomictic progenies thus pointing to the root as the main organ involved.The study of differential gene expression at the root of two salt tolerance-contrasting full-sibs, which differ in relation to the Cleopatra allele in LCl-6, has identified some transporter coding genes in this QTL as positional candidate genes to explain xylem chloride exclusion under salinity conditions, as well as other genes that could be related to the growth and lignification of the root; this suggests that other salt tolerance mechanisms involved in and underlying salt tolerance in the same genomic region may be present.Markers developed for candidate genes in LCl-6 constitute new tools to increase the efficiency and speed of salt tolerance breeding programs for citrus rootstocks and to set up core collections of Citrus and Poncirus accessions where new salt-tolerant materials might be more easily discovered.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

2 Figure 1 .
Figure 1.Means and standard deviations of R×Pr hybrids and parents (Cleopatra, Cleo, and trifoli ata Orange, Rich) under both control and salinity conditions for (A) leaf [Cl − ], (B) root [Cl − ], and (C their difference relative to root [Cl − ].The salt-tolerant (107) and salt-sensitive (90) hybrids are indi cated by an arrow (black and red, respectively).

Figure 1 .
Figure 1.Means and standard deviations of R×Pr hybrids and parents (Cleopatra, Cleo, and trifoliata Orange, Rich) under both control and salinity conditions for (A) leaf [Cl − ], (B) root [Cl − ], and (C) their difference relative to root [Cl − ].The salt-tolerant (107) and salt-sensitive (90) hybrids are indicated by an arrow (black and red, respectively).

Figure 2 .
Figure 2. Co-expression network for differentially expressed genes (DEGs) in AHC-1 (FigureS13) in the salt tolerance quantitative trait locus (QTL) LCl-6.Names of transcripts and genes are provided, together with their putative function, showing that processes related to cell wall biosynthesis (specifically, glucuronoxylan and xylan biosynthesis, as well as lignin metabolism) and secondary metabolism (alkaloid/phytoalexin biosynthesis) at the root could be involved in salt tolerance.

Figure 2 .
Figure 2. Co-expression network for differentially expressed genes (DEGs) in AHC-1 (FigureS13) in the salt tolerance quantitative trait locus (QTL) LCl-6.Names of transcripts and genes are provided, together with their putative function, showing that processes related to cell wall biosynthesis (specifically, glucuronoxylan and xylan biosynthesis, as well as lignin metabolism) and secondary metabolism (alkaloid/phytoalexin biosynthesis) at the root could be involved in salt tolerance.

Figure 3 .
Figure3.Log of odd ratio (LOD) profiles of quantitative trait loci (QTLs) for Cl-related traits in linkage group 4c (A) and R4c (B) after incorporating candidate genes into the linkage maps and using Kruskal-Wallis statistic K for the association between phenotype and Cleopatra marker/gene segregation (C).The red horizontal line indicates significant levels.Cl*_S corresponds to the salt tolerance QTL reported by[17].To distinguish between experiments, the prefix GP or NG has been added to the trait name (Cl_L_S).

4. 6 .
RNA-Sequencing Analysis of the Root of Selected Materials 4.6.1.RNA Isolation of Selected Materials

Author
Contributions: M.J.A. and A.B. (Andres Belver) conceived and designed the research and also obtained funding; V.R., M.J.A., M.R.R.-A., J.E., J.A.T. and A.B. (Andres Belver) conducted the experiments.E.A.C. and M.J.A. analyzed the phenotypic data and interpreted the results; the molecular marker analyses were carried out by V.R. and M.J.A.; G.P.B. and J.C.T. obtained and analyzed the genomic sequences; A.B. (Amanda Bullones) and M.G.C. performed the bioinformatic analysis of RNA sequences; and M.J.A. wrote the manuscript.All authors have read and agreed to the published version of the manuscript.

Table 2 .
Means and standard errors for parents Cleopatra, trifoliate orange, whole population (Pop and the salt-sensitive (SS) and tolerant (ST), reference hybrids (90 and 107, respectively) depending on the experiment and treatment (Exp/E).Trait codes are summarized in Tables

Table 2 .
Means and standard errors for parents Cleopatra, trifoliate orange, whole population (Pop) and the salt-sensitive (SS) and tolerant (ST), reference hybrids (90 and 107, respectively) depending on the experiment and treatment (Exp/E).Trait codes are summarized in TablesS1 and S2for GP and NG experiments, respectively.

Table 3 .
List of positions in centiMorgan (cM), log of odd ratio (LODs), nearest markers, or marker intervals of quantitative trait loci (QTLs) detected by IM in the NG experiment in the integrated Citrus reshni-Poncirus trifoliata genetic linkage map (LG) using the cross-pollination model.The QTLs that were also detected in the individual parental linkage maps are indicated by putting the parental linkage group in parentheses (R and Pr for C. reshni and P. trifoliata, respectively).The four genotypic means (ac,

ad, bc, and bd, being
C. reshni ab and P. trifoliata cd), as well as the percentage of explained variance (PEV), are also indicated.LOD values higher than genome-wide significant LOD scores are indicated in bold.