SNP Loci and Favorable Haplotype Mining for Alkalinity Tolerance in the Rice Bud Stage

: The mining of favorable SNP loci and haplotypes is of great signi ﬁ cance for the further cloning and molecular-assisted breeding of alkalinity-tolerance genes in rice. To improve the utilization rate of saline–alkaline lands and the yield of rice, we used the 1,322,884 SNPs obtained from the sequencing 173 of rice accessions in this study. Alkalinity-tolerance-related traits, including the germination energy (GE), germination rate (GR), seedling length (SL), root length (RL), relative alkaline damage rate (RADR) of the GERADR of the GR, RADR of the SL and RADR of the RL, were evaluated in 2019 and 2020 and revealed abundant phenotypes in the studied population. A genome-wide association analysis detected 10 quantitative trait loci (QTLs) related to alkalinity tolerance. In addition, a transcriptome sequencing analysis of the alkalinity-tolerant rice variety ‘Yuedao 9’ and the alkali-sensitive rice variety ‘Tijin’ under alkalinity stress and control conditions was performed. Three candidate genes that were predicted to be related to alkalinity tolerance in rice, namely LOC_Os06g06600, LOC_Os011g44680 and LOC_Os011g44600, were screened based on gene annotation, coding sequences and haplotype analysis. The results of this study provide important genetic information for the molecular improvement of rice.


Introduction
Salinity-alkalinity stress has become the main factor restricting agricultural food production because it reduces the crop yield, affects the crop quality and limits the utilization rate of agricultural land [1].Rice (Oryza sativa L.) is one of the most important food crops and is widely distributed worldwide, and the monitoring of rice planting areas is necessary to ensure national food security [2].Rice is considered a salt-sensitive species, especially during the seedling and reproductive stages [3].Salt-alkalinity stress inhibits the growth and development of rice at various growth stages by reducing the nutrient solubility, increasing the external osmotic pressure and destroying the ion balance, which growth and led to a serious loss in crop yield [8][9][10].Therefore, studying the salinityalkalinity tolerance of rice not only has important practical significance for the breeding of salinity-alkalinity-tolerant rice accessions, but can also improve saline-alkaline land to ensure the development of sustainable agriculture.
Rice alkalinity tolerance is a complex trait controlled by multiple genes [11].To date, approximately 182 quantitative trait loci (QTLs) related to alkalinity tolerance have been identified [12][13][14][15][16][17].Li et al. [1] analyzed 295 japonica rice varieties and found that qSNK3-1 (15.0 Mb), which is located on chromosome 3, was correlated with ion homeostasis under salinity-alkalinity stress.Wu et al. [15] utilized the upland rice variety 'Xiaobaijingzi' as the maternal plant and the japonica rice variety 'Kongyu 131′ as the paternal plant for hybridization and self-breeding for multiple generations, which yielded 200 F9 population families to serve as the test materials.A total of 11 saline-alkali-related QTLs were detected on chromosomes 1, 2, 3, 5 and 12, and qARSH1-1 and qARRN1 were identified as possible new loci related to alkalinity tolerance at the rice bud stage.Wang et al. [16] constructed recombinant inbred lines (RILs) with 'Luhui 99′ and 'Shen Nong 265′ as the experimental materials.Under alkalinity stress, 15 QTLs related to the germination rate, germination potential, seedling length, root length, root number and dry matter production at the seedling stage were identified on chromosomes 1, 2, 5, 6, 7, 8 and 9. Li et al. [17] constructed RIL and japonica natural populations using the hybrid 'Kongyu 131′ (alkalinity-tolerant variety) and the upland rice variety 'Xiaobaijingzi' (alkalinity-tolerant variety).Through linkage mapping and genome-wide association analysis, these researchers found that qAT11 with F-box and DUF domains on chromosome 11 enhances the tolerance to alkalinity.
At least six genes associated with cloned alkalinity tolerance in rice have been cloned [18][19][20][21][22][23][24][25].OsLOL5 is a rice LSD1-like-type zinc finger protein (ZFP) gene containing two C2C2 domains, similar to LSD1.In rice with high OsLOL5 expression, oxidative-stressrelated genes are significantly upregulated, which improves the salt-alkalinity tolerance of rice and regulates oxidation and salt-alkalinity stress through reactive oxygen detoxification pathways [20].OsCSLD4 is a 1,4-β-D-xylan synthetase gene in the hemicellulose biosynthesis pathway [21].Loss of function of OsCSLD4 induces the constitutive activation of the defense response in rice [22].OsCSLD4 overexpression increases the grain width and weight.In addition, this overexpression enhances the expression of ABA synthesis genes, increases the ABA content and improves salt tolerance in rice, which suggests that OsCSLD4 regulates salt tolerance in rice through abscisic acid (ABA) synthesis [23].OsCSLD4-overexpressing strains exhibit a significantly higher survival rate than wild-type and mutant strains and thus show stronger salinity-alkalinity resistance [24].ALT1 encodes the chromatin-remodeling ATPase of the Snf2 family.Studies on the alkali-resistant mutant alt1 in rice revealed that this mutant exhibits low hydrogen peroxide levels and tolerance to methyl purpurine treatment under alkalinity stress, which can reduce the production of ROS and play a positive role in alkalinity resistance by defending against oxidative damage to improve the alkalinity tolerance of rice [25].OsSAP is a stress-associated protein that regulates soda saline-alkalinity tolerance through reactive oxygen species (ROS) homeostasis.OsSAP6 overexpression in rice significantly enhances salinity-alkalinity tolerance [26].OsPK3 affects pyruvate kinase activity in rice, and OsPK3 interacts with OsSAP6 to positively regulate the tolerance of rice in response to alkalinity stress [27].The OsSDG721 gene encodes H3K4 methyltransferase, and increased expression of the OsSDG721 gene can enhance the tolerance of rice in response to saline-alkalinity stress.An RNA-Seq analysis has demonstrated that OsSDG721 positively regulates high-affinity potassium (Na + ) transporter 1;5 (OsHKT1;5), maintains the K + /Na + dynamic balance under salt stress and is involved in the regulation of the saline-alkalinity stress response in rice [28].The plant height, grain size, grain weight and leaf angle are also affected by OsHKT1;5, which regulates the expression of genes related to rice growth [28].
To date, many QTLs and genes related to rice salt tolerance have been identified in previous studies, but the mapping of rice alkalinity-tolerance genes has been very limited.In fact, the damaging effect of alkali stress on rice is more serious and complex than that of salt stress [29].Rice has different response mechanisms to saline-alkalinity stresses [30].To further study the mechanism of alkalinity resistance, it is necessary to mine more genes related to alkalinity resistance.Previous research on alkalinity resistance has mainly focused on the seedling stage, and few studies have investigated the bud stage.The present study aimed to fill this research gap.To that end, the main goals of this study were to identify the loci that are significantly related to alkalinity stress at the rice bud stage and to predict and screen candidate genes related to alkalinity resistance.This study utilized a genome-wide association study (GWAS), transcriptome sequencing and bioinformatics to obtain results efficiently and rapidly [31].The results provide a reference for research on the resistance of rice to alkalinity and a foundation for further cultivation of salinity-alkalinity-resistant rice varieties.

