A GWAS identified a major QTL for resistance to Fusarium wilt (Fusarium oxysporum f. sp. vasinfectum) race 4 in a MAGIC population of Upland cotton and a meta-analysis of QTLs for Fusarium wilt resistance

A major QTL conferring resistance to Fusarium wilt race 4 in a narrow region of chromosome D02 was identified in a MAGIC population of 550 RILs of Upland cotton. Numerous studies have been conducted to investigate the genetic basis of Fusarium wilt (FW, caused by Fusarium oxysporum f. sp. vasinfectum, FOV) resistance using bi-parental and association mapping populations in cotton. In this study, a multi-parent advanced generation inter-cross (MAGIC) population of 550 recombinant inbred lines (RILs), together with their 11 Upland cotton (Gossypium hirsutum) parents, was used to identify QTLs for FOV race 4 (FOV4) resistance. Among the parents, Acala Ultima, M-240 RNR, and Stoneville 474 were the most resistant, while Deltapine Acala 90, Coker 315, and Stoneville 825 were the most susceptible. Twenty-two MAGIC lines were consistently resistant to FOV4. Through a genome-wide association study (GWAS) based on 473,516 polymorphic SNPs, a major FOV4 resistance QTL within a narrow region on chromosomes D02 was detected, allowing identification of 14 candidate genes. Additionally, a meta-analysis of 133 published FW resistance QTLs showed a D subgenome and individual chromosome bias and no correlation between homeologous chromosome pairs. This study represents the first GWAS study using a largest genetic population and the most comprehensive meta-analysis for FW resistance in cotton. The results illustrated that 550 lines were not enough for high resolution mapping to pinpoint a candidate gene, and experimental errors in phenotyping cotton for FW resistance further compromised the accuracy and precision in QTL localization and identification of candidate genes. This study identified FOV4-resistant parents and MAGIC lines, and the first major QTL for FOV4 resistance in Upland cotton, providing useful information for developing FOV4-resistant cultivars and further genomic studies towards identification of causal genes for FOV4 resistance in cotton.


