A total of 8 tobacco accessions were used as parents of the MAGIC population. The construction of the MAGIC population used in this study is shown in Figure 1A . Beinhart1000-1 (BH) and Florida301 are cigar tobacco. Basma and Samsun are oriental tobacco. Xiaohuaqing (XHQ) and Tangpeng (TP) are sun-cured tobacco from China. Vam belongs to burley tobacco. Honghuadajinyuan (HD) is a flue-cured tobacco cultivar. The first stage was pairwise crossing of eight parents, and 4 two-way crosses were produced. The resulting 4 two-way crosses were incomplete diallel crossed to produce 6 four-way crosses. The 6 crosses were then intercrossed to produce 132 eight-way crosses. To increase the number of recombination events, the resulting eight-way crosses were intercrossed to produce 218 eight-way crosses, and the 218 crosses were randomly mated and produced 320 eight-way crosses. Three single plants of each 320 eight-way crosses were selected and propagated by single seed descent (SSD) method during five selfing generations. Finally, more than 800 eight-way tobacco MAGIC homozygous lines were constructed.

Five hundred and sixty MAGIC homozygous lines (MLs) and 8 parents were planted in Zhucheng, Shandong, China (119.40_E, 36.99_N, 61 m altitude), in 2019 and 2020 (ZC2019 and ZC2020) and in Chenzhou, Hunan, China (113.01_E, 25.78_N, 190 m altitude), in 2020 (CZ2020) with a randomized complete block design with two replications. In all trials, one replication of each line contained two 10-plant rows with a row length of 10 m, row spacing of 1.2 m and plant distance of 0.5 m. Two weeks after topping, two mature middle leaves from each plant were harvested and mixed as one sample. Three plants of each line were randomly sampled for nicotine content measurement, and the average value of three plants was used for further analysis. Nicotine content was measured using a Continuous Flow Analyzer (SAN++, Skalar Analytical B.V., Breda, the Netherlands) from ground cured leaves. Variance analysis was performed using QTL Icimapping version 4.2.53.

A total of 560 MAGIC lines and 8 parents was genotyped using the 430K tobacco SNP array with 432,362 markers at Zhengzhou Tobacco Research Institute of CNTC. Details of the genotyping and SNP calling procedures are described by Zhang et al. (2017). SNPs with no polymorphism among 8 parents, missing rate > 10% and minor allele frequency (MAF) < 0.05 were removed. Finally, a total of 93,766 SNP markers were used for subsequent studies ( Table S1 ).

The genome-wide LD between pairs of loci was performed using Tassel 5.0 software (Bradbury et al., 2007). The estimates of the LD were measured using the squared allele-frequency correlations (r²) for pairs of loci. The distance in base pairs that loci could be expected to be in LD or LD decay was computed by plotting r² onto physical distance using the threshold r^{2 =} 0.33 (one-third of the maximum value) as a cutoff. The LD contour plot was generated from PopLDdecay software (Zhang et al., 2019).

