Four plant materials were utilized in this study, including N. glutinosa (the primary source of N gene), Coker 176 (a flue-cured tobacco variety possessing the N gene), 0970A (a flue-cured tobacco inbred line possessing the N gene inherited from Coker 176), and CB-1 (a TMV-susceptible elite flue-cured tobacco variety, and one of the parents of 0970A).

To enhance the efficiency of TMV resistance identification and to avoid the impact of TMV inoculation on the growth of tobacco plants, we employed an in vitro identification method (Yin et al., 2023) in this study. Diseased tobacco leaves were collected, ground, and filtered to prepare a TMV solution. At the seven-leaf stage of the tobacco seedlings, one well-grown leaf from each plant was selected. A toothbrush was used to apply the TMV solution to the surface of the leaf, which was then placed in a self-sealing bag and incubated in a growth chamber at 26 °C with a 12-hour light and 12-hour dark cycle. After 72 hours, the leaf conditions were observed. Leaves containing the N gene (resistant to TMV) would form immune spots, while those without the N gene (susceptible to TMV) would not form such spots (Yin et al., 2023).

The genome of 0970A was sequenced and assembled by Berry Genomics (Beijing). Genomic DNA was extracted from approximately 300 mg of fresh leaf tissue using the CTAB method. DNA quality was assessed using a NanoDrop™ ND-2000 spectrophotometer (Thermo Fisher Scientific), a Qubit™ fluorometer (Invitrogen, Q33238), and the CHEF Mapper™ XA Chiller System (Bio-Rad, 1703672). A high-fidelity (HiFi) library was prepared with the SMRTbell^® Prep Kit 3.0 (Pacific Biosciences, 102-182-700) and sequenced on the PacBio Revio^® system (Pacific Biosciences, 102-090-600). The raw reads were processed with SMRT software (Ardui et al., 2018) to generate HiFi reads, which were then assembled using Hifiasm (Cheng et al., 2021). The resulting contigs were subsequently deduplicated using Purge dups (Guan et al., 2020). Contig quality and assembly completeness were assessed with Quast (Mikheenko et al., 2018) and BUSCO (Manni et al., 2021). All software was run using default parameters.

