Novel Genomic Regions of Fusarium Wilt Resistance in Bottle Gourd [Lagenaria siceraria (Mol.) Standl.] Discovered in Genome-Wide Association Study

Fusarium wilt (FW) is a typical soil-borne disease that seriously affects the yield and fruit quality of bottle gourd. Thus, to improve resistance to FW in bottle gourd, the genetic mechanism underlying FW resistance needs to be explored. In this study, we performed a genome-wide association study (GWAS) based on 5,330 single-nucleotide polymorphisms (SNPs) and 89 bottle gourd accessions. The GWAS results revealed a total of 10 SNPs (P ≤ 0.01, −log10P ≥ 2.0) significantly associated with FW resistance that were detected in at least two environments (2019DI, 2020DI, and the average across the 2 years); these SNPs were located on chromosomes 1, 2, 3, 4, 8, and 9. Linkage disequilibrium (LD) block structure analysis predicted three potential candidate genes for FW resistance. Genes HG_GLEAN_10001030 and HG_GLEAN_10001042 were within the range of the mean LD block of the marker BGReSe_14202; gene HG_GLEAN_10011803 was 280 kb upstream of the marker BGReSe_00818. Real-time quantitative PCR (qRT-PCR) analysis showed that HG_GLEAN_10011803 was significantly up-regulated in FW-infected plants of YD-4, Yin-10, and Hanbi; HG_GLEAN_10001030 and HG_GLEAN_10001042 were specifically up-regulated in FW-infected plants of YD-4. Therefore, gene HG_GLEAN_10011803 is likely the major effect candidate gene for resistance against FW in bottle gourd. This work provides scientific evidence for the exploration of candidate gene and development of functional markers in FW-resistant bottle gourd breeding programs.


