Quantitative Trait Locus Mapping and Candidate Gene Identification for Fruit Acidity in Chinese Dwarf Cherry (Cerasus humilis) Using a High-Density Genetic Map

1. Introduction

The Chinese dwarf cherry (Cerasus humilis) is a rare and endemic fruit tree resource in China. It has recently gained attention due to its unique biological characteristics and economic value [1,2,3]. The fruit of C. humilis is rich in vitamins, amino acids, carbohydrates, organic acids, and other substances [4]. Cerasus humilis ranks among fruits with the highest calcium content, earning the title, "the star of calcium supplementation" [5]. It has high stress resistance (cold- and drought-resistant, barren-tolerant) and is valuable for ecological restoration (soil and water conservation) and fruit product development [6,7]. With the development of characteristic ecological agriculture, its cultivation area has gradually expanded. However, current cultivated varieties are highly acidic and primarily used for processing. Consumer-preferred, low-acid, fresh-eating varieties are scarce, impeding commercialization of C. humilis. Therefore, elucidating the genetic regulatory mechanism of fruit acidity and directionally improving fruit quality via molecular means are crucial for commercial promotion of C. humilis [8].

Genetic linkage maps are fundamental tools for indicating the genetic basis of fruit tree traits and also provide crucial support for quantitative trait locus (QTL) mapping [9,10]. The construction of maps relies on efficient molecular marker technologies [11]. Single nucleotide polymorphisms (SNPs), due to their large quantity, wide coverage, and high detection efficiency, have become the mainstream choice for genetic map construction [12]. Many people have constructed genetic linkage maps using this method [13,14,15]. Restriction site-associated DNA-sequencing (RAD-seq) is based on genome-wide restriction enzyme sites, with high-throughput sequencing conducted on specific restriction fragments [16]. This method considerably reduces genome complexity and sequencing costs, and rapidly identifies high-density SNP loci. It is widely used in genetic variation detection, genetic map construction, gene mapping for important traits, and population genetic evolution analysis [17,18]. This sequencing technique proves useful in examining the genetic mechanisms underlying target traits and has been widely applied to non-model organisms, with one or more genetic maps constructed in various fruit trees species [19].

Fruit acidity is a quantitative trait with complex genetic and variation mechanisms [20,21]. QTL mapping identifies chromosomal regions linked to target traits and estimates their genetic effects using association analysis of phenotypic data and marker genotypes. This provides a basis for marker-assisted selection, thereby improving breeding efficiency [22]. Diaz-Garcia et al. mapped 16 QTLs associated with titratable acidity in five cranberry fruit traits [23]. Tang et al. identified eight QTLs related to total acidity in jujube fruit using composite interval mapping, with LODs of 3.05–4.01 and a phenotypic variation explained (PVE) of 17.70–22.60% [24]. Jiang et al. used resequencing data to map 11 QTLs associated with fruit acidity in apricot and identified 88 genes related to malic acid, 310 related to citric acid, and 40 related to total acidity within QTL intervals [25].

Advances in molecular marker technologies have promoted the formulation of genetic linkage maps and have enabled effective QTL mapping [26,27]. Currently, it has been widely in various fruit trees, including pear [28], peach [29], plum [30], and apple [31]. However, studies on the genetic map and QTL mapping of C. humilis are limited, hindering progress in molecular marker-assisted breeding of C. humilis. Cerasus humilis has fewer chromosomes (2n = 2x = 16) and a small genome (~229.21 Mbp) [32]. Constructing genetic maps for QTL mapping is essential for identifying genes related to fruit traits in C. humilis.

In the current study, 208 F₁ hybrid individuals from the Nongda 4 × DS-1 cross were used as mapping materials. RAD-seq was used to develop SNP markers to construct a genetic map. QTL mapping for fruit titratable acidity (TA) and the screening of candidate genes were conducted. This study aimed to provide theoretical and technical support for molecular marker-assisted breeding and accelerate the development of breeding high-quality C. humilis varieties.

2. Materials and Methods

2.1. Plant Materials

All experimental materials, including two parents and F₁ hybrids, were collected from the C. humilis Shixiang Experimental Field of Shanxi Agricultural University (37°26′ N, 112°32′ E). The female parent was from the high-acid Nongda4 cultivar and the male parent from the low-acid DS-1 strain. Their hybrid progenies were planted in 2019. The experimental field had uniform water and fertilizer conditions, and field management was conducted referencing conventional field production management to ensure consistent growth environments for all plants.

2.2. RAD-Sequencing and SNP Statistics