The N gene sequence (gene ID LOC107824161) was obtained from the public database GenBank. This sequence was then compared with all the contigs of 0970A through the BLAST comparison tool to identify the contig containing the N gene (termed N-contig). The collinearity between the N-contig and the CB-1 genome (unpublished) was analyzed to identify the location of the N-contig in the tobacco genome using Minimap2 software under default parameter -x asm20 (Li, 2018). Annotation of the N-contig was achieved by annotating all contigs. Initially, repetitive sequence elements within all contigs were identified using RepeatModeler v2.0.5 software (Flynn et al., 2020) and subsequently masked using RepeatMasker v4.1.2.p1–1 software (Tempel, 2012). Three strategies were then employed to predict protein-coding genes. The first strategy involved transcriptome-based prediction, utilizing the tobacco transcriptomes (NCBI accession numbers: SRR1199073, SRR1199122, SRR1199197, SRR1199202) and hisat2 v2.2.1–6 software (https://daehwankimlab.github.io/hisat2). The second strategy was homology-based prediction, leveraging the annotation sets from tomato (GCF_036512215.1), tobacco TN90 (GCA_000715135.1), and chili pepper (GCA_002878395.3). The third strategy was de novo prediction, employing GeneMark-ETP v1 software (Brüna et al., 2020) and AUGUSTUS v3.3.3 (Stanke et al., 2004). These results were integrated to generate an initial annotation file using Braker3 v3.0.8–0 software (Gabriel et al., 2024) and further refined using TSEBRA v1.2.5 software (Gabriel et al., 2021), resulting in the final protein annotation file.

The genomic DNA of CB-1, Coker 176, and N. glutinosa was extracted using the Plant Genomic DNA Kit (TIANGEN, Beijing) and the three DNA samples were subsequently sent to Berry Genomics (Beijing) for deep sequencing (at least 80 Gb each in depth) on the DNBSEQ-T7 platform. For each sample, adapter-trimmed reads were aligned to all contigs of 0970A using BWA-MEM2 v2.2.1 (Vasimuddin et al., 2019) under default parameters, with the exception of the -k option (minimum seed length), which was increased from 19 to 30 to enhance mapping stringency. The average depth of reads mapped to the N-contig was calculated every 10 kb using Samtools v1.19 (Li et al., 2009).

To introduce the N gene from 0970A into CB-1 and disrupt and eliminate the linkage drag associated with the N gene, we crossed 0970A with CB-1 and performed consecutive backcrosses up to the BC₄F₁ generation. In each generation, we selected a small number of plants that were morphologically most similar to CB-1 but retained the N gene and backcrossed them with CB-1. In the BC₄F₁ generation, we established a large population. We then employed a molecular marker (P371) developed based on the sequencing data to screen TMV-resistant plants for chromosomal recombination events within the target region (the linkage drag of the N gene from N. glutinosa). Such recombination events would lead to the partial elimination of the linkage drag. For clarity, we designated TMV-resistant plants containing the target recombination events as “recombinant plants,” while those without such events were termed “non-recombinant plants.” The identified recombinant plants were self-pollinated to produce the corresponding BC₄F₂ lines. Within these BC₄F₂ lines, we used a linked molecular marker (NR5) developed based on the sequencing data to identify plants with a homozygous N-gene genotype. These plants were further self-pollinated to obtain the corresponding N-gene homozygous lines, which are the recombinant lines with reduced linkage drag wanted. The above procedure is illustrated in Figure 1 .

To precisely determine the locations of recombination events occurred within the target region, the “homology analysis of N-contig” approach described earlier was applied to the identified recombinant lines. Specifically, each recombinant line was subjected to deep sequencing (approximately 40 Gb in depth), and the resulting reads were mapped to all contigs of 0970A. In addition, a non-recombinant line was analyzed in the same manner for comparative purposes.

A total of 6,892,630 HiFi reads were obtained from sequencing the genome of 0970A, with an average read length of 17.3 kb, a read N50 of 17,254 bp, and an average sequencing depth of 29 times the genome. These HiFi reads were then assembled and redundancies were eliminated, resulting in 554 contigs. Among these contigs, the N₅₀ length was 81,799,233 bp, and the N₇₅ length was 53,933,392 bp. Notably, approximately 80 contigs comprehensively covered the entire genome. The assessment using BUSCO revealed that the assembled genome sequence exhibited an impressive 99% completeness, underscoring the exceptional quality of the sequencing and assembly process.

Upon aligning the N gene sequence with all contigs of 0970A, it was determined that the N gene resided on contig PTG000094L, hereafter referred to as 094L (the detailed information of 094L, including full sequence and predicted genes, can be found at https://ngdc.cncb.ac.cn/genbase/search/gb/C_AA098042). The length of the N gene sequence used as the query was 7,345 bp, exhibiting a perfect match with the target sequence (E-value = 0, Score = 13,564). Notably, one end of contig 094L was identified as a telomere sequence with a repeat unit of TTTAGGG, suggesting that the contig represented a chromosomal terminus. This contig spanned 6.47 Mb in length, with the N gene located at 2.94 Mb from the non-telomere end. For convenience, we hereafter call the sequence of the contig 094L from the non-telomere end to the N gene and that from the telomere end to the N gene as the left side of N gene (LSN) and right side of N gene (RSN), respectively ( Figure 2A ).

A total of 167 genes were predicted on contig 094L, with 62 on the LSN and 104 on the RSN ( Supplementary Table S1 ). Interestingly, within the LSN, three genes were also annotated to encode the TMV resistance protein N, namely, ptg94lG0000310.t1 (abbreviated as G31) located at 1.94 Mb, ptg94lG0000390.t1 (G39) at 2.26 Mb, and ptg94lG0000440.t1 (G44) at 2.43 Mb ( Supplementary Table S1 ). However, they all were shorter than the full-length of the N gene sequence, measuring 6,419 bp (G31), 4,732 bp (G39), and 7,279 bp (G44), respectively, suggesting that they are only homologs of the N gene.

Since 094L was identified as the terminal sequence of a chromosome, we tried to locate the N gene on the chromosome by conducting collinearity analysis between 094L and the terminal sequences of all chromosomes of CB-1. The results revealed a significant collinearity between a segment of 094L and chromosome 11, referred to as Nt11 in accordance with the established chromosome numbering system of the K326 genome sequence ( Figure 2B ), suggesting that the N gene must be located on this chromosome. It is important to highlight that the collinearity between 094L and the corresponding CB-1 genome sequence is not consistent throughout. Specifically, the LSN exhibited a good collinearity with the CB-1 sequence, whereas the RSN displayed no collinearity with the CB-1 sequence ( Figure 2B ), implying that the RSN might not originate from N. tabacum.

To precisely determine the homology of the sequences flanking the N gene, we performed deep sequencing of the genomes of CB-1, Coker 176, and N. glutinosa. We then mapped their reads to all contigs of 0970A and analyzed the mapping results for 094L (the reads of CB-1, Coker 176, and N. glutinosa mapped to 094L can be found at https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA030530). The reads from Coker 176 uniformly covered the entire 094L with an average depth of approximately 20 ( Figure 3A ), indicating that the whole 094L is homologous to Coker 176. This finding aligns with the fact that Coker 176 is the donor parent of the N gene in 0970A. In contrast, the reads from CB-1 were predominantly mapped to the LSN and sparsely to the RSN ( Figure 3B ), suggesting that the LSN is homologous to CB-1. This result is consistent with the collinearity analysis ( Figure 2B ). Notably, the reads from CB-1 were evenly distributed across the LSN with an average depth of about 45, except for a small region (spanning from 1.55 Mb to 1.70 Mb) where the depth was significantly lower ( Figure 3B ). This indicates a structural variation between CB-1 and 0970A in that region, which the collinearity analysis suggests is an insertion in 0970A ( Figure 2B ). Conversely, the reads from N. glutinosa were primarily mapped to the RSN and barely to the LSN ( Figure 3C ), suggesting that the RSN is homologous to N. glutinosa. It is noteworthy that the depth distribution of N. glutinosa reads on the RSN was highly uneven. Although there was an overall depth plateau at around 30, reflecting the average sequencing depth of the N. glutinosa genome, there were many locations with extremely high depths (over 150), especially in the region beyond 4.92 Mb, where these high depths appeared almost continuously. The extremely high depths likely resulted from the presence of multiple copies of sequences dispersed throughout the N. glutinosa genome. Since these reads of multiple copies could not be mapped to other contigs of 0970A, they were exclusively mapped to the RSN, leading to the observed high depths.

The results presented above demonstrate that the linkage drag associated with the N gene from N. glutinosa in 0970A (as well as in Coker 176) is mainly restricted to the RSN of contig 094L. To facilitate removing the associated linkage drag while introducing the N gene from 0970A into CB-1, we developed a PCR-based molecular marker, designated as P371, located approximately 3.46 Mb from the N gene at the RSN end (close to the telomere) of contig 094L ( Figure 2A ). The P371 marker (forward primer: 5’-CTAAAAACGGGTGAAGGGGT-3’; reverse primer: 5’-ACATTCGGAGGTCAGGTCAAT-3’) is a dominant marker that produces a specific 162 bp amplification product in the donor parent 0970A but not in the recipient parent CB-1. Using the in vitro identification method, we identified 4,762 TMV-resistant plants (carrying the N gene) from a population of 10,013 BC₄F₁ plants. All resistant plants were subsequently genotyped using the P371 marker, and one plant was found to lack the specific amplification product ( Supplementary Figure S1 ). This finding indicates a recombination event between the P371 marker and the N gene, suggesting that a portion of the linkage drag, including the telomere, has been replaced with the corresponding segment from CB-1 in this plant. Therefore, this is a desired recombinant plant.

To identify individuals with a homozygous genotype for the N gene, we developed an SSR marker, designated as NR5, located approximately 740 kb from the N gene on the LSN of contig 094L ( Figure 2A ). The NR5 marker (forward primer: 5’-GGTCTATGTGCACCCTACCCT-3’; reverse primer: 5’-TTGGTTACGGCACGTAATGTT-3’) produces a 200 bp amplification product in 0970A and a shorter amplification product in CB-1 ( Supplementary Figure S2 ). Since NR5 is closely linked to the N gene, individuals exhibiting only the 0970A-specific band are likely to be homozygous for the N gene. Using the NR5 marker, we screened 30 progeny plants derived from self-pollination of the recombinant plant and identified 8 plants with a homozygous genotype for the N gene. These 8 plants were further confirmed to carry the N gene using the N gene-specific marker N1+N2 (forward primer N1: 5’-CGTCGACACATTATGCCATC-3’; reverse primer N2: 5’-GAGGGGTCTTACCCCATTGT-3’) (Lewis et al., 2005). Self-pollination of these plants yielded 8 recombinant lines with a homozygous genotype for the N gene. Using a similar approach, we also established 3 non-recombinant control lines (homozygous for the N gene) from a single non-recombinant plant.

To precisely characterize the recombination events in the target region of the recombinant lines, we conducted deep sequencing analysis on 3 recombinant lines and 1 non-recombinant line (the reads of these lines mapped to 094L can be found at https://ngdc.cncb.ac.cn/bioproject/browse/PRJCA030530). As expected, the non-recombinant line retained the complete 094L sequence ( Figure 4A ), whereas the 3 recombinant lines exhibited extensive replacement of the RSN sequence, with the crossover point located approximately at the position of 3 Mb, only ~60 kb away from the N gene ( Figures 4B-D ). This finding indicates that the residual linkage drag associated with the N gene in the recombinant lines is minimal.

Furthermore, based on SNP markers revealed by deep sequencing, we found that the 3 recombinant lines contained 97.85%-98.72% of the genetic components from the recipient parent CB-1 ( Supplementary Table S2 ), indicating a high degree of genetic background restoration. However, the non-recombinant line had a comparatively lower restoration rate of the CB-1 genetic background, containing only 93.42% of the genetic components from CB-1 ( Supplementary Table S2 ).