INTRODUCTION
Bottle gourd [Lagenaria siceraria (Mol.) Standl.] (2n = 2 × = 22), also known as calabash or long melon, is a member of the Legendaria genus, Cucurbitaceous family, and is an annual plant (Whitaker, 1971;Erickson et al., 2005). Its fresh young fruits are eaten as a vegetable in parts of Asia and Africa, while dry old fruits are used for containers, musical instruments, and crafts (Heiser, 1979;Morimoto and Mvere, 2004). Moreover, due to its strong resistance to disease and abiotic stress, bottle gourd is commonly used as grafting rootstock for other crops, such as watermelon and melon (Lee, 1994;Yetisir and Sari, 2003). According to records, bottle gourd has been cultivated in China for more than 7,000 years, covering a cultivation area of more than 2 million mu; this is an important vegetable crop in China, especially in the southern part.
Fusarium wilt (FW), which is caused by Fusarium oxysporum, is a typical soil-borne disease of economic crops worldwide (Katan, 1994;Wechter et al., 2012;Bodah, 2017). Since this pathogen can survive in the absence of host-infected plants, once the disease occurs in the field, F. oxysporum is likely to remain in the soil indefinitely, which seriously affects the yield of crops (Cha et al., 2016;Khan et al., 2017). FW has a broad host on cucurbit crops, including watermelon, melon, cucumber, luffa, and bottle gourd (Bodah, 2017). Usually, the infected plants morphologically show a constriction in the stem xylem, resulting in vascular bundle clogging, plant wilting, or death (Freeman et al., 2002;Singh et al., 2017). FW commonly occurs during the whole growth period of bottle gourd, especially the flowering to fruiting period. Its incidence is approximately 20-40%, and severe cases could cause devastating losses (data from Ningbo Institute of Agriculture), which severely restrict the sustainable development of the bottle gourd industry.
Breeding resistant varieties is one of the most effective and economic methods to control FW disease. At present, a series of commercial varieties that are highly resistant to FW has been grown for production, such as watermelon, melon, and cowpea (Zink and Gubler, 1985;Martyn and Bruton, 1989;Ehlers et al., 2000Ehlers et al., , 2009. To improve FW resistance, we need to exploit markers tightly linked to FW resistance using quantitative trait loci (QTL) and then generate germplasm by molecular marker assistant selection (MAS; Zhao et al., 2014;Li et al., 2017). To date, QTLs/genes conferring FW resistance have been thoroughly studied in many crops. For example, a dominant gene I-2 that confers resistance to race 2 of FW was cloned in tomato by mapbased cloning (Simons et al., 1998;Catanzariti et al., 2015). Using the same map-based cloning technique, Joobeur et al. (2004) identified two candidate genes of melon FW resistance in the physical range of Fom-2. Several genes associated with cowpea FW resistance were identified using QTL analysis in "California Blackeye 27" (Pottorff et al., 2012(Pottorff et al., , 2013. In cucumber, a major QTL Foc2.1 conferring resistance to FW was detected using recombinant inbred lines (Zhang et al., 2014), and a major QTL, Fo-1.1, associated with FW resistance to race 1 was identified by using selective genotyping in genetic populations derived from related watermelon lines (Lambel et al., 2014). However, research progress on the FW resistance of bottle gourd is relatively limited. Only its specialization of F. oxysporum f. sp. lagenariae has been reported, whereas the genetic mechanism of FW resistance and related genes/QTLs are unknown in bottle gourd.
To date, a high-density genetic map has been constructed, and a series of ISSR, SSR, and single-nucleotide polymorphism (SNP) markers has been exploited for bottle gourd (Xu et al., 2011(Xu et al., , 2014Bhawna et al., 2014), allowing the establishment of various marker-trait associations, such as association analysis for the free glutamate content of bottle gourd . Genomewide association study (GWAS), based on linkage disequilibrium (LD), has also been widely used in the study of plants, and various results have been reported (Joobeur et al., 2004;Wang et al., 2009;Sabbavarapu et al., 2013;Zhang et al., 2014). In bottle gourd molecular breeding, Wu et al. (2017) performed a GWAS for SNPs related to the free glutamate content of the umami factor. With the development of quantitative genetics, many researchers have proposed different analytical models, such as efficient mixed-model association (Kang et al., 2008), compressed mixed linear model (Zhang et al., 2010), restricted two-stage multi-locus GWAS (He et al., 2017), etc. Among them, general linear model (GLM) and mixed linear model (MLM) are still the common GWAS methods in plants Li et al., 2013;Fang et al., 2017). In this study, we initially genotyped 89 bottle gourd accessions using 5,330 SNPs and surveyed the disease index (DI) of FW resistance in two consecutive years. We then performed a GWAS to identify significant associated SNPs and potential candidate genes. Finally, three candidate genes associated with FW resistance were verified by quantitative realtime PCR (qRT-PCR). Our study is the first to use GWAS to identify genomic regions and candidate genes associated with FW resistance. The GWAS results can lay a foundation for MAS breeding and the genetic mechanisms of FW resistance in cucurbit crops.

Plant Materials
Germplasm consisting of 89 bottle gourd accessions was collected, consisting of 87 accessions from wide areas across China, one accession from Europe, and one accession from Mexico (Supplementary Table 1). All accessions (inbred lines) were evaluated for FW resistance in a glasshouse of the Haining Experimental Station (30 • N, 120 • E) in 2019 and 2020. According to a completely randomized block design, the plants were studied based on two replications in both years.

Inoculation System of FW Resistance in Bottle Gourd
During 2018-2019, bottle gourd FW fungus was isolated from wilted plants that were collected from severely affected areas such as Shaoxing and Haining (Supplementary Figures 1A,B). According to the conventional tissue separation method, FW strains with obvious antagonistic effects were obtained by using potato dextrose agar (containing 100 mg/ml Kana and Amp) for screening four to five times (Supplementary Figure 1C). Under a microscope with 10 × 40 magnification, small conidia were observed to be ovoid or kidney-shaped, and large conidia were spindle-shaped or sickle-shaped with hooked apex (Supplementary Figure 1D). PCR assays showed that the similarity between the sequence of FW isolates and the 16S rRNA sequence was as high as 99%. After cytological tests and PCR detection, the isolates were identified as F. oxysporum f. sp. strains and were stored at 4 • C at the Zhejiang Academy of Agricultural Sciences, Hangzhou, China.
Each FW strain was shake-cultured on potato sucrose broth for 3 days in the dark at 28 • C at 200 rpm. With the use of a hemacytometer, the conidial suspension was adjusted to a final concentration of 1.0 × 10 6 conidia/ml with sterile distilled H 2 O. The seeds of each accession were sown in mixed soil (nutritional soil/vermiculite = 3: 1) in plastic pots (6 by 6 by 5 cm) and were grown in a glasshouse set at 24 • C, 16h light/18 • C, 8-h darkness, 60% humidity. At the second true leaf of the seedling spreading stage, we used the root dipping method for bottle gourd FW resistance screening and testing. Each accession consisted of 10-12 seedlings, and two duplicates were set per environment.

Disease Assessment and Statistical Analysis
Leaf damage was used as a main index to evaluate resistant/susceptible phenotypic traits. The standard reported by Gao et al. (2015) and Xu et al. (2016) was further improved and implemented with a few modifications. We classified the phenotypes of plants according to a 0-4 scale as follows: level 0, no disease symptoms, i.e., immune (I); level 1, slight disease symptoms, featured by less than 25% of leaves with disease symptoms, with normal plant growth, i.e., highly resistant (HR); level 2, slight wilt symptoms, featured by 25-50% of leaves with disease symptoms, i.e., resistant (R); level 3, moderate wilt symptoms, featured by 50-90% of leaves with disease symptoms, i.e., susceptible (S); and level 4, completely wilted or dead plants, i.e., highly susceptible (HS; Supplementary Figure 2). After 10-12 replicates per material were evaluated individually, we calculated the mean value to determine the disease severity for each accession. The DI was calculated according to the following equation: where DI is the disease index, Pi is the grade of the DI, ni is the plant number of the corresponding DI grade, N is the total number of plants investigated, and P4 is the highest DI grade. According to the DI scores, the FW resistance of each material was determined following Shen and Li (2008) with a few modifications: immune (DI = 0, level 0), highly resistant (0 < DI ≤ 15%, level 1), resistant (15% < DI ≤ 45%, level 2), susceptible (45% < DI ≤ 75%, level 3), and highly susceptible (DI > 75%, level 4).

SNP Genotyping, LD, and Population Structure
The SNP markers used in this study were generated from RAD sequencing with paired-end sequencing (150 bp) on the Illumina HiSeq platform. We initially found 89 bottle gourd accessions that aligned to the Hangzhou gourd reference genome of ∼330 Mb (Wang et al., 2018) 1 and then removed those SNPs with a minor allele frequency (MAF) of ≤0.01 and a heterozygous rate ≥25% for data filtering. This left a total of 6,222 highquality SNPs. Of these SNPs, 85.66% were located on the 11 chromosomes of the bottle gourd, leaving 5,330 high-quality SNPs. These were used for the correlation analysis of traits . The density of SNP markers was estimated to be one SNP per 59.37 kb for the 11 bottle gourd chromosomes.
The LD parameters (r 2 ) for estimating the LD distance of the genome between pairwise SNPs (MAF > 0.01) were calculated using PLINK V1.07 (Purcell et al., 2007;Xu et al., 2021, unpublished), and the average LD decay was drawn with R. The population structure was constructed by STRUCTURE 2.3.4 software (Evanno et al., 2005). K (number of subgroups) values were set to 2-8, with 10,000 (MOMC, Markov chain Monte Carlo)-100,000 runs (MCMC, Monte Carlo Markov Chain) with four replications. Then, the best value of K was determined by Ln P(D) and Delta K according to the principle of maximum likelihood (Evanno et al., 2005). The neighborjoining tree was constructed using PHYLIP software. The kinship matrix was assessed based on the SNP dataset using TASSEL 5.2.14 to determine the relatedness among individuals (Anderson and Weir, 2007;Zhang et al., 2010). In previous studies, the population was divided into two subgroups depending on the markers used in the tests .

Genome-Wide Association Analysis and LD Block Construction
For natural populations, the population structure and relative kinship always lead to high levels of false positives in association maps (Yu et al., 2006). After assessment of the population structure (Q), relative kinship (K), and principal component analysis (PCA) of 89 accessions, four correlation analysis models including (1) a general linear model GLM (Q), GLM (PCA) and (2) a mixed linear model MLM (Q + K), MLM (PCA + K) were used to conduct a genome-wide correlation analysis of FW resistance using TASSEL 5.2.14 ( Anderson and Weir, 2007;Zhang et al., 2010). The significance threshold for SNP-trait associations was determined by 1/n, where n is the number of markers in the association panel , and P ≤ 1/5,330 or −log 10 P ≥ 3.7. Considering that population structure and kinship reduced the detection efficiency of SNPs associated with FW resistance, the −log 10 P value of significantly associated SNPs identified in this study was low, which has also appeared in previous studies (Atwell et al., 2010;Huang et al., 2010). In order to fully exploit the valuable genetic information in the bottle gourd germplasm population, the significant threshold for SNP-trait associations was set as −log 10 P = 2. This threshold has already been applied to other traits in an association analysis (Li et al., 2015;Zhang et al., 2018). The correlation analysis results were plotted using a Manhattan plot and Q-Q plot based on the "qqman" package in R software.
The genome-wide LD decay rate, defined as the LD block distance where the LD decays to half of its maximum value, was about 400 kb in a natural population of bottle gourd (from Xu et al., 2021, unpublished). We defined the average LD decay distance as the candidate region to select candidate genes associated with large-effect SNPs. The genome of "Hangzhou gourd" was used as a reference sequence (Wang et al., 2018). Based on the genomic annotations of GourdBase, 2 potential candidate genes for FW resistance were predicted.

Validation of Candidate Genes
The expression levels of the candidate genes were measured before and after infecting plants with FW by using qRT-PCR. Based on the phenotype data in 2019 and 2020, Hanbi (HR to FW, level 1), Yin-10 (HR to FW, level 1), and YD-4 (HS to FW, level 4) were chosen as extreme materials and were cultivated in the glasshouse. The leaves from healthy plants (CK) and treatment plants were collected 3 days after FW infection and stored in liquid nitrogen. Total RNA was extracted from Hanbi, Yin-10, and YD-4 leaves using an RNA Simple Total RNA kit (Tiangen, China). After the quality and concentration of total RNA were evaluated using 1% agarose gel and an Agilent 2100 Bioanalyzer, complementary DNA (cDNA) was synthesized by using a Script cDNA Kit (Tiangen, China) with a standard protocol. The CDS sequences of genes were obtained from the GourdBase website. 3 qRT-PCR primers (Supplementary Table 2) were designed using the Integrated DAN Technologies website 4 and were synthesized by Sangon Biotech (Shanghai) Co., Ltd. The bottle gourd TuB-α gene (BG_GLEAN_10019523) was used as the internal control gene. qRT-PCR was performed on a Bio-Rad CFX96 Touch q-PCR System (Bio-Rad, CA, United States) with SuperReal PreMix Plus/SYBR Green (Tiangen, China). Each reaction was replicated three times. The relative expression level of candidate genes was evaluated by the 2 − Ct method (Livak and Schmittgen, 2001); healthy plants (CK) served as the control.

Phenotypic Analysis of FW Resistance in the Natural Population
In this study, a total of 89 bottle gourd accessions were evaluated for resistance to Fusarium wilt in a glasshouse in 2019 (DI2019) and 2020 (DI2020), with two replicates per environment. DI, as an evaluation of FW resistance, had a wide range of phenotypic variation in the 2-year trials. The DI of all accessions ranged from 6 to 95%, with a mean value of 46% in 2019 (DI2019), and from 11 to 94%, with a mean value of 55% in 2020 (DI2020). The ANOVA results showed that the broad-sense heritability (h 2 ) was 87.19% across the 2 years (Table 1), suggesting that the genetic effect played a predominant role in the phenotypic variation of FW resistance in bottle gourd. We divided the DI into five levels (Supplementary Figure 2): immune (level 0), highly resistant (level 1), resistant (level 2), susceptible (level 3), and highly susceptible (level 4), according to relevant previous studies (Gao et al., 2015;Xu et al., 2016). Only a tiny percentage of accessions had DI values less than 15% (8 in 2019 and 1 in 2020), whereas the majority of the accessions were within the range of 15.01-45% (35 in 2019 and 28 in 2020) and 45.01-75% (37 in 2019 and 42 in 2020). When DI exceeded 75%, there were 9 accessions in 2019 and 18 accessions in 2020 (Figure 2). Unfortunately, we did not select for any material that was immune to FW in the 2-year trials; only a small amount of material had high resistance to FW (Figure 2), which showed that the bottle gourd germplasm resource of FW resistance is scarce. The correlation coefficient between the 2-year trials was as high as 0.62 (Supplementary Figure 3), and the frequency distribution of DI was approximately normally distributed, which indicated that this natural population could be suitable for correlation analysis for FW resistance.

SNP Marker Analysis
The SNP markers used in this work resulted from RADsequencing by using the Illumina HiSeq platform. After removal of SNPs with a MAF of ≤0.01 and a heterozygous rate ≥25%, a total of 5,330 high-quality SNPs were retained for GWAS of the FW resistance trait. These SNPs covered all 11 chromosomes, with an uneven distribution across the genome ( Table 2). The average density of SNP markers was about 59.37 kb/SNP. The maximum marker density was found on chromosome 11 (101.18 kb/SNP) followed by chromosome 6 (67.25 kb/SNP), whereas the minimum marker density was found on chromosome 1 (42.11 kb/SNP). Based on the SNP markers, we estimated the population structure of 89 bottle gourds using STRUCTURE software and cluster analysis. The   delta K reached a sharp peak when K was 2. Therefore, this association population was divided into two subgroups, namely, subgroup 1 and subgroup 2 (Figures 3A,C). Subgroup 1 contained 80 accessions, and subgroup 2 was small and included nine accessions. A neighbor-joining result also classified the population into two subgroups, consistent with the population grouping result (Figure 3D). Because all accessions have some distant relationship, there were no primary factors, such as geographic distribution, affecting the population structure of the 89 accessions. Genotype data were subjected to correlation analysis of the free glutamate content trait in bottle gourd .
To determine the mapping resolution of the GWAS, we estimated the genome-wide LD decay distance of the association population. The average LD decay distance was approximately 400 kb, where r 2 = 0.3 for all chromosomes [half of its maximum value, from Xu et al., 2021, unpublished] (Figure 3B).

Model Comparison for Correlation Analysis
To reduce a false positive association, we applied four kinds of association analysis models to GWAS for FW resistance in bottle gourd. Based on the mean value of DI across the 2 years, quantilequantile (Q-Q) plots were drawn (Supplementary Figure 4). The results showed that there was a large deviation in the −log 10 P value between the observed values and the expected values in GLM (PCA) and GLM (Q) models, which indicated that the two models might cause a high false positive correlation. Due to the introduction of the covariable K, the observed −log 10 P values fit well with the expected values in the MLM (PCA + K) and MLM (Q + K) models, indicating that those two models could effectively control the false positive of the association analysis results. Taking into account the Q-Q plots of each environment, the MLM (Q + K) model (red scatter plot in Supplementary  Figure 4) was selected for the follow-up association analysis for FW resistance.

Genome-Wide Association Analysis
A GWAS was performed to detect SNPs associated with FW resistance between 5,330 SNP markers and 89 phenotype data points from the mean across the 2 years (aDI) and within an individual year (DI2019 and DI2020). The Manhattan plots and Q-Q plots for the GWAS results are shown in Figure 4. The GWAS result showed that 20 SNPs (with a significance threshold of p ≤ 0.01, −log 10 P ≥ 2.0) significantly associated with FW resistance were detected in at least one environment (Supplementary Table 3), including 12 SNPs from the 2019 data, 11 SNPs from the 2020 data, and 11 SNPs from the mean data. Among these SNPs, 10 significantly correlated SNP sites were detected in at least two environments, which were located on chromosomes 1, 2, 3, 4, 8, and 9, indicating that the FW resistance of bottle gourd is controlled by multiple genes (Figures 4A-C). The phenotypic variation explained by these sites ranged from 8.82 to 15.03% (Table 3). Among them, markers of BGReSe_14212 and BGReSe_14202 were located on chromosome 9, and those two SNP markers were within the range of the genome-wide LD block (400 kb). BGReSe_14202 was detected in all three environments with relatively high significant levels (−log 10 P = 2.81/2.49/2.46) and an effect on FW (R 2 = 14.14%/13.90%/10.49%). Therefore, the region range of chromosome 9 may contain the major genes associated with FW resistance. On chromosome 8, two SNP markers were detected with a certain LD distance away. BGReSe_12911 was significantly correlated with FW resistance in all three environments, and BGReSe_12338 was detected in DI2019 and aDI. Two SNP markers, BGReSe_02569 and BGReSe_02108, were detected on chromosome 2. Of these two, BGReSe_02569 explained the largest phenotypic variation in DI2020 and aDI, i.e., 16.19 and 15.38%, respectively. BGReSe_02108 was detected in three environments, with a contribution rate for phenotypic variation of 12.60, 11.03, and 11.28%. BGReSe_01042 and BGReSe_00818 were located on chromosome 1. One of the markers, BGReSe_00818 (−log 10 P = 2.25/2.02), was significantly correlated with FW resistance in the two environments of DI2019 and aDI, and its contribution rate for phenotypic variation was 12.26 and 12.84%, respectively ( Table 3).

Prediction of Candidate Genes for FW Resistance
In this study, we were interested in the markers with the greatest effect, such as marker BGReSe_00818 (MAF = 1.16) on chromosome 1 and markers BGReSe_14202 (MAF = 1.09) and BGReSe_14212 (MAF = 1.07) on chromosome 9. To reduce the range of candidate regions, we performed LD block structure analysis. The results showed that BGReSe_14202 and BGReSe_14212 could form an obvious LD (±400 kb) block, meaning that these two SNPs were closely linked ( Figure 5A). The candidate region of chromosome 9 was narrowed down to approximately 140 kb. This region contained 15 genes, of which two candidate genes were significantly associated with FW resistance of bottle gourd ( Table 4). Both of them, HG_GLEAN_10001030 (ethylene-responsive transcription factor RAP2) and HG_GLEAN_10001042 (GDSL esterase), are involved in signaling pathways, such as resistance genes and hormone induction. LD block reduced the candidate region of BGReSe_00818 to about 415 kb, which contained seven genes ( Figure 5B). Among them, HG_GLEAN_10011803 encodes carboxylesterase and a CDPK-related kinase protein and plays a role in the signal transduction pathway. To confirm whether the potential candidate genes participated in the FW resistance pathway, the expression patterns of the three genes in both FW-infected and healthy bottle gourd plants were analyzed via qRT-PCR. The representative materials were selected from the association analysis population in this study. The DI of Yin-10 and Hanbi was 10.83 and 13.44%, respectively, and both were TABLE 3 | Significant markers associated with Fusarium wilt resistance in at least two environments.

Marker
Chromosome Positive Allelic 2019DI 2020DI aDI −log 10 P R 2 (%) −log 10 P R 2 (%) −log 10 P R 2 (%)  highly resistant (HR, level 1) to FW. The DI of YD-4 was 87.81%, i.e., highly susceptible (HS, level 4) to FW. The expression pattern of three potential candidate genes HG_10011803, HG_10001030, and HG_10001042 in materials YD-4, Yin-10, and Hanbi is presented (Figure 6). Compared to healthy YD-4 (HS material, level 4), the expression levels of the three candidate genes were all significantly higher (P < 0.001) in the FW-infected group (3 days after infection) ( Figure 6A). For Yin-10 and Hanbi (HR materials, level 1), the expression level of gene HG_10011803 showed a significant difference (P < 0.05 and P < 0.001) between FW-infected and healthy groups. However, the expression levels of HG_10001030 and HG_10001042 in infected plants showed a higher or lower expression level than those in healthy Yin-10/Hanbi plants, without a significant difference (Figures 6B,C). Combining the above-mentioned interesting results, we speculated that HG_10011803 is a major effect gene, whereas HG_10001030 and HG_10001042 might be the candidate genes involved in the FW resistance response in bottle gourd.

DISCUSSION
Fusarium wilt is one of the most important diseases throughout the world, which seriously affects the yield and quality of cucurbit crops (Miguel et al., 2004;Wechter et al., 2012;Oumouloud et al., 2013). The genetic mechanism of resistance to FW in cucurbit crops is complex, showing genetic diversity. However, there are no studies on the genetic effect and inheritance of genes governing FW resistance in bottle gourd, and molecular markers linked to FW resistance are also poorly reported. Genome-wide association study has emerged as a powerful tool to study complex traits and genetic variations in SNP loci and has been successfully applied to different crops in recent  years (Joobeur et al., 2004;Wang et al., 2009;Sabbavarapu et al., 2013;Zhang et al., 2014). Especially the general linear model and mixed linear model are still the common GWAS methods in plants Li et al., 2013;Fang et al., 2017). GLM, based on a linear regression model, is usually used for the analysis of quantitative traits and discrete resistance traits. MLM, based on population structure (Q) and kinship (K) as covariance, could better reduce the false positive association (Yu et al., 2006). Taking into account the deviation between expected −log 10 P and observed −log 10 P in Q-Q plots, we finally selected the MLM (Q + K) (Supplementary Figure 4) as GWAS model for FW resistance. In this study, we used a GWAS to evaluate a population of 89 accessions for FW resistance under glasshouse inoculation conditions. A total of 20 SNP S (P ≤ 0.01, −log 10 P ≥ 2.0)significantly associated with FW resistance were identified in at least one environment (Supplementary Table 3). These sites were distributed on seven chromosomes, which could explain the phenotypic variation up to 16.19%. Among them, 10 significantly correlated SNP sites were detected in at least two environments, which were located on chromosomes 1, 2, 3, 4, 8, and 9 (Figure 4). According to the reference genome sequence of "Hangzhou Gourd" (Wang et al., 2018), we preliminarily predicted three candidate genes in candidate regions or LD block regions of these 10 SNP markers (Figure 5). HG_GLEAN_10011803, a candidate gene, which was located 280 kb upstream of the BGReSe_00818 marker on Chr.1, encodes calcium-dependent protein kinase (CDPK) protein.
There have been increasing studies confirming the involvement of CDPKs in plant disease resistance defense responses (Boudsocq and Sheen, 2013). For example, Loss-AtCPK28 or overexpression-AtCDPK1 mutants displayed enhanced responses to antibacterial immunity in Arabidopsis (Coca and Segundo, 2010;Monaghan et al., 2014). SlCRK6 in tomato played a role in resistance to both Sclerotinia sclerotiorum and Pseudomonas syringae pv. tomato (Pst) DC3000 (Wang et al., 2016). StCDPK5VK in potato could increase resistance to late blight fungus through the production of ROS (Kobayashi et al., 2012). In addition, by conducting a qRT-PCR analysis, we found that the expression level of HG_GLEAN_10011803 in FW-infected plants was significantly higher than that in healthy plants (Figure 6). Therefore, we inferred that the candidate gene HG_GLEAN_10011803 might be related to the FW resistance of bottle gourd.
In the LD block region of another candidate marker BGReSe_14202, one candidate gene HG_GLEAN_10001030, located 50 kb upstream of this marker on chromosome 9, encoded the ethylene-responsive transcription factor (ERTF) RAP2 protein. ERTFs play an important regulatory role in plant signal transduction of disease resistance and stress resistance, and overexpression could improve plant disease resistance and stress resistance (Singh et al., 2002;Gutterson and Reuber, 2004). For example, OsRAP2.6-overexpressed plants showed improved resistance to rice blast fungus (Wamaitha et al., 2012). TERF1 and TSRF1 genes in tomato could be resistant to Ralstonia solanacearum and Botrytis cinerea Zhang et al., 2004Zhang et al., , 2008. Another candidate gene, HG_10001042, located 18 kb downstream of this marker on chromosome 9, is a member of the GDSL gene family. The GDSL gene family consists of a wide range of members and plays important roles in plant growth, development, and stress defense responses (Akoh et al., 2004;Chepyshko et al., 2012). Overexpressed GDSL genes, such as AtGLIP1 and CaGLIP1, could enhance the resistance to a variety of pathogenic fungi (Hong et al., 2008;Lee et al., 2009;Naranjo et al., 2010). The qPR-PCR results showed that the expression levels of these two candidate genes were significantly increased in FW-infected YD-4 (HS material, level 4), while their expression levels were not significantly different before and after infection of Yin-10/Hanbi (HR materials, level 1) (Figure 6). Thus, we postulate that these three genes were candidate genes for FW resistance; in particular, HG_GLEAN_10011803 might be a major effect gene. However, further evidence is needed to functionally validate this hypothesis. To our knowledge, this study is the first to perform GWAS for FW resistance in cucurbit crops. Our results provide the molecular tools for FW resistance selection and lay a foundation for candidate gene discovery.
The resistant materials and SNP markers that we identified will promote breeding programs for FW-resistant bottle gourd.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
YL and GL conceived and designed the research. YL, YW, and JW performed the experiments. YL wrote the manuscript. YL, YW, and XinW constructed the population and collected phenotypes. XiaW, BW, and ZL carried out the field work. All authors analyzed the data and approved the submitted version.