Introduction
Fusarium wilt of cotton, caused by soilborne pathogen Fusarium oxysporum f. sp. vasinfectum Atk. Sny & Hans (FOV), is a devastating disease that can lead to significant yield losses in most cotton-producing regions of the world. In the US, FW caused annual cotton (Gossypium spp.) yield loss ranging from 0.19 to 1.36% between 1952 and 2012 (Blasingame and Patel 2013), and the yield loss still remained at high levels in several states including Alabama, Arizona, Arkansas, Georgia and New Mexico at 0.4-1.0% and in California at 2.5% in 2020 (Lawrence et al. 2021). The pathogen infects cotton roots and colonizes the vascular system, causing plant wilting, leaf chlorosis and necrosis, stunting, vascular discoloration, defoliation, and eventual plant death (Zhang et al. 2015a;Sanogo and Zhang 2016). FOV has been characterized into eight pathogenic races (designated races 1 to 8) based on DNA sequence analysis and pathogenicity tests (Cianchetta et al. 2015;Davis et al. 2006). FOV race 4 (FOV4) is an emerging pathogen that threatens cotton production in the west and southwest and potentially the entire US Cotton Belt (Kim et al. 2005;Hutmacher et al. 2013;Halpern et al. 2018;Zhu et al. 2020). It is highly virulent to both Pima (G. barbadense L.) and Upland (G. hirsutum L.) cotton and infects cotton in the early seedling growth stages under field conditions, causing seedling wilt and deaths in the absence of nematodes in many cotton fields. It has been a significant problem in California for the past 20 years (Kim et al. 2005;Diaz et al. 2021) and has since been reported in El Paso, Texas in 2017 (Halpern et al. 2018;Bell et al. 2019;Diaz et al. 2021;Zhu et al. 2019Zhu et al. , 2021b and Las Cruces, New Mexico in 2019 (Zhu et al. 2020;Zhu et al. 2021a). Four different genotypes of FOV4 isolates were recently identified through a multiplex PCR method based on the presence of Tfo1 (T type), Tfo1/ MULE (MT type) and Tfo1/MITE (MiT type) or absence of above insertions (N type) in the PHO gene encoding for phosphate permease (Bell et al. 2019). Managing the disease is extremely difficult because FOV can persist in soil for many years or decades, and the most and only effective way to manage FOV4 is to breed and utilize resistant cotton cultivars (Zhang et al. 2015a;Sanogo and Zhang 2016).
Breeding for resistant cotton cultivars to FW is challenging. Many efforts have been made to enhance FW resistance in cotton germplasms or commercial cultivars in the US (Bird 1982;Kappelman and Bird 1981;Thaxton and El-Zik 1998, 2003, 2004aZhang et al. 2015a). For example, through the multi-adversity resistance (MAR) cotton system in Texas, progressive improvements in resistance to FW were made in the MAR germplasms , 2003, 2004a. In Arkansas, numerous germplasm lines and several cultivars with FW (caused by FOV race 1 and/or 2) resistance including Arkot 518 (Smith 1988) and H1330 (Bourland 1996) were released through a modified MAR system. Auburn University has also released cotton germplasms or commercial cultivars that possess some levels of FW resistance through testing in the National Fusarium Wilt Nursery. As a result, breeding for FW (caused by FOV race 1 and/or 2) resistance has been substantially enhanced in many cotton cultivars.
However, most commercial cultivars and breeding lines which were developed in non-FOV infested fields are susceptible to FW. The response to FOV infections in cotton has been treated as both a qualitative and a quantitative trait to study the genetic basis of resistance using traditional Mendelian and quantitative genetic approaches (Zhang et al. 2015a). In the 1930-40s, India reported that resistance to FOV4 was determined by two complementary dominant genes and a third inhibitory gene in G. arboreum L. and one dominant gene in G. herbaceum L. (cited by Smith and Dick 1960). Studies in other countries showed resistance to race 3 and 6 was controlled by a major gene (Fahmy 1934;Netzer 1982;Zhang et al. 2015a, b, c). In the US, a major dominant resistance gene (FOV1) to race 1 was confirmed in Pima S-7 (G. barbadense L.) and later mapped on chromosome D07 (c16) using simple sequence repeat (SSR) markers (Wang and Roberts 2006;Ulloa et al. 2011). A major resistance gene (FOV4) conferring resistance to FOV4 was identified in resistant Pima S-6 and mapped to chromosome D02 (c14) based on SSR markers . Also, two major resistance genes (dominant Mi1 and additive Mi2) to root-knot nematodes [RKN, Meloidogyne incognita Kofoid and White (Chitwood)] were identified and mapped to chromosomes A11 (c11) and D02 (c14), respectively, based on SSR markers (Niu et al. 2007;Gutierrez et al. 2010;He et al. 2014). Because RKN can increase the incidence of FW in both resistant and susceptible cultivars as FW caused by races 1 and 2 require the presence of RKN, it may be possible to breed for FW resistance by breeding solely for high RKN resistance. However, resistance to the soilborne FOV pathogen is not always associated with the RKN resistance. Furthermore, because FOV4 infections on cotton do not require the presence of RKN in the soil, it is currently unknown if RKN would reduce resistance to FOV4 in resistant lines and if one of the two RKN resistance genes is located within the same chromosomal region on D02 (c14) for FOV resistance. In China, two major dominant resistance genes (Fw 1 and Fw 2 ) to FOV race 7 were identified in resistant Upland cottons (Feng et al. 1998). Subsequently, a dominant resistance gene (Fw R , likely Fw 1 or Fw 2 ) on chromosome D03 (c17) was identified using SSR markers in two F 2:3 populations (Wang et al. 2009). A major resistance QTL (Fov7) on D03 (likely Fw R ) was recently identified to be Gh_ D03G0209 (named GhGLR4.8) encoding a GLUTAMATE RECEPTOR-LIKE (GLR) protein through a GWAS of 290 Chinese Upland accessions (Liu et al. 2021). Interestingly, a major resistance QTL to race 7 on D03 with two candidate genes (Gbar_D03G001430 and Gbar_D03G001910) was also recently reported in G. barbadense based on GWAS of 336 accessions (Zhao et al. 2021). However, the two QTLs and their candidate genes on D03 were different from one another and also different from another candidate gene for another race 7-resistant QTL and its candidate gene on the same chromosome [Gb_D03G0217, named GbCML encoding for a calmodulin (CaM)-like (CML) protein] identified in a G. barbadense RIL population of 110 lines through linkage mapping (Han et al. 2022). Previously, Zhu et al. (2010) mapped a single dominant resistance gene on chromosome D11 (c21) in a G. barbadense cross of susceptible Xinhai 21 × highly resistant HK 237 based on a limited number of SSR markers. The two major resistance genes for race 1 (on D07 from Pima S-7) and race 4 (on D02 from Pima S-6) are different from the two major resistance genes for race 7 (one on D03 from Upland and another on D11 from G. barbadense).
The response to FW in cotton could be treated as a qualitative or a quantitative trait depending on the specific cotton germplasm and race of the fungal pathogen. Quantitative trait loci (QTL) mapping is a powerful tool that has been widely used to dissect the genetic basis of complex quantitative traits in cotton. Identifying QTLs and dissecting the genetic basis of FOV resistance are important for marker-assisted selection (MAS) in cotton breeding. To date, many types of molecular markers, such as restriction fragment length polymorphism (RFLP), amplified fragment polymorphism (AFLP), simple sequence repeats (SSR), and single nucleotide polymorphisms (SNP), have played an important role in detecting resistance QTL locations or resistance genes for different races of FOV using biparental and association populations (Zhang et al. 2015a). In a metaanalysis, 33 QTLs mapped for FW resistance in cotton was summarized by Zhang et al. (2015b) from previous published papers based on early segregating populations and recombinant inbred lines (RILs). Based on a more recent meta-analysis, Abdelraheem et al. (2017) increased the number of QTLs reported for FW resistance to 47, and most FOV resistance QTLs were located on five chromosomes, i.e., A06 (c6), D02 (c14), D03 (c17), D04 (c22), and D06 (c25). Using a RIL population from an interspecific cross of Pima S-7 (resistant to FOV1 and susceptible to FOV4) × Acala NemX (susceptible to FOV1 and tolerant to FOV4 in the greenhouse) with 367 SSR markers, Wang et al. (2018) detected six major QTLs on chromosomes A01 (c1), A02 (c2), A12 (c12), D01 (c15), and D11 (c21) for FOV1 resistance and two major QTLs on D02 (c14) and D03 (c17) for FOV4 resistance. Most recently, 13 FOV4 resistance QTLs were identified on six chromosomes, i.e., A08 (c8), D02 (c14), D03 (c17), D05 (c19), D07 (c16), and D13 in the US Upland cotton with 25,677 SNP markers through a genomewide association study (GWAS) of 367 US Upland accessions (Abdelraheem et al. 2020a). Seven QTLs for FOV4 resistance were mapped on six chromosomes (D02/c14, D03/c17, D05/c19, D06/c25, D08/c24, and D11/c21) using 163 RILs of a G. barbadense cross from FOV4-resistant Pima S-6 × susceptible 89590 with 403 SSR markers (Abdelraheem et al. 2020b). Although many QTLs conferring FW resistance have been mapped on all except for a few chromosomes including A04 (c4) and D10 (c20), major effect resistance QTLs are located on only a few chromosomes (Zhang et al. 2015a, b;Abdelraheem et al. 2017).
Based on biparental populations, a large number of QTLs for FW resistance have been reported from the linkage analysis, and many QTLs may be commonly detected across multiple tests in the same population or across different populations (Zhang et al. 2015a, b;Abdelraheem et al. 2017). However, limited by low genetic diversity and marker density on a genetic map, linkage analysis can only detect and map a few QTLs in each genetic population (Huang et al. 2017;Li et al. 2017;Islam et al. 2016). Furthermore, limited numbers of markers and relatively small population sizes also reduce the resolution power to fine map FW resistance QTLs. In recent years, GWAS as a useful and powerful tool to detect QTLs has been widely used in QTL mapping in many plant species (Zhao et al. 2011;Atwell et al. 2010;Kump et al. 2011;Zhao et al. 2015). In cotton, GWAS has also been widely used in QTL mapping for important economic traits, such as disease resistance (Abdelraheem et al. 2020a;Zhao et al. 2021), fiber quality traits (Islam et al. 2015;Nie et al. 2016) and yield components Qin et al. 2015). Two recent GWAS reported identification of two major FOV7-resistance QTLs on D03 with each in Chinese Upland and G. barbadense using 290 and 336 accessions, respectively (Liu et al. 2021;Zhao et al. 2021). Abdelraheem et al. (2020a) conducted a GWAS using the CottonSNP63K array in 376 US Upland cotton accessions and detected two consistent QTL clusters on c16 (D07) and c19 (D05) for resistance to both Verticillium wilt (VW) and FOV4 in the greenhouse studies. However, reports of GWAS for FOV4 resistance in Upland cotton are still limited. As compared to linkage mapping, GWAS offers a higher resolution with relatively higher marker coverage on the genome and is more cost-effective for detecting important QTLs or genes associated with complex traits . With the use of high-throughput next-generation sequencing technology, such as genotyping by sequencing (GBS), high density of SNP markers were identified and used in Upland cotton (Islam et al. 2015) to improve the resolution and accuracy of QTL mapping. However, the mapping population size should also be sufficiently large to delineate a QTL to a narrow chromosome region. This will facilitate the identification of candidate genes underlying the natural variation in FW resistance, leading to the cloning of FW resistance genes.
In the US, early FW resistant cultivars developed were lower yielding with deficiencies in many other agronomic traits in non-infested fields than the best agronomically acceptable but less resistant cultivars. However, many 1 3 moderately or highly FW (not caused by FOV4) resistant cotton cultivars with high yield potentials and good quality have been developed since the late 1960s (Zhang et al. 2015a). Therefore, representative Upland cotton cultivars can be used to develop a multi-parent advanced generation inter-cross (MAGIC) population to capture more genetic variation for breeding and genetic and genomic studies. In this study, a GWAS was conducted for FOV4 resistance in a MAGIC population consisting of 550 RILs created from 11 US Upland cotton cultivars using a total of 473, 516 SNP markers developed through whole genome sequencing. The aim of this study was to identify candidate QTL regions and genes involved in FOV4 resistance and provide information on the genetic mechanisms of resistance to FOV4 in Upland cotton. This MAGIC population has been recently used to identify QTLs and candidate genes for fiber quality traits, resistance to RKN and VW, and tolerance to drought and salt (Islam et al. 2016;Thyssen et al. 2019;Naoumkina et al. 2019;Wubben et al. 2019;Zhang et al. 2020a;Abdelraheem et al. 2020a).

