The RILs were generated from two elite flue-cured tobacco parents Y3 and K326. Y3 is a backbone cultivated variety that originated from Zimbabwe with elite agronomic traits and complicated parental sources. K326, whose genome has been assembled (Edwards et al., 2017), was introduced from America with high commercial quality and disease resistance but moderate agronomic performance. A total of 274 RILs were employed in this study, consisting of two parents, one F₁ generation hybrid (YKF₁; Y3   ×   K326) and 271 F₇ generation individuals. The materials were planted at Shilin (N: 23.46; E: 103.17) field experiment stations using complete random design with 5 replications, and were cultivated according to local technical measures for quality tobacco production. The eight hazardous substances, including benzo[α]pyrene (B[a]P), carbon monoxide (CO)¸ crotonaldehyde (CRO), hydrogen cyanide (HCN), ammonia (NH₃), 4-(methylnitrosamino)-1-(3-pyridyl)-1-butanone (NNK), phenol (PHE) and tar (TAR) and were collected in the mainstream smoke of cigarettes produced using tobacco planted at Shilin in 2018, 2019 and 2020; and calculated hazard index (HI). Three combinations of location and year were treated as environments denoted as E1 (2018 Shilin), E2 (2019 Shilin) and E3 (2020 Shilin).

B[a]P was analyzed using the gas chromatography mass spectrometry (GC-MS) method, as described by the Chinese standard method GB/T 21130-2007 (Standardization Administration of China, 2007). CO was determined in the vapor phase using a non-dispersive infrared analyzer, as described by the Chinese standard method GB/T 23356-2009 (Standardization Administration of China, 2009). CRO was analyzed using the high-performance liquid chromatography (HPLC) method, as described by the Chinese tobacco industry standard method YC/T 254–2008 (Standardization Administration of China, 2008a). HCN was quantified using Ion Chromatography, the method as described by the Chinese tobacco industry standard method YC/T 403–2011 (Standardization Administration of China, 2011). NH₃ was analyzed using Ion Chromatography, the method as described by the Chinese tobacco industry standard method YC/T 377–2017 (Standardization Administration of China, 2017). NNK analysis was carried out using GC-TEA method, as described by the Chinese standard method GB/T 23228-2008 (Standardization Administration of China, 2008b). PHE was analyzed using HPLC, the method as described by the Chinese tobacco industry standard method YC/T 255–2008 (Standardization Administration of China, 2008c). TAR was analyzed using a routine analytical smoking machine, as described by the Chinese standard method GB/T 19609-2004 (Standardization Administration of China, 2004). HI was calculated by the following formula, H I = ( X C O 14.2 + X H C N 146.3 + X N N K 5.5 + X N H 3 8.1 + X B [ a ] P 10.9 + X P H E 17.4 + X C R O 18.6 ) × 10 / 7 (Xie et al., 2009), where X C O represents emission level of CO (mg/cigarette), X H C N for HCN (μg/cigarette), X N N K for NNK (ng/cigarette), X N H 3 for NH₃ (μg/cigarette), X B [ a ] P for B[a]P (ng/cigarette), X P H E for PHE (μg/cigarette) and X C R O for CRO (μg/cigarette). Three biological replicates were used for each assessment.

Variance components analysis and heritability estimation were performed based on the following linear model,

where y k h i is the phenotypic value of the i-th replication of the k -th line in the h -th environment; μ is the population mean; g k is the genotypic value of the k -th genotype, random effect, ɡ k ∼ N ( 0 , σ ɡ 2 ) ; e h is the effect of the h -th environment, random, e h ∼ N ( 0 , σ e 2 ) ; ɡ e k h is the interaction effect between the k -th genotype and the h -th environment, random, ɡ e k h ∼ N ( 0 , σ ɡ e 2 ) ; ε k h i is the residual effect of the individual, random, ε k h i ∼ N ( 0 , σ ε 2 ) .