All marker-trait associations were performed using TASSEL 5.0 software. A mixed linear model (MLM) was employed for GWAS analysis by incorporating a kinship matrix along with PCA. SNPs with minor allele frequency (MAF) <0.05, no polymorphisms and deletion rates higher than 10.0% were filtered out. Components of PCA were set as 5. Pedigree analysis (K) kinship method chose Centered_IBS with Max Alleles as 6. The compression level of the MLM is the optimum level, and P3D is selected for variance component estimation. Manhattan and Q-Q plots were derived from the GWAS results using the CMplot package (https://github.com/YinLiLin/R-CMplot) within R software (http://www.r-project.org/). For results with inflation, genomic control was applied to further control the confounding with population structure (Yurii et al., 2007). The Bonferroni correction method based on the FWER (Family-Wise Error Rate) standard was used to assess the significance threshold of association sites, and the significance level was determined as p = 1.0E-07 (0.01/93,766). When multiple loci were associated with a trait and were within a 3 Mb interval, they were considered as a single QTL.

Four MAGIC parent cultivars (BH, HD, Florida301, Samsun) and one flue-cured tobacco cultivar (XHJ) were used for resequencing. Genomic DNA of the 5 tobacco accessions was randomly fragmented to 350 bp for library construction. The constructed library was sequenced on the Illumina NovaSeq 6000 platform by Novogene, Beijing, China, and 150 bp paired-end reads were generated. Furthermore, quality-controlled data were generated after removing reads with ≥10% unidentified nucleotides (N); >10 nt aligned to the adapter, allowing ≤10% mismatches; >50% bases having Phred quality <5 and putative PCR duplicates generated in the library construction process. After quality control of the raw data, the clean reads were mapped to the reference genome.

The candidate regions (83 Mb - 86 Mb) of chromosome 7 from five sequenced cultivars were compared, and the average coverage depth was counted using Tablet software (James Hutton Institute) with a window size of 100 kb.

Total DNA was extracted from fresh tobacco leaves using a CTAB method. The quality and quantity of total DNA was determined using a NanoDrop 2000 spectrophotometer (NanoDrop Technologies Inc., Wilmington, DE, USA). PCRs were performed with 2× Phanta Flash Master Mix (Dye Plus) and gene-specific primers according to the manufacturer’s instructions. PCR products were analyzed by electrophoresis and purified using FastPure Gel DNA Extraction Mini Kit (Vazyme). The purified DNA fragment was ligated into the pEASY-T1 vector (Transgen) and then transformed into Escherichia coli Trans5a competent cells (Transgen) for sequencing. All primer information is listed in Table S2 .

CRISPR/Cas9-mediated target mutagenesis was used to knock out the candidate genes, qNIC7-1 (NtERF170, NtERF199, and NtERF890) and NtERF189, the homologous gene of NtERF199. The targeting sequences (listed in Table S2 ) and their reverse complement oligomer were chemically synthesized and annealed, and the dimers obtained were then inserted into the PHSE401 binary vector. For multigene editing, the dimers are connected with the tRNA sequence.

The constructs obtained above were individually transformed into competent Agrobacterium tumefaciens strain EHA105 cells using the freeze–thaw method. Leaf discs from 6-week-old XHJ were infected with A. tumefaciens strain EHA105 harboring the knockout vectors following the methods reported (Wang et al., 2022). The infected leaf discs were then plated onto generation medium (MS medium containing 3% (w/v) sucrose, 1 mg/L 6-benzylaminopurine, 0.15 mg/L 1-naphthaleneacetic acid, 50 mg/L cefotaxime sodium, and 8 mg/L hygromycin, pH 5.7). Hygromycin-resistant seedlings were obtained, and genomic DNA was extracted. NtERF170, NtERF189, NtERF199 and NtERF890 were amplified by PCR from genomic DNA, and the PCR products were Sanger sequenced to detect the mutations.

Single-gene (M170, M199, M890) and multigene (M-2, M-3) T2 generation knockout mutants were obtained, and homozygous screening was conducted. The DNA of each individual plant was extracted. Full-length primers were used to amplify target genes, and the PCR products were purified and ligated into pEASY-T1 vectors for sequencing. At least ten clones from each sample were sequenced. Only when all clones had identical mutations were the corresponding plants considered as homozygous mutants.

Mature middle leaves of 12-week-old plants were collected, dried in an electro-thermostatic blast oven, and ground using a tissue grinder (Jingxin Industrial Development Co., Ltd., Shanghai, China). Fifty milligrams of ground leaves were sampled for nicotine content measurement on a HITACHI Chromaster (autosampler Chromaster 5210, column oven Chromaster 5310, diode array detector Chromaster 5430, RI detector Chromaster 5450). Leaves collected from three independent plants and mixed together were considered as one replicate. Three replicates were measured for each mutant line.

Different accessions were obtained from the National Medium-term GenBank of Tobacco Germplasm Resources of China. A total of 1,187 flue-cured tobacco accessions and 114 cigar accessions were planted in the field in Jimo, Shandong, China (120°44’71’’E, 36°38’93’’N, 22 m altitude) in 2020. The nicotine content was measured as described above. Leaf samples were collected, and genomic DNA was extracted. NtERF189 and NtERF199 were amplified in 1301 DNA samples using both full-length primers and short segment primers. PCR products were analyzed by electrophoresis, and those with no target amplification were considered as target gene deletions.

Eight tobacco accessions (BH, Florida301, Vam, Basma, Samsun, XHQ, TP, HD) that belong to different tobacco types were used to produce the MAGIC population ( Figure 1A ). The effectiveness of random mating and its consequence on the structure of this population was evaluated using 93,766 SNP markers distributed on 24 chromosomes ( Figure 1B ). Population structure analysis by the neighbor-joining method indicated that no specific population structure exists among the 560 MAGIC lines ( Figure 1C ). Kinship analysis also showed no specific clustering pattern ( Supplementary Figure 1 ). To determine the mapping resolution for GWAS, the average extent of LD decay was quantified ( Figure 1D ). Using the whole set of SNPs, the LD decay rate of the population for the entire genome was estimated to be 1,500 kb, with r^{2 =} 0.33 (one-third of the maximum value). All data indicated that the tobacco MAGIC population was suitable for GWAS analysis with no obvious population structure and moderate LD level.

Nicotine content data of 8 parents and 560 MLs were collected from three environments (2019 in Zhucheng, 2020 in Zhucheng, 2020 in Chenzhou). In this study, the nicotine content of mature middle leaves was measured, and the results are shown in Figure 2 and Supplementary Table 3 . As shown in Table 1 , the mean values of parents and 560 MLs were similar. However, the range of variation was more extensive in MLs, possibly due to transgressive segregation. The variation in nicotine content was determined by genotype and environment. The heritability was approximately 0.5 by plot, demonstrating that the trait was quantitatively inherited.

GWAS was carried out using the best linear unbiased estimates (BLUEs) of MLs over three replicates in a mixed linear model (MLM) by employing Tassel 5.0 software. The quantile-quantile (Q-Q) plot showed that most of the observed p values followed a uniform distribution, but the few that were in LD with a causal polymorphism had significant p values in the tail ( Supplementary Figure 2 ). The Manhattan plot for nicotine content was generated from TASSEL software and the R package “CMplot” ( Figure 3 ). At genome wide significance (p value ≤ 1.0E-07), a total of 4 QTLs (qNIC7-1, qNIC7-2, qNIC7-3, qNIC5-1) distributed on chromosome 7 and chromosome 5 were identified ( Table 2 ). Among them, qNIC7-1 and qNIC5-1 were repeatedly detected in most environments. The other 2 QTLs (qNIC7-2, qNIC7-3) were only identified in ZC2019 and ZC2020.

The most significant locus, qNIC7-1, associated with nicotine content had 19 markers spanning 772 kb ( Supplementary Table 4 ). The SNPs with the strongest association (AX-117513847, AX-117513846, AX-117632223, AX-117522763, p = 6.85E-36) explained 32.6% of the phenotypic variation. Based on the results of LD decay analysis and association analysis, we estimated the candidate region of qNIC7-1 to be 83 Mb - 86 Mb, which contained 28 SNPs and 41 annotated genes ( Supplementary Table 5 ). In qNIC5-1, the SNP with the strongest association (AX-117702271, p = 3.01E-12) explained 9.2% of the phenotypic variation. The candidate region of qNIC5-1 was estimated to be 159.5 Mb-162.5 Mb, which contained 64 SNPs and 46 annotated genes ( Supplementary Table 5 ).

We noticed that eight tandem ERF transcription factors were located on the region covering the strongest-associated SNPs (84 Mb - 85 Mb) in qNIC7-1. Resequencing data analysis showed that three of the ERFs (NtERF170, NtERF199, and NtERF890) were deleted in one parent cultivar, Beinhart1000-1, which had the lowest nicotine content, whereas the remaining seven showed no sequence or expression differences ( Figure 4 ). Therefore, we focused on the three deleted ERFs. To examine whether their deletions are responsible for nicotine content reduction, we used CRISPR/Cas9 technology to generate loss-of-function mutants of three genes separately (M170, M199, M890) and collectively (M-3) in the flue-cured tobacco cultivar XHJ and evaluated their phenotypic effect in the T2 generation. The genotypes of wild-type (XHJ), single-gene homozygous mutants (M170, M199, M890) and multigene homozygous mutant (M-3) plants were determined by DNA sequencing. The average nicotine content of XHJ was 3.50 mg/g. The average nicotine contents of M170 and M890 were 3.37 mg/g and 3.26 mg/g, respectively, showing no obvious differences compared to the wild type. However, M199 and M-3 exhibited significant decreases with the nicotine contents of 1.19 and 0.97 mg/g, approximately 1/3 of that in wild-type XHJ, corresponding to the nicotine level in BH (1.11 mg/g) ( Figure 5 ).

BLAST analysis showed that NtERF189 on chromosome 19 is homologous to NtERF199, sharing an identical binding domain. To investigate the effect of both genes on nicotine synthesis, we also generated a two-gene mutant (M-2) in which NtERF189 and NtERF199 were collectively knocked out ( Figure 5 ). As expected, nicotine was barely detectable in M-2 plants (0.03 mg/g), approximately 1/100 of that in WT plants, 1/30 of BH, M199 and M-3.

The candidate region of qNIC5-1 contained 46 annotated genes. According to the annotation and gene function analysis, we predicted two candidate genes (NIC5-1 and NIC5-2) that might be associated with nicotine content. NIC5-1 was annotated as a MYB35 transcription factor and located 304 kb downstream of AX-117702271. NIC5-2 was located 1,188 kb downstream of AX-117702271 and homologous to the MYC2 transcription factor. The functions of the two genes in nicotine synthesis need to be further studied in the future.

To explore the practical utility of the deletion mutation in qNIC7-1, we designed an insertion-deletion (InDel) marker (NIC7001) (forward primer: 5’-GATGTGGGTGTCAACCTTTTCG-3’ and reverse primer: 5’-CGAACCCTTCCTTTCACATAA-3’) and genotyped 560 MLs and 1301 tobacco accessions. The InDel marker NIC7001 locus has two types: deletion type and nondeletion type (839 bp amplification). The two types were compared using one-way ANOVA using IBM SPSS Statistics software version 23 (IBM, New York, USA). In the MAGIC population, the lines of the nondeletion type had a mean nicotine content of 3.95 mg/g, which was significantly higher than the mean nicotine content of the deletion type lines (1.77 mg/g). In 1301 tobacco cultivars, which included 1187 flue-cured tobacco accessions and 114 cigar accessions, only 11 cultivars (8 flue-cured tobacco types and 3 cigar tobacco types) were detected as deletion types ( Table 3 ). Five out of the 11 cultivars, together with six nondeletion-type cultivars, were phenotyped. In accordance with the results in MAGIC lines, the nicotine content of deletion type cultivars was significantly lower than that of nondeletion type cultivars ( Figure 6 ).

The distribution of NtERF189 (NIC2) deletion was also investigated using the InDel marker NIC7002 (forward primer: 5’-GCGTGCAATAGGGCAAAGC-3’ and reverse primer: 5’-TGAAAGCTTCCTTCCTTTCACA-3’) in 1301 tobacco germplasm cultivars ( Table 3 ). The deletion of NtERF189 was detected in 14 flue-cured tobacco accessions (1.18%) and 15 cigar accessions (13.16%), with a total deletion rate of 2.23%. NtERF199 was deleted in only 8 flue-cured tobacco accessions (0.67%) and 3 cigar accessions (2.63%), with a total deletion rate of 0.85%. No accession was identified that carried deletions in both genes.

To identify genetic variation controlling complex traits variation, many biparental populations have been widely used for QTL mapping in tobacco (Drake-Stowe et al., 2017; Sun et al., 2018; Agacka-Mołdoch et al., 2021). However, the main limitations of biparental population for QTL mapping are the low genetic diversity of the mapping population and mapping resolution. In MAGIC population, multiple founders are intercrossed several times in a well-defined order to combine the genetic material of all the founders in a single line (Cavanagh et al., 2008). Therefore, compared to biparental mapping population, MAGIC population contains more recombination events and genetic diversity, which enhances mapping resolution and the ability to analyze several alleles simultaneously. In this study, eight founders in the MAGIC population enriched population with higher allelic diversity compared to those derived from the traditional biparental populations. Furthermore, we narrowed down the target region, and NtERF199 was considered exactly as the candidate gene by MAGIC population. In addition, the distribution of NtERF199 in tobacco germplasm resources was analyzed. The deletion of NtERF199 was considered as a rare allelic variation (RAV) with a deletion rate of 0.85%. Sun et al. (2017) performed GWAS using 219 flue-cured tobacco accessions and identified 2 QTLs associated with nicotine content by GWAS. However, the most major-effect QTL (qNIC7-1) in this study was not identified using the 219 flue-cured tobacco accessions, probably due to the filteration of RAV. Because the crossover technique elevates the frequency of all parental alleles to be almost equal, MAGIC population have higher power to identify low frequency alleles and can better evaluate allelic effects between founders than genome-wide association panels (Dell'Acqua et al., 2015; Scott et al., 2020). Another main advantage of MAGIC population over association panels is the lack of an underlying unknown structure, which enhances the risk of detecting false-positives. Therefore, MAGIC population represents an intermediate resources between biparental crosses and association panels regarding the number of traits that can be investigated, resolution, allele diversity and population structure (Rakshit et al., 2012; Pascual et al., 2016). We believe that this tobacco MAGIC population provides a good resource for studying the genetic basis of complex traits in tobacco.

GWAS is a forward genetic approach to identify underlying causal genes, mutations and putative functional markers that affect complex quantitative traits. In this study, we identified four QTLs distributed on two chromosomes associated with nicotine content. The candidate gene NtERF199 on qNIC7-1 and its homologous gene NtERF189 were cloned and experimentally verified through loss-of-function mutant analysis. The results are consistent with previously reported studies. Shoji et al. (2010) identified that the NIC2 locus comprises a cluster of ERFs (including NtERF189) that regulate nicotine biosynthesis in tobacco. Qin et al., (2021) validated NtERF199 as the NIC1 gene. The detection of qNIC7-1 highlighted the strength of using tobacco MAGIC population in GWAS.

As in the case of NIC2 locus in tobacco, eight ERFs are clustered together in the qNIC7-1 locus. However, both the previous study by Qin et al., (2021) and our investigation confirmed that NtERF199 is the only causal gene for nicotine biosynthesis within the qNIC7-1 locus. Same case also exists in tomatoes, with GAME9 locates within an ERF cluster and works as the only functional regulator of steroidal glycoalkaloid biosynthesis (Cardenas et al., 2016). These results suggest potential functional redundancy and divergence of clustered ERFs.

Allopolyploid plants usually display higher adaptability than their diploids ancestors since their preexisting homeologous loci from the progenitors show a predominantly additive gene expression pattern (Fu et al., 2016). As a natural allopolyploid plant, N. tabacum was formed by hybridization of Nicotiana sylvestris (S-genome) and Nicotiana tomentosiformis (T-genome) (Clarkson et al., 2017). Blast analysis and functional studies have demonstrated that NtERF199 and NtERF189 are a pair of homologous genes originated from N. sylvestris and N. tomentosiformis, respectively, and showed dose-dependent effects on nicotine levels (Shoji et al., 2010; Qin et al., 2021). It has also been proven through transcriptome analysis that most homologous genes from the S- and T-genomes of tobacco have been maintained and are expressed in a similar way (Bombarely et al., 2012; Edwards et al., 2017). Therefore, the effect of homologous genes must be taken into consideration in the study of gene function in allotetraploid plants, especially N. tabacum.

We also identified a significant QTL on chromosome 5 (qNIC5-1), which explained 9.6% of the phenotypic variation. Among the 46 annotated genes in qNIC5-1, we predicted two candidate genes (NIC5-1 and NIC5-2) likely associated with nicotine content based on the annotation and gene function analysis. NIC5-1 was annotated as a MYB35 transcription factor. MYB35 has been reported to regulate methyl jasmonate (MeJA) and wound-responsive G10H-1 in Bacopa monnieri (Jeena et al., 2021). Nicotine biosynthesis in tobacco is highly regulated by jasmonic acid (JA) (Baldwin et al., 1996; Shoji et al., 2008). Therefore, MYB35 might participate in nicotine synthesis by regulating JA-responsive and nicotine biosynthesis-associated genes. NIC5-2 was annotated as a MYC2 transcription factor, which has been described to control nicotine biosynthesis genes by directly binding the G-box in the target promoters and upregulating the NIC2-locus ERF genes (Shoji and Hashimoto, 2011; Zhang et al., 2012; Sui et al., 2021).

To investigate the distribution of NtERF189 and NtERF199 deletion, we manually designed two InDel markers (NIC7001 and NIC7002) and screened a total of 1301 cultivated tobacco accessions, including 1187 flue-cured tobacco accessions and 114 cigar tobacco accessions. In flue-cured tobacco, the deletion mutation rates of NtERF189 and NtERF199 were 1.18% and 0.67%, respectively, while in cigar tobacco, they reached 13.16% and 2.63%, respectively. The results demonstrated that the two ERF deletion mutations displayed higher frequency in cigar tobacco. In fact, the low-nicotine trait was initially found in strains of Cuban cigar tobacco varieties and subsequently introduced into cigarette varieties through a series of backcrosses (Valleau, 1949). This might explain the higher mutation frequency in cigar tobacco. Overall, the deletion mutation of NtERF189 and NtERF199 was detected in only 2.23% and 0.85% of the tested tobacco cultivar, demonstrating them as rare allelic variations. Deletion of both genes was detected in no accessions. This may be attributed to the importance of nicotine in resisting biological stresses, and those with both ERFs mutated have already been naturally eliminated.

Moreover, we validated the effectiveness of NIC7001 using a panel of 11 tobacco cultivars (5 deletion types and 6 nondeletion types). The results provide promising evidence of InDel marker application for future breeding programs of low nicotine tobacco cultivars.