Young leaves of all samples were collected in spring and treated with liquid nitrogen. Genomic DNA extraction was performed with a DNA extraction kit (Solarbio, Beijing, China). Library construction and high-throughput sequencing were performed by Lianchuan Biological Company using the Illumina HiSeq PE150 platform (Illumina, San Diego, USA). The quality-controlled sequencing data were aligned to the C. humilis reference genome (https://doi.org/10.6084/m9.figshare.11669673) using BWA software under default parameter settings, and the number of aligned reads was counted. Subsequently, variant calling was conducted using GATK and SAMtools [33]. SNPs consistently identified using both software were retained as reliable loci. Polymorphic SNPs were selected to form the final SNP dataset for subsequent genetic map construction.

2.3. Genetic Map Construction and Quality Evaluation

The double pseudo-test-cross method was used for the construction of a genetic map [34]. To ensure genetic map quality, polymorphic SNPs were filtered and genotyped [35]. SNPs with a parental sequencing depth <5× were excluded. Selected markers whose genotypes covered at least 75% of all offspring individuals were retained, with thresholds adjusted appropriately based on marker data quantity. Markers with a chi-square test (P < 0.001) were filtered out to avoid impacting QTL mapping. As a cross-pollinated population was used, only polymorphic markers other than the aa × bb type (homozygous) were retained to match population characteristics.

SNP markers were divided into eight linkage groups based on their positions in the reference genome. JoinMap was used to determine the linear order of markers within linkage groups and calculate genetic distances between neighboring markers to produce the genetic map [36]. Genetic map quality was assessed by analyzing mapped marker linkage relationships. Visualization of these relationships were performed by generating plots using an in-house R script.

2.4. QTL Mapping and Screening of Candidate Genes

Using the constructed genetic map of C. humilis and phenotypic data collected over two consecutive years [37], QTL mapping for fruit titratable acidity in C. humilis was conducted using the interval mapping method in MAPQTL software [25]. A threshold of LOD ≥3 was established to identify the presence of a QTL. QTLs with an LOD >3.5 were considered significant. The proportion of PVE for each QTL was computed. QTL mapping intervals were matched against the physical positions of the C. humilis reference genome to obtain detailed information on genes within these regions. All genes within QTL intervals were functionally annotated using databases (NR, Swiss-Prot, COG, GO, KEGG) to screen candidate genes related to titratable acid metabolism in C. humilis fruit.

3. Results

3.1. RAD Sequencing Data Analysis

RAD-seq was conducted on two parents and 208 F₁ individuals, generating 241.86 Gb of data. The maternal parent, Nongda 4, produced 0.79 Gb of data, with a base quality Q30 ratio of 91.15% and an average GC content of 41.61%. The paternal parent, DS-1, produced 1.78 Gb of data, with a base quality Q30 ratio of 90.87% and an average GC content of 41.7%. The 208 offspring produced 239.29 Gb of data, with an average base quality Q30 ratio of 90.65% and an average GC content of 40.82%. The average Q30 value of all samples was 90.89%, indicating that sample sequencing data were of good quality with high accuracy (i.e., low error rate) (Table S1). When aligned to the reference genome of C. humilis (Table S2), the maternal parent, Nongda 4, and paternal parent, DS-1, had an average sequencing depth of 5.89×, whereas that of the hybrid offspring was 5.25×; the alignment rates with the reference genome were 95.83%, 95.03%, and 93.71%, respectively.

3.2. Genetic Map Construction

After multiple rounds of rigorous filtering and genotyping, 2,491 SNP markers were finally confirmed to be usable. These included 1,537 "lm × ll," 668 "nn × np," and 286 "hk × hk" marker types. A genetic map was constructed that included 2,491 SNP markers, with eight linkage groups (LGs) covering a total genetic distance of 672.71 cM (Table 1; Figure 1). The mean interval between markers was 0.27 cM. Linkage group lengths ranged from 54.6 cM (LG8) to 142.71 cM (LG2), with an average length of 84.09 cM per group. The number of SNP markers ranged from 153 (LG8) to 465 (LG1). LG2 was the longest linkage group, containing 142.71 SNP markers. The shortest was LG8, with 153 SNP markers. LG2 had the longest average genetic distance (0.37 cM), and LG1 had the shortest (0.17 cM). The average proportion of gaps <5 cM in each linkage group was 99.49%.

3.3. Quality Evaluation of the Genetic Map

For assessing the quality of the genetic map, recombination plots were generated for each linkage group (Figure 2). A positive correlation was found between genetic distance and linkage strength across all linkage groups. This indicates that the constructed genetic map was of high quality.

3.4. QTL Mapping and Candidate Gene Identification

QTL mapping results for titratable acidity in C. humilis fruit showed that 18 QTLs were detected in the 2021 or 2022 datasets (Table 2; Figure 3), with 11 being significant. In 2021, 13 QTLs were detected, which were mapped to linkage groups LG2 (21TA-1 to 21TA-7) and LG5 (21TA-8 to 21TA-13), with LOD thresholds of 3.25–11.28 and PVE of 10.19–31.13%. In 2022, five QTLs were identified, fewer than in 2021. These QTLs were each mapped to one of the linkage groups (LG1, LG2, LG3, LG5, and LG6), with LOD thresholds of 3.14–3.86 and PVE of 7.1–8.7%.

Based on the results of QTL mapping and database annotation, two candidate genes were identified. Both genes were located on LG2: a phosphoester phosphatase gene (MSTRG.15867) and a MATE transmembrane transporter (MSTRG.17104). Transcriptome data from Nongda 4 and DS-1 [8] were integrated to analyze the expression levels of candidate genes identified using QTL mapping in fruit of two C. humilis germplasms at different developmental stages. Simultaneously, a comparison was made with fruit titratable acid content at the corresponding stages (Figure 4).

As C. humilis fruits developed, the relative expression levels of the phosphoester phosphatase gene of both Nongda 4 and DS-1 increased and subsequently decreased, aligning with the changes in titratable acid content. Moreover, the expression level in Nongda 4 was higher than that in DS-1. According to the correlation analysis (Table 3), the relative expression level of the phosphoester phosphatase gene was strongly correlated with titratable acid content (P <0.01), with a correlation coefficient of 0.93. This indicates that the gene is closely related to the regulation of fruit acidity in C. humilis.

4. Discussion

Cerasus humilis not only has strong stress resistance but its fruit is rich in nutrients and holds considerable value in both ecological and economic aspects. However, high fruit acidity in existing cultivars limits their suitability for fresh consumption. Therefore, elucidating the genetic mechanism underlying fruit acidity and performing targeted improvements are essential for progress in commercial utilization. In this study, the experimental material consisted of an F₁ hybrid population derived from crossing the high-acid Nongda 4 with the low-acid DS-1 varieties of C. humilis. A genetic linkage map was built utilizing RAD-seq technology. Two years of phenotypic data were used for QTL mapping of fruit titratable acidity and candidate gene identification. Genetic linkage maps are fundamental tools to indicate the genetic basis of fruit tree traits, and their density and accuracy directly affect QTL mapping efficiency [38].

In the current study, a genetic map was built using 208 F₁ hybrids from a cross between Nongda 4 and DS-1, offering advantages in parental selection and population size. The marked differences in acidity and other fruit traits between parents provided abundant materials for genetic variation, leading to extensive fruit trait segregation in the F₁ generation, ensuring sufficient variation for genetic map construction. The size of the F₁ segregating population critically affects the map’s accuracy, resolution, saturation, and application. A larger segregating population is conducive to constructing a more accurate genetic map [39]. Large-scale cultivation of fruit tree populations requires a large amount of land area. Subsequent genetic linkage mapping also demands substantial investment in human resources, materials, and finances. This results in relatively small population sizes used for constructing fruit tree genetic maps. Currently, most constructed fruit tree genetic maps are based on populations of <150 individuals [23,24], which may reduce mapping precision. However, the 208 F₁ individuals used in this study produced more molecular markers, substantially improving the map’s accuracy, resolution, saturation, overall representativeness, and reliability.

Traditional molecular marker technologies have single marker types and limited numbers, resulting in large gaps in constructed genetic maps, thus compromising map quality and subsequent research [40,41]. The advent of third-generation molecular SNP markers, genotyping by sequencing, and restriction enzyme digestion sequencing technologies has led to major advances in fruit tree genetic mapping [42]. RAD-seq, a reduced-representation genome sequencing approach based on genome-wide restriction enzyme digestion sites, reduces genome complexity and sequencing costs while enabling rapid identification of high-density SNP loci [43].

In the current study, RAD-seq technology generated 241.86 Gb of data, with an average Q30 value for all samples reaching 90.89%, indicating high sequencing quality data and low error rates, ensuring accurate SNP marker development. The constructed C. humilis genetic linkage map contained 2491 SNP markers. The eight linkage groups had a total genetic distance of 672.71 cM, with an average inter-marker distance of only 0.27 cM. Moreover, the average proportion of gaps <5 cM per linkage group was 99.49%. This map provides a solid platform for subsequent QTL mapping and gene discovery. The mapping results align with successful RAD-seq application in other fruit trees, including apple [42], pear [43], and hawthorn [44].

QTL mapping is a key method for elucidating the genetic mechanisms of quantitative traits; results are affected by multiple factors, with LOD score thresholds being particularly important [45]. The current study set an LOD score of 3.0 to detect QTL effects and 3.5 to define significant QTLs. These thresholds balance mapping accuracy with minimizing false-positives and the omission of minor-effect QTLs. Interval mapping identified 18 QTLs associated with fruit titratable acidity in C. humilis, of which, 10 were major-effect QTLs concentrated in LG2 and LG5. The distribution of QTLs showed similarities to and distinctions from fruit acidity-related QTLs in other fruit trees. For example, QTLs controlling acidity traits in peach have been mapped to LG5 [46], which concurs partly with the findings of this research. This suggests that the genetic regulation of fruit acidity may involve partially-conserved genetic mechanisms across different fruit tree species. Due to differences in genetic background, genome structure, and evolutionary history in different fruit tree species, QTL mapping results for fruit acidity also demonstrates obvious variation. For example, in studies on apricot [47], pear [28], and apple [48], the distribution and number of QTLs associated with fruit titratable acidity vary. The QTLs for C. humilis fruit acidity identified in the current study showed inconsistent intervals between 2021 and 2022, potentially due to environmental sensitivity to quantitative traits: differences in climatic conditions, pest and disease incidence, and harvest periods between 2021 and 2022 may have altered fruit acidity, affecting QTL mapping results. This also highlights the necessity of conducting multi-year and -location experiments in quantitative trait research. Yao [49] also found that quantitative traits were highly complex and environmentally-influenced. Therefore, stable QTLs represent a small proportion of all identified QTLs, with few being practically useful. Despite QTL instability across years, they still provide valuable insight into genetic regulatory mechanisms underlying C. humilis fruit acidity.

Three genes were screened within the QTL mapping intervals. After database annotation, two candidate genes on LG2 were identified: a phosphoester phosphatase gene (MSTRG.15867) and a MATE transmembrane transporter (MSTRG.17104). The relative phosphoester phosphatase gene expression level aligned with the variation in titratable acid content and was strongly correlated (0.93; P <0.01). This suggests that this gene is crucial for regulating fruit acidity in C. humilis. It is purported to influence fruit acidity by participating in organic acid metabolic pathways, such as synthesis regulation, decomposition, or transport of organic acids. Phosphatases may be involved in phosphorus metabolism-related signal transduction, influencing organic acid accumulation or degradation [50]. However, the specific regulatory mechanism remains unclear and requires further comprehensive verification using techniques such as gene editing and transgenesis.

MATE transmembrane transporters may regulate fruit acidity by mediating organic acid intracellular transport and influencing their distribution among different organelles or cells [51]. The gene identified in the current study was located within the QTL interval. Although the correlation between its expression level and fruit titratable acid content has not been verified, its potential role in regulating C. humilis fruit acidity remains plausible. Further studies, including gene expression analysis and protein function verification, are needed to elucidate its specific regulatory mechanism in C. humilis fruit acidity.

However, despite its contribution, this study had several limitations. Despite the high marker density of the genetic map, some linkage groups (e.g., LG4) had relatively large gaps, which may lead to the omission of some QTLs. Candidate gene functions were only inferred based on expression correlations without in vivo and in vitro validation. Non-coding RNAs or regulatory elements within QTL intervals were not analyzed, potentially ignoring epigenetic regulatory mechanisms.

Future research can be improved in three respects: applying higher-density sequencing and fine mapping to close gaps between and refine candidate gene regions; verifying candidate gene functions using gene editing or transgenic technologies to clarify their specific roles in acidity regulation; and integrating metabolome data to analyze the associations between candidate genes and organic acid metabolites (e.g., citric and malic acid) to elucidate related metabolic pathways.

5. Conclusions

A highly saturated and precise genetic linkage map was constructed for C. humilis, forming a strong basis for genomic research in the plant. Using this map and two years of phenotypic data, 18 QTLs associated with fruit titratable acidity were identified, including 11 major-effect loci concentrated on LG2 and LG5. These findings elucidate the genetic mechanisms underlying fruit acidity. A phosphoester phosphatase gene and a MATE transmembrane transporter gene were selected as candidate genes, with the former strongly correlated to acidity variation. The study provides theoretical and technical support for marker-assisted breeding and accelerates the development of high-quality, low-acid C. humilis cultivars.