The mmer module of sommer R package (Covarrubias-Pazaran, 2016) was applied to estimate the variances of random effects ( σ ^ ɡ 2 , σ ^ e 2 ,   σ ^ ɡ e 2 , σ ^ ε 2 ) and to predict the random effects by BLUPs (best linear unbiased predictions, g ^ k , e ^ h , g e ^ k h ) by solving the mixed model equation (MME). Heritability was estimated with the formula h g 2 = σ ^ ɡ 2 / ( σ ^ ɡ 2 + σ ^ ɡ e 2 + σ ^ ε 2 ) and h g e 2 = σ g e 2 / ( σ ɡ 2 +   σ ɡ e 2 + σ ε 2 ) , where σ ^ ɡ 2 was the estimated genotypic variance, σ ^ ɡ e 2 was the estimated variance due to gene-by-environment interaction, and   σ ^ ε 2 was the estimated residual variance. The rcorr module of Hmisc R package (https://cran.r-project.org/web/packages/Hmisc/index.html) was employed to calculate Pearson correlation coefficients between six studied traits: (1) phenotypic correlation coefficients with y k h i for each environment, respectively; (2) genetic correlation coefficients with y ^ k . ( y ^ k = μ ^ + g ^ k , where y ^ k is the adjusted genotypic value of the k -th line by the environment effects, μ ^ is the estimated population mean, and g ^ k is the genotypic value of the k -th line predicted by BLUP).

The high-resolution linkage map (Tong et al., 2023) that contained 46,324 markers, which were classified into 7,107 bins, a group of markers with least genotype missing rate and same genetic distance, distributed on the 24 linkage groups (LGs) and covered 3334.88 cM with an average genetic distance of 0.469 cM, was employed in this study.

A QTL full model was adopted for modeling the genetic architecture of complex traits from multi-environment trials, which includes additive effect (a) of each QTL, additive-by-additive epistatic effect (aa) of each pair of epistatic QTL, treated as fixed effects, and their corresponding environment interaction effects (ae and aae) as random effects. Suppose a trait is controlled by s segregating QTLs, of which t pairs of QTLs are involved in epistasis. Then, the phenotypic value of the m-th replication of the k -th genotype in the h -th environment ( y h k m ) can be expressed by the following mixed linear model (Yang et al., 2007),

where, μ is the population mean; a i is the additive effect of the i -th QTL with indicator variable x i k , fixed effect; a a i j is the additive-by-additive epistatic effect of the i -th QTL and the j -th QTL with indicator variable x i k x j k , fixed effect; e h is the effect of the h -th environment, random effect, e h ∼ ( 0 , σ E 2 ) ; a e h i is the additive-by-environment interaction effect of the i -th QTL and the h -th environment with observation u h i k ( = x i k ), random effect, a e h i ∼ ( 0 ,   σ A i E 2 ) ; a a e h i j is the interaction effect of the a a i j and the h -th environment with observation u h i j k ( = x i k x j k ), random effect, a a e h i j ∼ ( 0 , σ A A i j E 2 ) ; and ϵ h k m is the residual effect of the individual, random, ϵ h k m ∼ ( 0 ,   σ ϵ 2 ) .

QTLNetwork 2.0 software were employed to detect QTLs by the mixed-linear-model-based composite interval mapping (MCIM) method (Yang et al., 2008). One- and two-dimensional genome scans for QTLs were performed with configurations of 10 cM testing window, 1 cM walking step and 10 cM filtration window size. To control the experiment-wise type I error rate, a critical F-value based on the Henderson III method was determined by the permutation test with 1,000 times for each tested locus at the significance level of 0.05. Based on the significant QTL, a QTL full model was established and used to estimate each parameter based on the samples generated by Markov chain Monte Carlo (MCMC) with 20,000 Gibbs sampler iterations.

The physical positions of the marker interval with QTL were determined using Nucleotide BLAST module of NCBI (https://blast.ncbi.nlm.nih.gov/Blast.cgi), which utilized sequence information of two adjacent bin markers in the linkage map. Variants including SNPs and Indels located within the QTL regions were selected for subsequent filtration. These variants were annotated by software SnpEff (http://pcingola.github.io/SnpEff/) based on K326 reference genome (https://solgenomics.net/ftp/genomes/Nicotiana_tabacum/edwards_et_al_2017/assembly/Nitab-v4.5_genome_Chr_Edwards2017.fasta ), and those showing HIGH or MODERATE impact on related protein effectiveness were retained, according to annotation results. Then, the eligible variants with P value less than 0.05 were identified by performing single-marker regressions with y ^ k as response variable (as mentioned in the Statistical analysis of phenotypes). Before performing the enrichment analysis of genes with variants above, the protein sequences of the K326 reference genome (https://solgenomics.net/ftp/genomes/Nicotiana_tabacum/edwards_et_al_2017/annotation/Nitab-v4.5_proteins_Edwards2017.fasta) should be uploaded to the eggNOG-mapper website (http://eggnog-mapper.embl.de/) for functional annotation. Gene Ontology (GO) and KEGG pathway enrichment analyses were carried out using the clusterProfiler R package (https://github.com/YuLab-SMU/clusterProfiler) as criteria for predicting candidate genes.

For the nine smoke-related traits, the estimated heritability of genetic effects h ɡ 2 ranged from 11.03% for CRO to 20.00% for CO, displaying a limited stability across environments for these traits ( Table 1 ). Furthermore, gene-by-environment effects were found to contribute significantly more to phenotypic variation with most estimated interaction heritability ( h ^ ɡ e 2 ) nearly over 80%, indicating selection of these traits should design specific strategy for different environment conditions. Most phenotypic correlations of the traits were positive and reached statistical significance ( α = 0.05 ), for example, with coefficients over 0.8 between CO and CRO in three environments ( Figures 1A-C ). Nevertheless, strong negative correlations were observed between NNK and other three constituents: B[a]P, CO, and CRO, while, weak negative correlations between B[a]P and PHE, and between NNK and TAR. In general, the phenotypic correlations in three environments showed similar pattern as well as the genetic correlations ( Figure 1D ), indicating a strong and reliable underlying genetic relationship among these traits that make further exploration necessary. However, inconsistencies in the performance of the traits across three environments could be observed, with E2 exhibiting a higher frequency of exceptions. For instance, the phenotypic correlations between CO and PHE were significantly positive in E1 ( Figure 1A ) and E3 ( Figure 1C ), but significantly negative in E2 ( Figure 1B ). The genetic correlation coefficients, calculated using the estimated genotypic values, exhibited the highest correlation of 0.84 between CO and CRO, followed by 0.79 between NH₃ and HI, and several other strong correlations above 0.7, such as NH₃ and PHE ( Figure 1D ). The existence of high underlying genetic correlations between these traits might help to explain their strong phenotypic correlations. Further, correlations between HI and the other traits were all displayed significantly positive except B[a]P, which is consistent with the fact that HI was a composite index comprising seven representative smoke traits.

A total of 19 QTLs were identified responsible for the nine smoke-related traits. Among them, NH₃ had four QTLs; CO, NNK, PHE and TAR each had three QTLs; CRO, HCN and HI each only detected one QTL ( Table 2 ). These QTLs distributed on seven linkage groups (LGs), with LG15 containing the highest number of QTLs (7 QTLs), followed by LG6 (6 QTLs) and LG1 (2 QTLs), while each of LG4, LG5, LG10 and LG19 only harbored one QTL.

It is noteworthy that qPHE1 and qTAR1 located in the same marker interval with flanking marker PT61201 and SNP_0912209_288, similarly, qPHE6 and qTAR6 in the same marker interval ranged by SNP_0002499_215080 and SNP_0011326_16206, qNH₃15-2 and qTAR15 in the same interval ranged by SNP_0002539_133 and SNP_0000535_1549 ( Table 2 ). These co-location indicated the potential pleiotropic effect of QTL on PHE and TAR or on NH₃ and TAR, which requires further investigation for revealing the molecular mechanism of genetic correlation between traits. Our inference on existence of pleiotropic QTL was enhanced by the significant and relatively high genetic correlations estimated between the studied traits ( Figure 1 ), which were 0.74 between NH₃ and TAR, and 0.56 between PHE and TAR as we mentioned above.

Totally, 17 QTLs with additive (a) main effects were detected for nine traits, of which six QTLs also showed additive-by-environment interaction (ae) effects ( Table 3 ). Most of the QTLs exhibited small additive effects, which were regarded as minor-effect QTLs and accounted for approximately 2% phenotypic variance. The average proportion of phenotypic variance explained by individual QTL ( h a 2 ) was 2.28%. Moreover, around 70% of all QTLs each explained less than 3% of the phenotypic variance, the rest QTLs taking relatively larger effects on five traits, qHCN6 ( h a 2 =4.99%), qHI6 ( h a 2 =3.39%), qTAR6 ( h a 2 =3.14%), qCO15 ( h a 2 =3.22%), qNH₃15-2 ( h a 2 =3.21%), located mainly on LG6 and LG15. Most QTLs, whose homozygous genotypes of the alleles from the parent K326 (QQ), contributed positive additive effects (i.e., increasing the trait value) ( Table 4 ); in contrast, the other QTL genotype with alleles from the parent Y3 (qq) contributed negative additive effects.

Significant additive by environment interaction effects (ae) were found for less than half of QTLs, and their contributions to phenotypic variation ( h a e 2 ) were around 1%, with most being lower than their corresponding additive effects. However, there was an exception that the interaction between qNNK6 and two environments (E1, E2) accounted for 1.17% of the phenotypic variance, which was higher than 0.80% explained by corresponding additive effects. Moreover, the ae effects of QTL could take same or opposite effect direction as their main effects (a), and also showed different effect direction across environments. However, the ae effects of different QTLs for same trait may show consistent effect direction in same environment. For example, the ae effects of qNNK6 and qNNK15-1 both contributed negative effects in E1. It is widely recognized that QTLs with no significant ae effects have crucial application potential in breeding new varieties with strong environment adaptability. Of significance, the absence of ae effects of QTLs for HCN, NH₃, TAR, and HI suggested that these four traits may exhibit stable performance across various environments.

Two-dimensional genome scan found a digenic epistatic QTL pair (qNH₃5- qNH₃10) with epistasis-by-environment interaction effects (aae) for NH₃, but no paired epistatic QTLs for other traits ( Table 3 ). This epistatic QTLs didn’t contribute additive-additive epistatic effects (aa), only positive aae effects in E1 and negative in E2, accounting for 1.51% of total phenotypic variation. Additionally, the epistatic QTL pair qNH₃5 on LG5 and qNH₃10 on LG10, each had no individual additive effect.

Overall, environmental effect (E) explained the largest part of the phenotypic variance which ranged from 26.23% to 53.83%, followed by genetic main effect (G) and gene-by-environment interaction effect (GE) ( Supplementary Table S1 ). For all traits, the G effects, which were entirely composed of additive effects of one to three QTLs, explained 1.04% to 8.00% of the phenotypic variance. Besides, for CO, CRO, NH₃, NNK, PHE, the contribution of GE to the total variation was small, with some in the form of ae effects and others in the form of aae effects. In summary, our analysis showed that these traits were primarily controlled by a single major gene or polygenes.

To explore the potential of the identified QTLs in improving smoke-related traits through genetic and molecular manipulation, we undertook QTL genotype design and assessed the trait values potentially achieved by the general superior homozygous line (GSL(-)) and the superior homozygous line for each year (SL(-)) based on the genetic effects of these QTLs. Notably, for the HCN, TAR, and HI, there are no QTL involved in interaction with environment, the designed superior genotypes GSL(-) and SL(-) of each trait remained consistent across all environments ( Table 4 ). On the CO, CRO, NH₃, NNK, and PHE, although they were influenced by gene-environment effects, their superior lines (GSL, SL) have identical homozygous genotypes for each trait across three environments, respectively, except the SL(-)1 with QQ at qNNK6 for NNK. In particular, the superior lines of HCN, NNK and PHE were constructed by all homozygous genotypes with alleles from the parent Y3 (qq) at QTLs except qNNK6 for NNK, which indicated that the pyramid of genes from Y3 at these QTLs can effectively decrease the levels of hazardous smoke constituents, HCN, NNK, and PHE. Meanwhile, this advantageous genotype remained valid under all environments. This discovery provided valuable indication for efficient simultaneous improvements of HCN, NNK, and PHE by utilization of elite genes in Y3.

Annotation of variants (SNP/Indel) was conducted using SnpEff based on the K326 reference genome, and 559,513 variants were identified within the putative physical positions of the 17 additive-QTL regions. Among these, only variants with a HIGH or MODERATE impact on protein function were retained, as they were considered potentially functional in candidate genes, as a result, 5,268 variants were selected in total. Then simple regression analysis of the predicted genetic values on these variants was performed, and 600 significant variants in 76 genes were selected ( P < 0.05 ) ( Supplementary Table S2 ). These prioritized genes demonstrated enrichment in five GO biological processes, four GO molecular functions, and two KEGG pathways. Specifically, the enrichment was observed in cellular processes (biological process), catalytic activity (molecular function), and signaling molecules metabolic pathways ( Table 5 ). For instance, “callose deposition in cell wall”, “polysaccharide localization” and “callose localization” were associated with carbohydrate biosynthesis.

The genes Nt21g04598.1, Nt21g04600.1 and Nt21g04601.1 were predicted as potential pleiotropic candidate genes for qPHE1 (PHE) and qTAR1 (TAR), BLASTP results indicated their putative functions. Nt21g04598.1 was shown to encode a protein highly homologous to non-specific phospholipase C4 in Arabidopsis thaliana, which functioned as plasma membrane bound and promoted tolerance to phosphate deficiency and hyperosmotic stress (Nakamura et al., 2005; Peters et al., 2010; Wimalasekera et al., 2010; Kocourková et al., 2011; Pejchar et al., 2015; Yang et al., 2021). Nt21g04600.1 and Nt21g04601.1 were predicted to code a protein with high homology to transcription factor bHLH91 in Arabidopsis thaliana, which regulated the transcriptional expression, thereby regulating the plant’s adaptive responses (Qian et al., 2021). The gene Nt20g03473.1 was pinpointed as a candidate gene for qCO6 (CO), which encoded a protein highly homologous to phospholipase D (PLD) delta in Arabidopsis thaliana. PLD delta has been proposed to play a role in many cellular processes such as signal transduction, membrane trafficking, cytoskeletal rearrangements, and membrane degradation (Distéfano et al., 2012; Guo et al., 2012; Uraji et al., 2012; Jia et al., 2013); for example, it was involved in H₂O₂ and abscisic acid (ABA)-induced stomatal closure, nitric oxide (NO) signaling and ABA-promoted senescence. Moreover, gene Nt16g00273.1 and Nt16g00284.1 have been identified as candidate genes for qNNK6 (NNK), which have been predicted to encode uncharacterized serine-rich protein C215.13 and type II inositol polyphosphate 5-phosphatase 15 isoform X2 in Nicotiana tabacum, respectively, awaiting further annotation in the future.