Plant Materials and Field Planting
The population consisted of 173 rice accessions with a wide variation in alkalinity tolerance from multiple regions and countries in Asia (Table S1, Figure S1); these accessions were collected, preserved and provided by the State Key Laboratory of Crop Genetics and Germplasm Innovation of Nanjing Agricultural University.All 173 accessions were planted in the regular growing season (May to October) at the Jiangpu Experimental Station of Nanjing Agricultural University in Nanjing in 2019 and the Lujiang Experimental Station of Anhui Agricultural University in Hefei in 2020.A randomized block design with two replicates was used for both two-year experiments.Each plot contained 40 plants planted in five rows with a hill spacing of 17 cm × 20 cm.The plots were managed based on conventional agricultural management measures.The main stem panicles of the 10 plants in the middle three rows of each plot in the randomized block design were harvested, threshed and dried under natural sunshine.

Identification of Alkalinity Resistance
The seeds of 173 rice accessions collected in 2019 and 2020 were used to identify their traits of alkalinity resistance.Thirty clean seeds of 173 rice accessions were placed in Petri dishes with pure water as the control and 15 mmol/L Na2CO3 as the stress treatment concentration.Three replicates were established for both the control and treatment groups, and 10 mL of pure water and Na2CO3 was added.On the fourth day, the number of germinated seeds was investigated, and the relative alkaline damage rate (RADR) of the germination energy (GE) was calculated using the following formula.On the tenth day, the seed germination number was investigated, and the RADR of the germination rate (GR) was calculated.At this time point, the seedling length (SL) and root length (RL) of the materials were measured, and the RADR of the RL and RADR of the SL were calculated using the following formulas.RADR of the SL (%) = {(control seedling length − treatment seedling length)/control seedling length} × 100%