Materials and methods
In this study, a MAGIC population, consisting of 550 RILs designated MG1 to MG550, was developed from a random mated Upland population of 11 parents (Jenkins et al. 2008). The 11 parents represented diverse Upland cotton genetic backgrounds from major cotton seed companies and two public breeding programs in the US and also differed in FW resistance (Table 1). The MAGIC population used in this study was created through 5 generations of random-mating followed by 6 generations of self-pollination (Jenkins et al. 2008), and the development of this MAGIC population was described in detail by Islam et al. (2016).
To evaluate the population for FOV4 resistance, two tests (Test 1 and 2) each with 2 replications (with 10 seedlings for each genotype per replication) were arranged in a randomized complete block design (RCBD) in a greenhouse at the Fabian Garcia Research Center, New Mexico State University, Las Cruces, NM, USA. The seed of each genotype was planted in a 10-cm plastic pot with 10 seed per pot (five hills per pot, two seed per hill) at two different planting and evaluation dates for the two tests (May to July for Test 1 and September to October for Test 2). The temperature in the greenhouse ranged from 22 to 35 °C. The potting soil used was Miracle-Gro Moisture Control Potting Mix 2 CF (Scotts Co., Marysville, OH, USA). An established inoculation method was used to inoculate cotton seedlings with spores from a virulent local FOV4 isolate, FW-JF-1A Zhu et al. 2021a, b;Zhang et al. 2020c). FW-JF-1A was identified to be FOV4 based on a FOV4specific PCR, and it belonged to the MT genotype based on a multiplex PCR (Bell et al. 2019). Virulence on cotton between MT and other two FOV4 genotypes (N and T) was correlated (Zhu et al. 2021a, b). A total of 10 ml of 1 × 10 6 spores ml −1 conidial suspension of FOV4 were inoculated to emerged seedlings with 1-2 true leaf in each pot (10 plants pot −1 ). The inoculum was pipetted onto the soil surface with no root wounding, and pots were lightly watered immediately after inoculation. Fungal culture, spore quantification, and inoculation techniques followed Zhang et al. (2021a, b).
The response of the cotton seedlings to FOV4 infections was monitored daily and evaluated at 30 days after inoculation (DAI) using our previously established protocol with a 0-5 rating scale for foliar disease severity ratings (DSR) on an individual plant basis (Zhang et al. , 2020b, as follows: 1 No symptom 2 One wilted cotyledon 1 3 3 Two wilted cotyledons or two cotyledons abscised 4 First true leaf wilted or three leaves (including two cotyledons) abscised 5 Whole plant wilted or more than three leaves abscised 6 Dead plant The disease incidence (DI, percentage of infected plants with symptoms) and average DSR (the sum of the DSR divided by the number of plants) were calculated on a replication basis for each genotype. Since the two tests were from two different periods in the greenhouse, plot data were then subjected to a combined analysis of variance (ANOVA) for the two tests using SAS 9.3 (SAS Institute Inc., Cary, NC, USA). PROC MIXED statistical procedure was used to determine the statistical significance of various sources of variation, in which genotype was treated as a fixed effect, and test and replication as a random effect. Broad-sense heritability (h 2 ) in each test and across tests was estimated from variance components (Abdelraheem et al. 2018).
The genome of the MAGIC population was sequenced for SNP identification and GWAS for DSR-based FW resistance, and the details of methodology were described elsewhere for the same MAGIC population . Briefly, to genotype the MAGIC population, a whole genome 100-150 bp paired-end sequencing (Illumina) was used to sequence parents at 20 × and RILs at 3 × genome coverage by Novogene Corporation (Chula Vista, CA, USA) using Illumina HiSeq 2500. Sequence reads were aligned to the NBI G. hirsutum cv TM-1 reference genome (Zhang et al. 2015c) with updates using GSNAP software (Wu and Nacu 2010) for SNP identification using samtools and bcftools software (Li et al. 2009). As a result, 473,516 SNPs were identified. For GWAS, GAPIT software with default parameters (Lipka et al. 2012) was used to determine the association between SNPs and FOV4 resistance using a mixed linear model . In this study, a QTL was declared when 4 significant SNPs at 1 × 10 −3 were detected within a 6 Mb region. For a multiple comparison adjustment, a false discovery rate (FDR) adjusted P-value (i.e., the P value at 0.05 divided by the total number of 473,516 SNPs) was also calculated (Benjamini and Hochberg 1995). A QTL cluster was declared when multiple QTL were overlapped or were within a 15-40 Mb region.
To identify candidate genes for FOV4 resistance, annotated genes that were within the QTL regions from the GWAS were retrieved from the GFF files of the G. hirsutum cv TM-1 reference genome (Zhang et al. 2015c;Hu et al. 2019;Wang et al. 2019). To search for candidate genes, the annotated genes within the QTL regions were filtered based on relevant gene ontology terms of Arabidopsis orthologs, including response to fungus or biotic stimulus (Berardini et al. 2015). These candidate genes were then examined for non-synonymous variants. A SNP-index for each SNP within a QTL region was further calculated for the most FOV4-resistant RIL group and the most FOV4susceptible RIL group separately as the frequency of a SNP allele in each group. The difference in the SNP frequency between the two groups is called delta (SNP-index), i.e., D(SNP-index) (Takagi et al. 2013). If a SNP is not linked to a QTL, D(SNP-index) is close to zero. To assess if the functionally relevant genes with sequence variation within the QTL regions were differentially expressed, root tissue or germinating seed-specific RNA-seq data for these genes were retrieved from ccNET (http:// struc tralb iology. cau. edu. cn. gossy pium) . Then, genes that had D (SNP-index) equal to or higher than 0.3 and were expressed in relevant tissues and showed responses to the biotic stress stimuli were identified as candidate genes.

Genetic difference of FOV4 resistance among the 11 founder parents
On a genotype basis, DI ranged from 16.7 to 100% in Test 1 and from 11.8 to 87.5% in Test 2, while DSR ranged from 0.2 to 3.3 in Test 1 and from 0.2 to 2.6 in Test 2, among the 11 founder parents. Significant genotypic variation in FOV4 resistance was detected (Table 2). Tests 1 and 2 showed relatively consistent results in FOV4 resistance among the 11 founder lines, as indicated by a significant positive correlation in DSR between the two tests (r = 0.715, P <0.05). As shown in Table 2 , and therefore, these four genotypes were considered as the susceptible parent group. The existence of the quantitative genetic variation in FOV4 resistance among the 11 founder parents allowed the following analysis of the genetic and genomic basis of FOV4 resistance in the MAGIC population derived from them.

Analysis of variance and broad-sense heritability estimates for FOV4 resistance in the MAGIC population
As shown in Table 3, a combined ANOVA for the entire MAGIC population indicates that DI and DSR were significant (at P <0.001) in test, genotype, and also genotype × test interaction. Therefore, a further analysis was performed on an individual test basis. In both Tests 1 and 2, DI and DSR showed a significant genotypic variation (at P <0.001, except for DI in Test 1 at P <0.05). Test 1 averaged significantly higher DI and DSR than Test 2, due likely to relatively lower temperatures in the greenhouse during the test. Furthermore, the two tests were grown in the same greenhouse during different time periods with two different inoculum preparations and inoculations, in addition to other uncontrollable factors in greenhouse and crop management. These factors may have also contributed to significant differences between the two tests and genotype × test interactions, as seen from different results between the two tests for some parental genotypes (Table 2). Because DSR takes into consideration of incidence and leaf symptom severity including plant death, it was most reliable and therefore further used to quantitatively measure the overall FOV4 resistance in this study. Broad-sense heritability estimates for DSR in the two tests were similar, ranging from 0.60 in Test 1 to 0.63 in Test 2, and the estimate was 0.62 in a combined ANOVA from the two tests (Table 3). The results indicated that approximately 60% of the phenotypic variation in FOV4 resistance was due to genetic variation in this MAGIC population. Therefore, a further detection of genetic factors for FOV4 resistance was warranted.

Genetic difference of FOV4 resistance in the MAGIC population of 550 RILs
In the MAGIC population (Fig. 1), the DI ranged from 0 to 100% in both tests with an average of 89.4 and 83.6% in Tests 1 and 2, respectively; and the DSR ranged from 0 to 4.2 with an average of 1.8 in Test 1 and from 0 to 4.9 with an average of 1.5 in Test 2. An average of only 1.0 and 5.0% plant mortality was observed in Tests 1 and 2, respectively. Averaging across the two tests, many RILs had significantly and consistently  (Table 2), while the most susceptible RILs were more susceptible than the most susceptible parent (Deltapine Acala 90). The above results suggested that transgressive segregation toward susceptibility occurred during the intermating process among the progeny from the hybrids of the 11 founders. However, because responses of some of the RILs were likely different to FOV4 infections under optimal (low) temperature conditions (21-23 o C) in that higher plant mortality would be expected based on Zhang et al. (2020bZhang et al. ( , 2021a, RILs with significantly higher FOV4 resistance than the parents might be identified.

QTLs mapped for FOV4 resistance in the MAGIC population
Mean DSRs for each MAGIC line from Tests 1 and 2 and the average from the two tests were used for GWAS based on 473,516 polymorphic SNP markers in this study. Based on means across the two tests, 8 QTLs including 2 (one each on A10 and A12) on the A subgenome and 6 QTLs (5 on D02 and 1 on D04) on the D subgenome were detected (Table 5, Fig. 2). In Test 1, 10 QTLs were detected, including 3 QTLs with one each on A06, A10 and A13 of the A subgenome and 7 QTLs on the D-subgenome (5 on D02 and one each on D03 and D11). Notably, qFOV4-D02-3a detected between the 0.7 and 1.2 Mb region represented a major QTL with the highest delta(SNP-index) value and LOD score (Table 5, Fig. 2). In Test 2, 13 QTLs were detected, including 4 QTLs with one each on A05, A08, A12, and A13 of the A-subgenome and 9 QTLs with 6 on D02 and one each on D01, D11, and D13 of the D-subgenome. Three common QTLs were detected between Test 1 and the overall means across the two tests, while 3 QTLs were common between Test 2 and the overall means across the two tests. Therefore, a total of 6 common QTLs were detected between an individual test and the overall means across the two tests. Taken together, 22 QTLs (7 on the A-subgenome and 15 on the D-subgenome) for FOV4 resistance were detected on 12 of the 26 pairs of Upland cotton chromosomes. Two QTLs were detected on each of A12 and D11. However, the common QTL (qFOV4-A12-1) on A12 between Test 2 and the overall means across the two tests was not in proximity (9 Mb apart) to the QTL (qFOV4-A12-2) detected in Test 1. The 2 QTLs on D11 detected from Tests 1 and 2 were also different, as they were 37 Mb apart. Interestingly, two QTL clusters with a total of 9 QTLs were declared on D02. The two QTLs (qFOV4-D02-3 and qFOV4-D02-4) were overlapped within a 1.4 Mb interval (at the chromosome position between 0.4 and 1.8 Mb), which were only less than 35 kb away from qFOV4-D02-2 on the left and 400 kb apart from qFOV4-D02-5 on the right. Therefore, these 5 QTLs were clustered. The other 4 QTLs were 5 Mb apart (within the chromosome position 18-23 Mb) and were considered another QTL cluster on the same chromosome.

A meta-analysis of QTL reported for FW resistance in cotton
A summary of all reported FW resistance QTLs including the ones detected in this current study on chromosomes is shown in Table 6  On an individual chromosome basis, A04 and A11 had no FW resistance QTLs detected, and A02, A05, A07, A09, A10, and A13 had only 1-3 QTLs reported, while A03, A06, A08, and A12 had the most FW resistance QTLs (6-9) within the A subgenome. The resistance to FOV4 infections in the two A genome species-G. arboreum and G. herbaceum is inherited as a Mendelian trait (Zhang et al. 2015a); however, thus far, no major gene or QTL for FOV4 resistance has been reported on the A subgenome of Upland or Pima cotton. Du et al. (2018) reported a maineffect QTL on chromosome c11 (A11) for FOV7 resistance based on GWAS using 215 Chinese G. arboreum accessions. However, it is currently unclear whether this QTL confers resistance to FOV4, and no FOV resistance QTL has been reported on its homologous chromosome in tetraploid cotton. On the D subgenome, only D10 did not have reported QTLs for FW resistance, and all other chromosomes carried 3 or more FW resistance QTLs including D01, D02, D03, D05, D07, D11, D12, and D13 with 6 or more QTLs detected. Therefore, overall, the D subgenome harbors more QTLs for resistance to FW, and the FW resistance QTLs were not evenly distributed on the tetraploid cotton chromosomes. Furthermore, the chromosomes with more QTLs detected are not necessarily the ones possessing a major gene/QTL for FW resistance including FOV1 on D07 (c16), Fw R on D03 (c17), and an unnamed gene on D11 (c21) with the exception of FOV4 on D02 (c14).
On a homeologous chromosome pair basis, the two pairs (A03/D03 and A12/D12) had more QTLs detected, while A10/D10 had almost no QTL reported among the 13 pairs. A06 and A08 had 2.3-3.0 times more QTLs detected than their homeologous counterparts (D06 and D08), respectively. However, D01, D02, D03, D04, D05, D07, D11, and D13 carried 1.8 times or more QTLs than their homeologous A subgenome counterparts, respectively. Overall, there was 1 3 no significant positive correlation (r = 0.018, P >0.05) in the number of FW resistance QTLs between the 13 pairs of homeologous chromosomes, suggesting different mechanisms evolved from the two subgenomes to fight FOV infections in tetraploid cotton.

Identification of candidate genes for FW resistance QTL
In this study, a total of 532 genes were identified within the 22 QTL regions declared, of which 198 genes were found to have non-synonymous SNPs among the 11 parents. Based on the concept that these parental SNP should be reflected in RILs with different responses to FOV4, SNP indexes were further calculated between two groups (resistant vs. susceptible) of the selected RILs, as shown in Table 4. An average SNP index within each QTL region is shown in Table 5. Only those genes with a SNP index of 0.3 or higher were considered as candidates for the QTL. As a result, 48 genes met all the criteria, most of which were located on D02 and 15 genes were related to defenses to fungi (Supplementary Table 2). Interestingly, the major FOV4-resistance QTL locus qFOV4-D02-3 (with peak LOD scores greater than the FDR of 1.056 × 10 −7 at the 0.92-0.98 Mb region) contained 14 candidate genes, three of which (Gh_D02G0105, Gh_D02G0115, and Gh_ D02G0119) responded to fungal infections (Supplementary Table 2). More studies are needed to further narrow candidate genes for the QTL (see Discussion).

Discussion
FW caused by FOV is one of the most important soilborne fungal diseases affecting world cotton production. Numerous studies have investigated the genetic basis of FW resistance using bi-parental or association mapping populations, resulting in more than 100 QTLs reported, as summarized in this study. However, no multi-parent populations were developed and used for identification of genes or QTLs for resistance to FW until now. In this study, we evaluated a MAGIC population of 550 RILs derived from intermating of 11 American Upland cotton parents representing major US seed companies and public cotton breeding programs. These 11 cultivars were developed under non-FOV4 infested field conditions in a span of 24 years between 1979 and 2002. Based on a pedigree analysis, seven parents can trace their parentages to Deltapine 15 and/or Stoneville 7, while three parental lines had Acala or New Mexico germplasm in their parentages. However, Tamcot Pyramid was developed in the MAR program and its parentage could not be traced to any specific germplasm lines. It appeared that the genetic difference in FOV4 resistance detected in this study was not due to the inclusion of the backbone germplasm lines (Deltapine 15 and Stoneville 7) in their pedigrees. For example, the three more resistant lines-M-240RNR, Stoneville 474, and SureGrow 747 and the four more susceptible lines-Coker Through artificial inoculation in two replicated tests in the greenhouse, significant genotypic variation among the parents and within the MAGIC RIL population was detected, allowing estimation of broad-sense heritability (0.62). The result was consistent with our previous greenhouse studies (Abdelraheem et al. 2020a;Zhang et al. 2020b). Based on almost half a million polymorphic SNPs, 22 QTLs (7 on the A-subgenome and 15 on the D-subgenome) including one major QTL on D02 for FOV4 resistance were detected on 12 of the 26 pairs of Upland cotton chromosomes. The narrow chromosomal region at an 0.2 Mb interval (at 0.79-0.98 Mb with the peak at 0.92-0.98 Mb) for the major resistance QTL on D02 allowed identification of 14 possible candidate genes for FOV4 resistance. This study represents the first report using a MAGIC population with currently the largest number of RILs for mapping QTLs conferring resistance to FW in cotton. The results will be useful to both markerassisted selection and further studies towards cloning of candidate genes for FW resistance QTLs. Furthermore, the results from this current study allow discussion of several important issues facing genomic studies of resistance to FW in the cotton community.

Quantitative or qualitative resistance to FOV4 in cotton
FW caused by FOV4 under field conditions is an early season disease with high seedling mortality (Zhang et al. 2020b). In this study, no parents including the most susceptible ones achieved high levels of mortality and DSR because of the high temperature conditions in the greenhouse when the parents and MAGIC RILs were evaluated for FOV4 resistance. This is consistent with Zhang et al. (2020b) that very low seedling mortality was observed in more than 2,000 lines evaluated under high temperature ranges (26-32 °C) in the greenhouse. However, high mortality including 100% mortality (i.e., with DSR of 5) in 50% of the 1,200 germplasm lines evaluated were observed under low temperature (20-21 °C) settings (Zhang et al. 2020b). Responses of the 11 parents to FOV4 infections under the low temperature conditions may be different, and they may have high seedling death rates. Zhang et al. (2020b) suggested that the low temperature conditions allow detection of high levels of qualitative FOV4 resistance, while the high temperature conditions detect low levels of quantitative resistance to FOV4. Zhang et al. (2021a, b) confirmed that FW caused by FOV4 is a temperature dependent disease with an optimal temperature of 22-23 °C to cause mortality. Equally importantly, FW caused by FOV4 is also inoculum density dependent. Hao et al. (2009) showed that FW symptoms and plant growth reductions began at the inoculum density of 10 3 FOV4 conidia g− 1 of potting soil mix in highly susceptible Pima cultivar DP 744 and moderately susceptible Upland Acala Ultima. However, they increased more rapidly in the former and reached the plateau at the inoculum level of 10 5 FOV4 conidia g− 1 of soil, whereas the increase in disease severity and decrease in plant growth did not reach plateau in the latter even at the inoculum level of 10 6 FOV4 conidia g− 1 of soil. Therefore, optimal disease development conditions play an important role in determining if a genotype possesses a qualitative resistance to FOV4.
Qualitative or quantitative resistance is genotype dependent. In many studies on resistance to FW in cotton, evaluations of plant responses to FOV infections were conducted under natural infested field conditions. In other cases, artificial inoculation and evaluation may not be under optimal environmental conditions for FW disease development. Furthermore, low FOV inoculum density and non-uniformity of FOV inoculum in the soil under field conditions did not allow reliable separation of true resistance from escapes and differentiation of quantitative differences in resistance on an individual plant basis. As such, resistance to FW is often treated as a quantitative trait with low to moderate heritabilities estimated (Zhang et al. 2015a, b, c). An example is the Chinese Upland resistant cultivars to FOV7. When high and uniform inoculum density of FOV7 spores was provided, heritabilities for resistance were > 0.9 with only one resistance gene estimated based on a diallel cross mating design (Zhang et al. 1994;Feng et al. 1998). When the same inoculum method was used to evaluate F 2 and backcross generations, a typical Mendelian segregation ratio (3:1 or 1:1) was observed. Therefore, for genotypes with high and qualitative resistance, identification of a major resistance gene is dependent on both high and uniform inoculum density under optimal environmental conditions (low temperatures for FOV4) for FW disease development.

Subgenome bias of QTL distribution for disease resistance
Through a meta-analysis of reported 133 QTLs for FW resistance in cotton, this study showed a significant D-subgenome bias of FW disease resistance QTLs, which was consistent with the findings of resistance QTLs for bacterial 1 3 blight (BB, caused by Xanthomonas citri pv. malvacearum (Wright et al. 1998) and Verticillium wilt (VW, Zhang et al. 2015b;Abdelraheem et al. 2017). Wright et al. (1998) mapped 5 (83%) of the 6 resistance QTLs for BB resistance to the D subgenome. In a meta-analysis, Abdelraheem et al. (2017) found that significantly more reported QTLs (118) for VW resistance were mapped to the D subgenome than the A subgenome (83). However, the reverse was true for resistance to RKN and reniform nematodes (Rotylenchulus reniformis) with 55 and 30 QTLs on the A and D subgenomes, respectively. Therefore, whether the D subgenome evolved new R genes or resistance QTLs more rapidly than the A subgenome or vice versa depended on pathogen after the two ancestral diploid genomes were merged into the same nucleus during the origin of the tetraploid cotton.

Comparison of major QTLs mapped for FOV resistance
In this study, we mapped a major FOV4-resistance QTL close to the telomere on one chromosome arm of D2. This QTL region (at 0.7-1.2 Mb) was close to the RKN-resistant Mi2 gene (Gh_D02G0276 at 3.45 Mb) (Wubben et al. 2019). Therefore, it appeared that the FOV4-resistance QTL and the RKN-resistance QTL locus were located toward the telomere of the chromosome and closely linked (ca. 2.5 Mb apart). Recently, a major effect QTL for FOV4 resistance in NemX was also mapped to the same chromosome (D02) but on the other terminal (at 69.48-69.58 Mb) anchored by two SSR markers (BNL3932 and DPL0473 with 1.2 cM and 0.1 Mb apart) and another SSR marker MUSB0850 within a 30.7 cM region using a RIL population of 90 lines ). However, due to the limited number of SSR markers and RILs used, the mapping resolution was too rough to delineate the QTL region. This QTL was in a QTL cluster with another major-effect QTL for FOV4 resistance within the 58-63 Mb region identified by Abdelraheem et al. (2020a) using 376 US Upland cotton accessions and 25,677 SNP markers. However, this current study did not detect any FOV4 resistance QTL within this chromosome region but detected two QTL clusters at the terminal region (0.25-2.14 Mb) and another region (18.04-23.25 Mb). Previously, a major FOV4-resistance gene in Pima S-6 was mapped to D02 by Ulloa et al. (2013), but the DNA sequence for the anchoring SSR marker BNL0834 was contained in D03 (c17) at 8.2 Mb. Therefore, the exact chromosomal location of this major FOV4 resistance gene is currently unknown and should be further investigated. Interestingly, Liu et al. (2021) have recently reported the identification of Gh_D03G0209 (named GhGLR4.8) on D03 (at 1.97-2.37 Mb) encoding a GLUTAMATE RECEPTOR-LIKE (GLR) protein as the candidate gene for the resistance gene Fov7 with resistance to FOV race 7 through a GWAS using 290 Chinese Upland accessions. CRISPR/Cas9-mediated knockout of the Fov7 gene confirmed that a SNP (C/A) in the gene, leading to an amino acid change (L/I) in the enzyme, was responsible for FOV7 resistance. However, neither Gh_D03G0209, nor its A-subgenome homeolog Gh_A02G1505 was located in the vicinity of the QTLs reported in this study. This small gene family did not appear to have any other obvious members. The Fov7 gene was very close to another major race 7-resistance QTL (at 1.4-1.6 Mb) identified in 336 Chinese G. barbadense accessions by Zhao et al. (2021), both of which were in a close proximity of a FOV4-resistance QTL (at 2.86-2.89 Mb) identified by Abdelraheem et al. (2020a) in Upland accessions and a race 7-resistance QTL (at 0.99-3.07 Mb) identified in a G. barbadense RIL population by Han et al. (2022). In addition, two QTLs detected on A05 (at 90.43-90.75 Mb) and D13  in this study were also in a close proximity of two race 7 resistance QTLs reported by Han et al. (2022) Han et al. (2022), but far away from the two QTLs (at 51.1-54.2 Mb and 55.7-60.4 Mb) reported by Abdelraheem et al. (2020a).
Based on GWAS using 215 Chinese G. arboreum accessions with 1.4 million resequencing-based SNP markers, Du et al. (2018) identified the candidate gene Ga11G2353 encoding for a glutathione S-transferase (GaGSTF9) for a major-effect QTL on chromosome A11 (within 1.0 Mb interval in the 102.5-103.5 Mb region) for FOV7 resistance. Its homologous gene in Upland cotton is Gh_A11G1459; however, no QTL on A11 for FOV4 resistance was identified in our current study and for other FOV races in previous studies. In addition, similar to previous studies, we detected fewer QTLs for FOV resistance on the A subgenome than on the D subgenome. Therefore, the genetic bases (in terms of QTL numbers and locations) for FOV resistance between the diploid A genome and the tetraploid AD genome appeared different. However, interestingly, the D-subgenome homeolog, Gh_D11G1615 is located within 1-Mb of the QTL on D11 we named qFOV4-D11-1. Three other genes with significant homology to Ga11G2353 (Gh_ A11G1459) were also observed near the presently described QTLs: Gh_D11G1496 at qFOV4-D11-1, Gh_D13G1398 at qFOV4-D13-1 and Gh_D01G1754 at qFOV4-D01-1 and are worthy of further investigations.

Difficulty in delineating major QTLs for FW resistance and solutions
It appeared that neither Ga11G2353 nor Gh_D03G0209 belongs to typical plant resistance genes (i.e., R-genes) coding for molecules essential for pathogen defense including pattern recognition receptors (PRRs), wall associated kinase (WAKs), receptors with nucleotide-binding domain (NLRs) and leucine-rich repeats (LRRs). It should be recognized that the resolution power from 200 to 300 accessions or genotypes in a genetic population is usually not high enough for high resolution mapping to pinpoint a candidate gene, because a 0.3-0.5 cM chromosome region may still contain 5-10 or more genes, regardless of the number of SNP markers. Based on almost half a million SNPs and 550 MAGIC lines, this current study identified a major resistance QTL within the 0.72-0.98 Mb region which still contains 14 genes. Furthermore, it is often true that numerous germplasm accessions used were not homozygous as evidenced from various percentages of heterogeneous and heterozygous SNP markers and segregation in phenotypes including FOV resistance within the same accessions (Zhang et al. 2020b). The lack of purity of germplasm lines used and high experimental errors in phenotyping cotton for FOV resistance further compromised the accuracy and precision in locating resistance QTLs and subsequent identification of candidate genes. Therefore, to reliably map and clone a FOV-resistance QTL or gene, homozygous lines (accessions) from a large genetic population (1,000-2,000 genotypes) should be developed and employed in multiple replicated tests; and a reliable artificial inoculation and screening method for FOV resistance including optimal temperature conditions and inoculum density should be established and used. Based on our extensive work using artificial inoculations in numerous replicated tests (Zhang et al. 2020a(Zhang et al. , b, 2021bAbdelraheem et al. 2020aAbdelraheem et al. , 2021, coefficients of variation ranged between 10 and 30% in evaluating cotton for resistance to the two soilborne fungal diseases-VW and FW. A reliable experimental design with more plants line −1 in each replication and more replications will certainly reduce experimental errors (because standard error = standard deviation or mean squared error/n, here n is the number of replications or plants). It is important to use blocking control when growing and managing plants on a replication basis to ensure that plants within each replication are exposed to similar environmental conditions (e.g., temperature, light, soil type, fertilizer, and moisture). It is also equally important to ensure that the same volume of the same inoculum density is used in each plant or pot using the same inoculation method, and other soilborne pathogens should be avoided.
Author contribution statement JZ conceived and designed the study. YZ, GT, and JZ. drafted the manuscript. YZ, AA and ZT. performed the phenotyping experiment and analyzed the data. GT and DF. developed the SNP markers. GT performed the GWAS and identified QTLs and candidate genes. JJ and JM developed the MAGIC population. TW and HH supported the study. All authors read and approved the final version of the manuscript.
Funding This research was supported in part by Cotton Incorporated and New Mexico Agricultural Experiment Station.
Data availability Except for the supplementary tables, no other data is available.

Conflict of interest
The authors declare that they have no conflict of interest.
Ethical approval Not applicable.