SNP Genotyping
For each of the 173 accessions to be sequenced, two blades from a single plant were collected at the tillering stage (1 month after seedling transplantation) for the extraction of genomic DNA using the standard cetyltrimethylammonium bromide protocol.Sequencing was performed based on a single genotype.Library construction, sequencing and sequence cleaning were performed by staff at Mega Genomics Beijing (http://www.genomics.cn/mobile.php/,accessed on 10 April 2019).A total of 1,322,884 SNPs were identified.The sequencing data presented in the study were deposited in the NCBI Sequence Read Archive repository with the accession number PRJNA554986.

Genome-Wide Association Analysis
Our previous structural analysis results showed that this population can be divided into two subgroups [32].In the present study, SNPs and phenotypic data were used for a GWAS based on a general linear model (GLM) and a mixed linear model (MLM) [33] in the R package "RMVP".A GLM can be used to test the genetic marker (S) [34], but when performing the analysis, the site effect value is overestimated such that false-positive results are obtained.An MLM can effectively correct data related to the population structure and the complex genetic relationships within populations.The p values were generally corrected by the Bonferroni correction method.We adjusted the threshold values according to the results from the actual Manhattan map and selected p = 10 −4 as the threshold for a significant association between SNP markers and traits.The SNPs in the same linkage disequilibrium (LD) region were regarded as QTLs, and the smallest SNPs were regarded as lead SNPs.The LD analysis of the population showed that the LD attenuation distance of the whole genome was approximately 112 kb [32].

Transcriptome Analysis
In this study, transcriptome sequencing was performed using seedlings and roots of the alkalinity-tolerant 'Yuedao 9' rice accession and the alkalinity-sensitive 'Tiejin' rice accession after 3 days of growth under control and stress conditions.Three replicates of each treatment were sampled, resulting in a total of 12 sequencing samples.Immediately after sampling, the samples were frozen in liquid nitrogen.
Part of the experiment was performed by Gidio Biotechnology Co., Ltd.(San Diego, CA, USA) Total RNA was extracted using a TRIzol kit.An Agilent2100 bioanalyzer was used to assess the RNA quality, and agarose gel electrophoresis was used to detect the contamination and degradation of RNA.After total RNA extraction, eukaryotic cell mRNA was enriched with oligo (dT) magnetic beads.The mRNA was then cut into short fragments using fragment buffer and reverse-transcribed into cDNA using the NEBNext Ultra RNA Library Kit (NEB, Ipswich, MA, USA).The purified double-stranded gene fragments were end-repaired, and a base was added and connected to Illumina sequencing connectors.AMPureXP (1.0×) was used to purify the ligation reaction.Agarose gel electrophoresis and PCR were utilized for amplification, and sequencing was then conducted with an Illumina NovaSeq 6000.
The 'Yuedao 9 (YD)' samples under the control and stress conditions were labeled YD and YDT, respectively, and the 'Tijin (TJ)' samples under the control and stress conditions were labeled TJ and TJT, respectively.For these four samples, we constructed four comparison schemes: YD vs. YDT, TJ vs. TJT, YD vs. TJ, and YDT vs. TJT.Fragments per kilobase million (FPKM) values were calculated to compare the gene expression differences between the samples, and DESeq2 (R4.1.2) software was used to analyze expression differences between any two groups.Differentially expressed genes were identified based on the following criteria: false discovery rate (FDR) < 0.05 and | log2Fold Change | > 1. Afterward, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed.

Candidate Gene Prediction and Haplotype Analysis
Based on the analysis of genome-wide LD decay, we estimated the candidate regions of chromosomes in the candidate area using the Rice Genome Annotation Project database (http://rice.uga.edu/index.shtml,accessed on 19 October 2023) to identify all the genes in this region.Combined with the differentially expressed genes identified from the transcriptome sequencing results, possible candidate genes were screened.Utilizing the Nipponbare reference genome sequence for comparison, we focused on the nonsynonymous SNPs that could cause amino acid changes in the GWAS results.Functional annotation, gene annotation and GO enrichment analysis of the candidate differentially expressed genes identified from the different combinations were performed.These genes were then selected as candidate genes related to alkalinity resistance in rice.The SNP loci on exons of candidate genes were extracted from the sequencing data for haplotype analysis, and the differences in phenotypic values between candidate gene haplotypes were compared by Tukey's test.

Identification of Alkalinity-Tolerant Phenotypes of Rice at the Bud Stage
The frequency data for alkalinity-tolerance-related traits were normally distributed in both years, which revealed the genetic characteristics of quantitative traits controlled by major genes or multiple genes (Figure 1).The mean, standard deviation (SD), range, coefficient of variation (CV) and generalized heritability (HB 2 ) for the GE, GR, RADR of the GE, RADR of the GR, SL, RL, RADR of the SL and RADR of the RL were derived for the natural populations.The r analysis of variance showed that the phenotypic data of alkalinity tolerance between the two blocks did not significantly differ from the randomized block design experiment in the field, whereas significant differences were found among accessions.Among the 173 accessions, the RADR of the GE ranged from 0 to 6.74% (2019) and from 0 to 7.95% (2020) with CVs of 61.35% (2019) and 67.92% (2020).The RADR of the GP ranged from 0 to 5.56% (2019) and from 0 to 6.67% (2020) with CVs of 81.83% (2019) and 80.72% (2020).The ranges of the RADR of the SL in 2019 and 2020 were 2.98-58.12%and 1-70.84%, with CVs of 56.28% and 50.86%, respectively.The RADR of the RL in 2019 and 2020 ranged from 29.91 to 90.8% and from 30.65 to 90.96%, with CVs of 81.20% and 89.68%, respectively (Table 1).The above-described results indicated that the effects of alkalinity stress on the SL, RL, RADR of SL and RADR of the RL were greater than those on the GE, GR, RADR of the GE and RADR of the GR.The response of the rice plants to alkalinity stress was mainly reflected in roots and seedlings.
We selected 12 different alkali-tolerant varieties at the bud stage from 173 accessions to show the wide variation in the population (Figure 2), and the actual variation in the population was mainly assessed using the CV (Table 1).The results showed that the 173 rice accessions exhibited wide variation in all traits associated with alkalinity tolerance at the rice bud stage and are thus suitable for association mapping.A correlation analysis showed a significant positive correlation between the GE and the GR.The correlation between the RADR of the GR and the RADR of the GE was significantly positive.A significant negative correlation was detected between the GE and the RADR of the GE.The correlation between the GR and the RADR of the GR was significantly negative (Figure 3).

Genome-Wide Association Analysis
The genome-wide association results are shown in Figures 4 and 5 and Table 2.  2 shows the GWAS results that overlapped between the two years and the two models.We found 23 SNPs associated with alkalinity resistance in rice, and among these, one locus was associated with the GE and located on chromosome 11.Two loci located on chromosomes 3 and 4 were associated with the GR.One locus was associated with SL and located on chromosome 5.A total of 15 SNP loci on chromosomes 1, 4, 6, 8, 9, 10 and 12 were associated with the RL.One locus on chromosome 2 was associated with the RADR of the GE.One locus on chromosome 5 was associated with the RADR of the GR, and two loci on chromosomes 1 and 12 were associated with the RADR of the RL.The LD decay distance was used to determine the candidate region.The QTLs related to alkalinity stress detected in this study were compared with salinity-alkalinity tolerance QTLs identified in previous studies.All nine QTLs discovered in this study, namely qGR4, qSL5, qRL1-1, qRL4-1, qRL4-2, qRL6, qRADRGR5, qRADRRL1 and qRADRRL12, are newly discovered QTLs (Figure 6).
The average GC content (ratio of guanine to cytosine) in the sequencing results was approximately 51.32%.Q20 and Q30 represent the percentage of bases with mass values greater than or equal to 20 and 30 in high-quality sequences.In this study, the average Q20 and Q30 percentages reached 97.34% and 92.85%, respectively.
As shown in Figure S4, the correlation between biological replicates of all tested samples was high, and this high correlation indicated good biological repeatability.The correlation coefficient between different treatments for the same accession was higher than 0.863, indicating a strong correlation.Combined with the above-described data, these results demonstrated that the quality of the transcriptome sequencing was excellent and that these data were suitable for further bioinformatics analysis.
Based on the FPKM values of each gene, we visualized the expression distribution of different sample genes or transcripts (Figure S5).Gene expression profiles were used to assess the differences in database construction, sequencing, comparison or quantification between some samples and other samples.Figure S5a reflects the concentration intervals of gene expression in different samples.Figure S5b shows the degree of dispersion in all gene expression levels, which more directly demonstrates the gene expression levels among different samples.
We used the following criteria to identify significantly differentially expressed genes: | log2Fold Change | > 1 and FDR < 0.05.As shown in Figure 7, under non-alkalinity-stress conditions, 4898 differentially expressed genes (3567 upregulated and 1331 downregulated genes) were found between YD and TJ.Under alkalinity stress, 234 differentially expressed genes (121 upregulated and 113 downregulated genes) were identified between YD and YDT.Moreover, 1753 differentially expressed genes (579 upregulated and 1174 downregulated genes) were found between TJ and TJT, and 5638 differentially expressed genes (3347 upregulated and 2291 downregulated genes) were detected between YDT and TJT.As shown in Figure 8, 96 differentially expressed genes were shared between the YD vs. YDT and TJ vs. TJT comparisons, 2998 differentially expressed genes were shared between the YD vs. TJ and YDT vs. TJT comparisons, and 34 differentially expressed genes were shared between the four comparisons.Based on the GO enrichment analysis results, we focused on the DEGs that may affect the saline-alkalinity tolerance of rice.In the YD vs. YDT comparison, 11 genes were potentially related to the salt-alkalinity tolerance of rice, and these were significantly enriched in three GO terms (Figure 9a).In the TJ vs. TJT comparison, 71 genes related to the salt-alkalinity tolerance of rice were significantly enriched in eight GO terms (Figure 9b).To further identify genes that affect rice growth and development under alkalinity stress, a KEGG metabolic pathway analysis of the differentially expressed genes was performed.From the combinations YD vs. YDT, YD vs. TJ, YDT vs. TJT, and TJ vs. TJT, the first 20 pathways with the lowest q value were selected for mapping, in which the ordinate and abscissa show the pathways and enrichment factors.The number of differentially expressed genes in a pathway was divided by all the numbers in the pathway, and the results represented the quantity; a redder color reflects a lower "q value" (Figure 10).

Analysis of Alkalinity-Tolerant Candidate Genes
A total of 221 genes with nonsynonymous mutations were identified (Table S2).A comparison of the differentially expressed genes identified from these four comparisons with the nonsynonymous mutant genes revealed that 30 genes overlapped between the GWAS and transcriptome analysis (Table S4).These DEGs were involved in a wide spectrum of biological processes and metabolic pathways (Figures 9 and 10).
Functional genes related to kinase and calmodulin were detected by GO enrichment analysis of these 30 genes (Table S3).Afterward, the candidate genes were further screened based on the Rice Genome Annotation Project database.We mainly focused on the DEGs that are potentially involved in the alkalinity-tolerance pathway.Three candidate genes, LOC_Os06g06600, LOC_Os11g44680 and LOC_Os11g44600, which encode the F-box domain and calmodulin-binding protein, were identified.These three genes are predicted to be likely involved in the regulation of the alkalinity tolerance of rice during the bud stage.

Haplotype Analysis of Candidate Genes
Among the three candidate genes, the LOC_Os06g06600 gene was located at the qGL6 locus and had three haplotypes, as revealed from the SNP analysis of cDNA from all 173 accessions (Figure 11).A statistical analysis of the relative alkaline damage rates of four traits in 2019 and 2020 showed no significant difference among the haplotypes of three traits (RADR of the GE, RADR of the GR and RADR of the SL).The average values of the three haplotypes of the RADR of the GE ranged from 2.2 ± 1.57% to 2.8 ± 1.64%.For the RADR of the GR, the average values ranged from 1.61 ± 1.59% to 1.82 ± 1.33%, and those of the RADR of the SL were ranged from 25.83 ± 13.43% to 32.17 ± 11.40%.Significant differences were found among the three haplotypes of the RADR of the RL.The average values of the three haplotypes of the RADR of the RL ranged from 65.69 ± 7.75% to 72.36 ± 7.73%.The value for HapC was significantly lower than that obtained for HapA and HapB, which indicates that HapC is the favorable haplotype of LOC_Os06g06600 for alkali tolerance.Therefore, LOC_Os06g06600 may be a candidate gene, and further validation through subsequent experiments is needed.The analysis of the LOC_Os11g44680 gene showed that all 173 rice accessions can be divided into seven haplotypes based on the SNPs in the cDNA (Figure 12).Statistical analysis of the relative alkaline damage rates of four traits in 2019 and 2020 showed no significant difference among the seven haplotypes of the RADR of the GR, and the average value of the seven haplotypes of the RADR of the GR ranged from 0.28 ± 0.40% to 2.04 ± 1.35%.Significant differences were detected among the seven haplotypes of the RADR of the GE, RADR of the SL and RADR of the RL.The average value of the seven haplotypes of the RADR of the GE ranged from 0.26 ± 0.40% to 3.09 ± 1.49%, and the value of HapG was significantly lower than that of the other six haplotypes.The average of the seven haplotypes of the RADR of the SL ranged from 7.55 ± 4.13% to 30.82 ± 8.11%, and that of HapG of the RADR of the SL was significantly lower than that of the other six haplotypes.The seven haplotypes of the RADR of the RL ranged from 52.11 ± 0.57% to 74.51 ± 4.67%, and the value of HapG was significantly lower than that of the other six haplotypes.Therefore, HapG is the favorable haplotype of the LOC_Os11g44680 gene.The mechanism of alkalinity tolerance based on HapG of LOC_Os11g44680 will be investigated in future studies.We analyzed the SNPs of cDNA exons at the qGE11 locus and found that LOC_Os11g44600 included nine nonsynonymous SNPs (Figure 13).The statistical analysis of the relative alkaline damage rates for the four traits in 2019 and 2020 revealed no significant difference among the nine haplotypes of the RADR of the GR, and the average value of the nine haplotypes of the RADR of the GR ranged from 0.28 ± 0.40% to 2.21 ± 1.34%.However, significant differences were observed among the nine haplotypes of the RADR of the GE, RADR of the SL and RADR of the RL.The average values of the nine haplotypes of the RADR of the GE ranged from 0.29 ± 0.40% to 3.32 ± 1.48%; for the RADR of the SL, the average values of the haplotypes ranged from 7.55 ± 4.13% to 33.79 ± 11.15%; and for RADR of the RL, the average values of the nine haplotypes ranged from 52.11 ± 0.57% to 75.66 ± 3.62%.For the RADR of the GE, the value of HapG was significantly lower than that of the other eight haplotypes.For the RADR of the SL, HapG had the lowest average RADR of the SL among all other haplotypes.HapG was an elite haplotype for the RADR of the RL.These results indicated that LOC_Os11g44600 might be a candidate gene for alkalinity tolerance in rice, and this finding needs to be verified by transgenic testing in the future.

Discussion
Although only 173 rice accessions were used in this study, these are all core germplasms.The coefficient of variation analysis results indicate that the population exhibits a wide range of variation in alkalinity-tolerance-related phenotypes during the budding stage and is suitable for association analysis, which is similar to the results reported in previous studies.Li et al. [1] used 295 japonica rice varieties and found that the CV of phenotypic data ranged from 27.82% to 45.50% under alkaline treatment at the seedling stage.Zhao et al. [35] used a Na2CO3 + NaHCO3 alkaline mixture at a concentration of 0.20% for the treatment of germinating seeds of 173 natural japonica rice populations and found that the CV of each trait ranged from 2.6% to 40.6%.Wu et al. [15] used 40 mmol•L −1 NaHCO3 for the treatment of the RIL population of "little white japonica/Kongyu131" (200 strains) in a germination experiment and found that the CV of phenotypic data ranged from 8.55% to 86.21%.Xing et al. [36] used 180 japonica rice varieties and measured the CV of each trait, which ranged from 0.28% to 58.23% under alkalinity stress (0.15% Na2CO3) at the seedling stage.Although the number of rice accessions included in the natural population of the 173 rice accessions in this study was not large, the natural population exhibited a wide variation in the alkalinity-resistance phenotype at the bud stage in contrast to the natural populations used in previous studies.
In this study, 14 QTLs were found to be similar to those previously reported.Zou et al. [37] used 180 families of the F2:3 generation to identify qRN1-2 between RM529 and RM302 on chromosome 1, which was found to be related to the root number through QTL detection of the root number (RN), root length (RL), chlorophyll content and relative alkaline damage rate (RADR) at the early seedling stage, and qRL1-2 in this study was in the qRN1-2 region.Lin et al. [38] used the F2 and F3 generations obtained by crossing Nona Bokra, an indica rice variety with strong salt tolerance, with Koshihikari, a disease-susceptible japonica rice variety, as genetic materials and identified three QTLs for seedling survival days (SDSs) under salt stress on chromosomes 1, 6 and 7 and one QTL related to Na + concentration (RNC) in roots under salt stress on chromosome 9.In this study, qRL1-3 on chromosome 1 was in the qSDS-1 (between C813 and C86) region.qRL9 overlaps with qRNC-9 (between R1751 and R2638) on chromosome 9. Li et al. [39] used the F2:3 population obtained after the hybridization of "Caidao" and "WD20342" as genetic materials, 151 SSR markers were used to demonstrate that qSDS1 on chromosome 1 was related to SDS), and the location of qRL8-6 identified in this study was similar to that of qSDS1.Liang et al. [40] used 200 recombinant inbred lines (F8) and detected seven additive QTLs and three pairs of AA epistatic QTLs related to the dead leaf rate (DLR) and dead seedling rate (DSR) under salt-alkalinity stress.These researchers also constructed genetic linkage maps of 155 SSR markers covering the 1541.5-cMwhole rice genome and detected qDSRs8-1 with a phenotypic variation of 15.96%.qRL8-1, qRL8-2, qRL8-3 and qRL8-4 were in the qDSRs8-1 region in this study.Singh et al. [41] used recombinant inbred lines (RILs) hybridized from "Cocodrie" and "Dular" and found that qRL10 was located at qSNC10.16 (aboveground Na + concentration), qGE11 was located at qSNK11.27 (aboveground Na+:K+ ratio) and qRADRGE2 was located at qSKC2.19 (aboveground K + concentration).Qi Dongling [42] used 200 F2:3 lines of "106/Changbai 9" as the mapping population to evaluate the GR and RADR of rice under alkalinity stress.The results also showed that the identified qGC-11-1 related to the GR overlapped with the qGE11 region in this study.Singh et al. [43] used 193 RILs hybridized from "Cocodrie" and "N22" to identify populations at the seedling stage under alkalinity stress.These researchers found that qRSR8.17 was related to the root-shoot ratio and that qSNC12.19 was related to the aboveground Na + concentration.In this study, qRL8-5 was found to be close to qRSR8.17, and qRL12 was close to qSNC12.19.Sabouri et al. [44] used "Tarommahali" and the salt-sensitive variety "KZAAR" to construct a linkage map based on the F2 population.qDM-3 and qDM-8 detected at the seedling stage contributed 20.90% and 17.72% to the total phenotypic variation, respectively.qGR3 is located in the qDM-3 region.All other nine QTLs found in this study (qGR4, qSL5, qRL1-1, qRL4-1, qRL4-2, qRL6, qRADRGR5, qRADRRL1 and qRADRRL12) were newly discovered QTLs (Figure 6) .In this study, 23 colocalized QTLs were selected as candidate genes.The GWAS results were compared with the transcriptome analysis results, and three alkalinity-resistance candidate genes were identified: LOC_Os06g06600, LOC_Os11g44680 and LOC_Os11g44600.LOC_Os06g06600 belongs to the F-box protein family.F-box proteins constitute a large family in eukaryotes with conserved F-box motifs consisting of approximately 40 amino acids and are critical for the controlled degradation of cellular proteins [45].F-box domain proteins, including OsMsr9, MAIF1, OsPP12-A13 and OsFBT4, participate in the regulation of biological stress responses.The F-box gene OsMsr9 in rice is induced under different stresses and improves plant salt tolerance [46].Under salt stress, transgenic rice plants overexpressing OsMsr9 exhibit an increased root length, a greater plant height and improved survival rates compared with wild-type plants.MAIF1 encodes a protein containing the F-box domain, controls root growth and plays a negative regulatory role in the response to abiotic stress [47].ABA and abiotic stress can induce MAIF1 expression.Plants overexpressing MAIF1 exhibit decreased sensitivity to ABA and reduced resistance to drought, salt and low temperature [48].OsPP12-A13 is a nucleus-localized protein that is strongly upregulated under salinity stress in rice.Studies have shown that OsPP12-A13 may improve the salt tolerance of rice by affecting Na + transport and the ROS dynamic balance [49].OsFBT4 belongs to a small subclass of rice F-box proteins that contain the conserved N-terminal Fbox domain.In rice, OsFBT4 transcripts are strongly upregulated under different abiotic stresses, including exogenous ABA.Plants overexpressing OsFBT4 show morphological differences compared with wild-type plants, which suggests that OSFBT4 plays a role in plant development [50].LOC_Os11g44680 and LOC_Os11g44600 encode calmodulinbinding proteins that can regulate cell function by binding to calcium.OsMSR2 is a new calmodulin-like gene in rice, and transgenic plants overexpressing OsMSR2 exhibit greater salt tolerance [51].Xu et al. [52] analyzed the recombinant OsMSR2 protein to demonstrate its potential ability to bind Ca 2+ in vitro.The regulation of OsMSR2 through the ABA-mediated pathway results in greater salt tolerance in Arabidopsis.OsCam1-1 is a calmodulin isolated from the rice genome by PCR amplification [53].Compared with wild-type and control plants, transgenic rice overexpressing the OsCam1-1 gene shows an upregulated expression of two genes involved in ABA biosynthesis and a higher ABA content [54].OsCam1-1 signaling may play an important role in ABA biosynthesis and affect the salt tolerance of rice [55].Calmodulin plays not only a role in the hypersaline response but also a regulatory role in drought resistance in rice.Transgenic plants overexpressing OsCML4 showed better growth and greater survival than wild-type plants under drought conditions, and these transgenic plants show significantly increased activities of ROS, superoxide dismutase (SOD) and catalase (CAT) and a significantly higher concentration of proline.The expression levels of the ROS clearance genes APXI and Cat-B and the stress-related gene OsP5CS1 are also enhanced under drought conditions.These results demonstrate that OsCML4 enhances drought resistance by partially clearing ROS and inducing other stress-related genes in an ABA-independent manner [56].
Through a haplotype analysis of the three candidate genes in this study, we found that HapC of the RADR of the RL in the LOC_Os06g06600 gene was significantly superior to HapA and HapB.These materials had a lower RADR of the RL and showed stronger resistance to alkalinity stress.Therefore, HapC is an excellent haplotype of LOC_Os06g06600, and HapC contains 'Nannongjing 1R', 'Longdao 6′, 'Longdao 8′, 'Zhendao 10′ and 'Yuedao 3′.Rice accessions containing the elite haplotype HapC are elite parents for improving the alkalinity resistance of rice.HapG of the RADRs of the GE, SL and RL in the LOC_Os11g44680 gene was significantly superior to the other six haplotypes.The RADRs of the GE, SL and RL of these materials were lower, showing strong alkalinity stress resistance.Therefore, HapG is an elite haplotype of LOC_Os06g06600, and HapG is found in dry grain accessions in the Nuo and Lincang Wa regions.Moreover, HapG is found in waxy dry valley and dry valley materials in the Lincang Wa ethnic region.Rice accessions containing the elite haplotype HapG are elite parents for cultivating alkalinity-resistant rice.LOC_Os11g44600 could divide 173 rice accessions into nine haplotypes.According to the phenotypic data, HapH of RADR of GE, RADR of SL and RADR of RL was significantly better than that of the other eight haplotypes.The RADRs of the GE, SL and RL were low, showing strong alkali stress resistance.Therefore, HapH is an elite haplotype of LOC_Os06g06600.HapH was observed with dry grain materials in the Nuo and Lincang Wa regions, and these accessions can be used as elite parents for cultivating alkalinity-resistant rice.

Conclusions
A total of 23 significant SNP loci were detected on chromosomes 6 and 11 by GWAS analysis during 2019 and 2020, and these include 9 newly discovered QTLs.A comparison of the results of the genome-wide association analysis and transcriptome sequencing identified three candidate genes, LOC_Os06g06600, LOC_Os011g44680 and LOC_Os11g44600, and the favorable haplotypes of these genes were identified.These results provide a good molecular basis and insight into germplasm resources for the breeding of alkalinity-tolerant rice varieties.

Supplementary Materials:
The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/agronomy13122954/s1, Figure S1.Geographic distribution map of 173 rice accessions used in this study.S1.Names and origins of 173 rice accessions used for association mapping.Table S2.Genes with nonsynonymous single nucleotide variation in GWAS.Table S3.GO enrichment analysis of genes in four combinations.Table S4.Functional annotation of candidate genes.

Figure 1 .
Figure 1.Histogram of the phenotypic frequency distribution in 173 rice accessions.(a) Germination energy; (b) germination rate; (c) seedling length; (d) root length; (e) relative alkaline damage rate of the germination energy; (f) relative alkaline damage rate of the germination rate; (g) relative alkaline damage rate of the seeding length; (h) relative alkaline damage rate of the root length.

Figure 2 .
Figure 2. RL and SL phenotypes of 12 rice accessions at the germination stage under alkalinity treatment.Bar = 2 cm.

Figure 3 .
Figure 3. Correlation analysis of alkalinity-tolerance-related traits in 173 rice accessions.** indicates significance at the 0.01 level; * indicates significance at the 0.05 level.
Figures 4, 5, S2 and S3 show the Manhattan plots of the GWAS results obtained with the GLM and MLM.Table

Figure 4 .
Figure 4. Manhattan plots for alkalinity-tolerance-related traits.(a) Manhattan plot for the GE in 2019; (b) Manhattan plot for the GR in 2019; (c) Manhattan plot for the SL in 2019; (d) Manhattan plot for the RL in 2019; (e) Manhattan plot for the GE in 2020; (f) Manhattan plot for the GR in 2020; (g) Manhattan plot for the SL in 2020; (h) Manhattan plot for the RL in 2020.The red arrows represent the QTLs identified in this study.

Figure 5 .
Figure 5. Manhattan plots for alkalinity-tolerance traits based on relative alkaline damage rates.(a) Manhattan plot for the RADR of the GE in 2019; (b) Manhattan plot for the RADR of the GR in 2019; (c) Manhattan plot for the RADR of the SL in 2019; (d) Manhattan plot for the RADR of the RL in 2019; (e) Manhattan plot for the RADR of the GE in 2020; (f) Manhattan plot for the RADR of the GR in 2020; (g) Manhattan plot for the RADR of the SL in 2020; (h) Manhattan plot for the RADR of the RL in 2020.The red arrows represent the QTLs identified in this study.

Figure 6 .
Figure 6.QTLs detected in natural populations in 2019 and 2020.The red and blue colors represent genes and QTLs reported in previous studies, respectively.

Figure 7 .
Figure 7. Volcano plots of gene expression in the four comparisons: (a) volcano plot of the YD vs. YDT comparison; (b) volcano plot of the YD vs. TJ comparison; (c) volcano plot of the YDT vs. TJT comparison; (d) volcano plot of the TJ vs. TJ comparison.

Figure 8 .
Figure 8. Distribution of differentially expressed genes among the four comparisons.

Figure 9 .
Figure 9. GO annotation of DEGs identified from the YD vs. YDT and TJ vs. TJT comparisons.(a) GO enrichment analysis of DEGs related to the saline-alkalinity tolerance of rice identified from the YD vs. YDT comparison; (b) GO enrichment analysis of DEGs related to saline-alkalinity tolerance

Figure 10 .
Figure 10.KEGG enrichment bubble diagram of differentially expressed genes identified from the four comparisons: (a) bubble diagram of the YD vs. YDT comparison; (b) bubble diagram of the YD vs. TJ comparison; (c) bubble diagram of the YDT vs. TJT comparison; (d) bubble diagram of the TJ vs. TJT comparison.

Figure 11 .
Figure 11.Haplotype analysis of the candidate gene LOC_Os06g06600.(a) Local Manhattan plot (top) and linkage disequilibrium heatmap (bottom).The candidate region lies between the red solid lines.(b) Schematic representation of the LOC_Os06g06600 structure and single-nucleotide polymorphisms in LOC_Os06g06600 cDNA between the two alleles.The blue boxes indicate exons.(c) Box plots for the RADR of the GE, RADR of the GR, RADR of the SL and RADR of the RL in the two alleles of LOC_Os06g06600 in all accessions in 2019 and 2020.The number of accessions (n) represents the number of accessions among all the accessions.The boxes show the medians and upper/lower quartiles.The whiskers extend to 1.5× the interquartile range, and any remaining points are indicated with dots.The differences between the haplotypes were statistically analyzed based on Tukey's test; NS, not significant.

Figure 12 .
Figure 12.Haplotype analysis of the candidate gene LOC_Os11g44680.(a) Local Manhattan plot (top) and linkage disequilibrium heatmap (bottom).(b) Schematic representation of the LOC_Os11g44680 structure and single-nucleotide polymorphisms in LOC_Os11g44680 cDNA.The blue boxes indicate exons.(c) Box plots for the RADR of the GE, RADR of the GR, RADR of the SL and RADR of the RL of LOC_Os11g44680 in all accessions in 2019 and 2020.The number of accessions (n) represents the number of accessions among all accessions.The boxes show the medians and upper/lower quartiles.The whiskers extend to 1.5× the interquartile range, and any remaining points are indicated with dots.The differences between the haplotypes were statistically analyzed based on Tukey's test; NS, not significant.

Figure 13 .
Figure 13.Haplotype analysis of the candidate gene LOC_Os11g44600.(a) Schematic representation of the LOC_Os11g44600 structure and single-nucleotide polymorphisms in LOC_Os11g44600 cDNA.The blue boxes indicate exons.(b) Box plots for the RADR of the GE, RADR of the GR, RADR of the SL and RADR of the RL of LOC_Os11g44680 in all accessions in 2019 and 2020.The number of accessions (n) represents the number of accessions among all accessions.The boxes show the medians and upper/lower quartiles.The whiskers extend to 1.5× the interquartile range, and any remaining points are indicated with dots.The differences between the haplotypes were statistically analyzed based on Tukey's test; NS, not significant.
Figure S2.Manhattan plots for alkalinity-tolerancerelated traits.(a) Manhattan plot for the GE in 2019; (b) Manhattan plot for the GR in 2019; (c) Manhattan plot for the SL in 2019; (d) Manhattan plot for the RL in 2019; (e) Manhattan plot for the GE in 2020; (f) Manhattan plot for the GR in 2020; (g) Manhattan plot for the SL in 2020; (h) Manhattan plot for the RL in 2020.Figure S3.Manhattan plots for alkalinity-tolerance-related traits.(a) Manhattan plot for the GE in 2019; (b) Manhattan plot for the GR in 2019; (c) Manhattan plot for SL in 2019; (d) Manhattan plot for the RL in 2019; (e) Manhattan plot for the GE in 2020; (f) Manhattan plot for the GR in 2020; (g) Manhattan plot for the SL in 2020; (h) Manhattan plot for the RL in 2020.

Figure S4 .
Correlation analysis between samples. 1, 2 and 3 represent biological replicates.Figure S5.Gene expression of the samples.(a) Gene expression abundance distribution map; (b) Gene expression violin plot.Table

Funding:
This research was funded by the National Natural Science Foundation of China, grant numbers 32101768 and U21A20214; the Talent Project of Anhui Agricultural University, grant number rc312002; the Natural Science Foundation of Anhui Province, grant number 2108085MC97; the Natural Science Foundation of Anhui Universities, grant number KJ2020A0118; the Natural Science Research Project of Colleges and Universities in Anhui Province, grant number YJS20210250; and the Key Research and Development Program of Anhui Province, grant number 202004a06020024.Data Availability Statement: Data are contained within the article and supplementary materials.

Table 1 .
Descriptive statistics for alkalinity-tolerance-related traits in 173 rice accessions at the germination stage.

Table 2 .
Loci significantly related to alkalinity-tolerance-related traits in 173 rice accessions at the germination stage.