Based on our previous research (Tian et al., 2018; Pi et al., 2022), we selected the hybrids and their parents with different nicotine and potassium contents, and 12 hybrids were made from 7 parents according to the incomplete diallel crossing (NCII) method( Supplementary File 1, Table 1 ). All materials were provided by Guizhou Key Laboratory of Tobacco Quality Research. From 2021-2022, in the tobacco research base of Guizhou University, Yangwu Township, Anshun City, Guizhou Province, China, Seeds were sown in greenhouses using the floating seedling method, and the seedlings with six real leaves were sampled. The methods were as follows: Firstly, 3 seedlings with similar growth vigor were selected for root phenotype determination. After that, fresh biological samples were frozen with liquid nitrogen and stored in an ultra-low temperature refrigerator at - 80 ℃. These fresh samples were used for transcriptome sequencing analysis and RT-qPCR.

The cleaned complete root system was scanned with the root scanner LA2400 to obtain a clear image of the root structure. Use WinRhizo Pro 2021a software (Regent, Canada) to count tobacco root tips.

Total RNA was extracted from the root tissue using TRIzol^® Reagent (Plant RNA Purification Reagent for plant tissue) according to the manufacturer’s instructions (Invitrogen), and genomic DNA was removed using DNase I (TaKara). Then RNA quality was determined by 2100 Bioanalyser (Agilent) and quantified using the ND-2000 (NanoDrop Technologies). Only high-quality RNA sample (OD260/280 = 1.8~2.2, OD260/230≥2.0, RIN≥6.5, 28S:18S≥1.0, >1μg) was used to construct sequencing library. Paired-end RNA-seq sequencing library was sequenced with the NovaSeq 6000 sequencer (2 × 150bp read length).

By preprocessing the raw reads, short sequences with a length <25 nt and low quality sequences were removed. After preprocessing, the obtained reads were mapped to the Nicotiana tabacum sequenced cultivar K326 genome (Edwards et al., 2017) using the splice-aware mapping tool, Tophat2 (Kim et al., 2019). RSEM was used to quantify gene abundances (Li and Dewey, 2011). Essentially, differential expression analysis was performed using the DESeq2 (Love et al., 2014), DEGs with |log2FC| ≥ 1 and p ≤ 0.05 were considered to be significantly different expressed genes. Gene ontology (GO) functional-enrichment analysis were carried out by Goatools (Klopfenstein et al., 2018). In addition, we use PlantTFDB (http://planttfdb.gao-lab.org/) Transcription factor analysis of genes.

In order to divide the differentially expressed genes into 12 expression modes, according to the definition of Rapp et al. (2009). we used STEM software (Ernst and Bar-Joseph, 2006) to perform additional analysis on the differentially expressed genes. To allow clustering on a reasonable number of possible model profiles, the parameter for “STEM clustering method,” model profiles “was set to 50 and 2 was selected as the” maximum unit change between time points “, and” Minimum Absolute Expression Change “was set to 0 to classify the expression patterns of all differentially expressed genes.

The typical genes of Arabidopsis root development were used to identify core genes in tobacco. BLASTP identified 6 core genes with E-values truncated at 1E-10. To further classify root development genes in tobacco, phylogenetic analyses were performed using homologous protein sequences from Arabidopsis thaliana and Triticum aestivum L. The phylogenetic tree was constructed by MEGA 11.0(https://www.megasoftware.net/), The ultrafast bootstrap with 1000 replicates was conducted to obtain the supporting values for each bootstrap of the tree. Finally, gene trees were visualized and colored using iTOL(https://itol.embl.de/).

For the validation of the transcriptome data, we randomly selected 6 DEGs for RT-qPCR (PCR quantitative real-time) analysis. Supplementary File 1: Table 16 lists the genes and corresponding primers used in RT-qPCR. The total RNA used in the sequencing was reversely transcribed to obtain cDNA, which was used as a template to amplify the target genes, and the RT-qPCR experiment was conducted. The RNA was reversely transcribed into cDNA using the FastKing reverse transcription kit (Tiangen, China) as per the manual’s instructions. qPCR was performed using the BIO-RAD CFX96 Real-Time qPCR system. The relative expression of each gene was calculated using 2^−△△Ct (Livak and Schmittgen, 2001).

SPSS 25.0 software was used for statistical analyses. The variance analysis of root phenotype was carried out using the Duncan’s new multiple range method (P < 0.05). Based on the root phenotype, the values of mid-parent heterosis (MPH) were calculated according to the following method, MPH  ( % ) = ( F 1 − MP MP ) × 100 , where F₁represents the value of the first generation of hybrid and MP represents the average value of parents ( P 1 + P 2 2 ) .

Based on the screening results of tobacco root number heterosis ( Supplementary File 1, Table 1 ), the root numbers of the 12 F₁ hybrids all showed strong heterosis values, indicating that tobacco F₁ hybrids have general heterosis. Two highly dominant F₁ hybrids, GJ (G70 × Jiucaiping No.2) and KJ (K326 × Jiucaiping No.2), were selected from these F₁ hybrids for further analysis (their median parent values were not significantly different). At the tobacco seedling stage, we observed that the root development of the two hybrids were significantly better than that in their parents ( Figure 1A ). The number of root systems in GJ was 623, whereas G and J had 515 and 343, respectively; the number of roots in KJ was 878, whereas K and J had 711 and 343, respectively. The results showed that the average number of roots of both hybrids were significantly higher than that of their parents ( Figure 1B ), and the heterosis was 45.66% (GJ) and 67.12% (KJ) ( Figure 1C ).

Using two hybrids, GJ and KJ, with strong heterosis and their three parents as sequencing materials, 15 cDNA libraries were constructed and sequenced on the Illumina Novaseq 6000 platform. Each sample generated 76-98 million clean reads in RNA-Seq, which were then used for further analysis. Clean reads were mapped to the K326 genome, and about 95% of clean reads could be mapped to the genome ( Supplementary File 1, Table 2 ). A total of 39,345 genes were expressed in 15 analyzed samples ( Supplementary File 1, Table 3 ). Among them, 8.30% of the genes were very highly expressed (FPKM ≥ 50), 64.19% and 27.51% of the genes were moderately (1 ≤ FPKM ≤ 10) and highly expressed (10 ≤ FPKM ≤ 50), respectively. These results reflect that there is no significant difference in expression coverage between hybrids and parents during development. However, according to hierarchical clustering analysis( Figure 2A ) and principal component analysis of parents and hybrids ( Figure 2B ), both parents are closely clustered, and hybrids can be clearly divided into two different expression situations. These results show that there are great differences in gene expression between parents and hybrids during root development.

In the case of significance level P ≤ 0.05 and | log₂ fold change | ≥ 1, we identified 1430 genes up-regulated and 1697 genes down-regulated between G (female parent) and GJ (hybrid), 2015 genes up-regulated and 1809 genes down-regulated between J (male parent) and GJ (hybrid) ( Figures 3A, B ; Supplementary File 1, Table 4 ), 1613 genes up-regulated and 1922 genes down-regulated between the middle parent and GJ hybrid, and 302 genes were up-regulated and 291 genes were down-regulated between G and J ( Supplementary File 1, Table 5 ). We identified 3082 genes up-regulated and 3007 genes down-regulated between K (female parent) and KJ (hybrid), 4118 genes up-regulated and 3810 genes down-regulated between J (male parent) and KJ (hybrid) ( Figures 3C, D ), 3683 up-regulated genes and 3874 down-regulated genes between middle parent and KJ hybrid, and 283 genes up-regulated and 342 genes down-regulated between K and J ( Supplementary File 1, Table 6 ). This demonstrates that there are relatively few differential genes between the parents of the two strong heterosis hybrids, but there are a large number of differential genes between the parents and the hybrids, which may be the reason for root number heterosis.

To further analyze the DEGs of hybrids and parents, the genes were divided into 12 expression patterns (P1-P12, Figure 4A ) according to the definition of Rapp et al. (2009). The genes in the P1 and P2 modes are additive, P3-P6 modes are dominant, and P7-P12 modes show overdominant expression, in which the genes in P7-P9 are down-regulated overdominant and P10-P12 are up-regulated overdominant. Among the overdominant genes, 238, 376, and 1564 genes of GJ hybrid showed up-regulated overdominant expression patterns( Supplementary File 1, Table 7 ), and 516, 857, and 733 genes showed down-regulated overdominant expression patterns; there were 373, 1160, and 2750 genes in KJ hybrid that showed up-regulated overdominant expression patterns( Supplementary File 1, Table 8 ), and 1375, 1995, and 772 genes that showed down-regulated overdominant expression patterns ( Figure 4B ). Among these non-additively expressed genes (P3-P12), GJ hybrid had the highest proportion (84.22%) in the expression pattern of overdominant (P7-P12), and KJ hybrid had the highest proportion (90.25%) in the expression pattern of overdominant (P7-P12) ( Figure 4C ). Therefore, our results show that the overdominant expression advantage is the main reason for the formation of heterosis of tobacco roots.

Gene Ontology (GO) enrichment analysis was performed on the overdominant expression gene sets in GJ and KJ hybrids. We found that most overdominant genes were involved in “biological process” (55.3 – 56.1%). In GJ hybrids, most of them were significantly enriched in the jasmonic acid biosynthetic process, riboflavin metabolic process, riboflavin biosynthetic process, aminoglycan metabolic process, and jasmonic acid metabolic process ( Supplementary File 1, Table 9 ; P < 0.05). In KJ hybrids, the overdominant genes were significantly enriched in the glucosamine containing compound catabolic process, chitin metallic process, amino sugar catabolic process, chitin catabolic process, and aminoglycan catabolic process ( Supplementary File 1, Table 10 ; P < 0.05). According to the GO enrichment analysis of GJ and KJ co overdominant expression genes, most of the genes involved in cell differentiation, proliferation, cell growth and other processes were significantly enriched. It indicated that the overdominant expression of genes in the process of cell differentiation and proliferation was related to the formation of root heterosis.

According to the analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) metabolic pathways, 4284 overdominant genes in GJ were annotated to 117 pathways ( Supplementary File 1, Table 11 ; P < 0.05), whereas 8425 dominant genes in KJ were annotated to 123 pathways ( Supplementary File 12, Table 12 ; P < 0.05). As shown in Figures 5A, B , the overdominant gene enrichment pathways of the two hybrids are mainly “phylpropanoid biosynthesis” (9.72 and 5.27%, respectively), “plant hormone signal transduction” (4.97 and 4.93%, respectively), “MAPK signaling pathway - plant” (4.91 and 4.06%, respectively), and “start and cross metadata” (2.48 and 2.56%, respectively). In addition, 84 and 158 overdominant genes of the two hybrids were related to “plant hormone signal transduction”, and 42 and 82 genes were related to “starch and sucrose metabolism”. Based on KEGG analysis, we further studied specific functional genes related to root development.

Combined with functional analysis, in order to strengthen our understanding of hormone signaling during heterosis development in tobacco root growth, we focused on the overdominant genes of hormone signal transduction related to the biosynthesis of AUX, CTK, ABA, ET, and SA ( Figure 6 ; Supplementary File 1, Table 13 ). In the AUX signal transduction pathway, a total of 11 genes were overexpressed in F₁ hybrids (GJ and KJ), including 3 that were annotated to the AUX1 gene, 5 that were annotated to the GH3 gene, and 3 that were annotated to the AUX response gene SAUR. In the CTK signal transduction pathway, one receptor histidine kinase (CRE1) gene was negatively overexpressed in F₁ hybrids (GJ and KJ). One gene was annotated to A-type response regulator (A-ARR), which showed up-regulated overdominant expression in hybrids. In the ABA signal transduction pathway, 9 genes were found to be overexpressed in F₁ hybrids, including 5 that were annotated to protein phosphatase (PP2C) and 4 encoding Subcrosse non-fermenting 1-related protein kinase 2 (SnRK2). Interestingly, the genes related to PP2C showed up-regulated overdominant expression in both hybrids, of which 2 genes were highly expressed (10 ≤ FPKM ≤ 50). In the ET signal transduction pathway, 4 genes were annotated to ET receptor sensor (ETR2), one was annotated to EIN3, three were annotated to bind F-box protein 1/2 (EBF1/2), and four were annotated to ET responsive transcription factor 1/2 (ERF1/2). These ET signal transduction pathway results show that the genes related to ETR2, EIN3 and EBF1/2 are up-regulated and overdominantly expressed in both hybrids. For the SA signal transduction pathway, a total of 11 genes were overdominantly expressed in F₁ hybrids (GJ and KJ). Among them, 1 gene was annotated to TGA transcriptional regulator, which showed down-regulated overdominant expression in the roots of hybrids, and 7 genes were annotated to protein PR-1. It is worth noting that these 7 PR-1 related genes not only showed an overdominant expression pattern in the hybrid, but also had high differential expression multiples and highly expressed (FPKM ≥ 50).

Among the common overdominant expression genes of the two hybrids, various GO categories such as cell population promotion, cell differentiation, multiple cell growth and primary meristem tissue development are significantly enriched. In the process of cell differentiation, most genes show up-regulated overdominant expression ( Figure 7 ), in which Nitab4.5_ 0002714g0010、Nitab4.5_ 0008353g0010、Nitab4.5_ 0002918g0030 showed high expression in hybrids. In the process of cell population promotion and cell population promotion, 3 and 5 genes were up-regulated overdominant expression in the hybrid, indicating that cell population promotion and cell population promotion were important processes of root heterosis. In the process of primary meristem issue development, the two genes are both expressed down-regulated in the hybrid. In general, a lot of changes have taken place in cell population promotion, cell differentiation and multiple cell growth of hybrids, which may be the reason for the formation of the heterosis of root.

A weighted gene co-expression network analysis (WGCNA) was performed on all genes to thoroughly characterize the expression of genes related to root number. WGCNA identified 10 vivid modules (including 19338 genes) ( Figure 8A ; Supplementary File 1, Table 14 ). The module analysis revealed that there was a significant correlation between the brown module and the heterosis of the number of hybrid roots ( Figure 8B ). The genes in the “brown” module are mainly related to cysteine biological process, protein serine/threonine kinase activity, ethylene binding, etc. Further analysis of 6 genes with degree>700 in the brown module, Nitab4.5_ 0000119g0080 gene is Lysine tRNA ligase, Nitab4.5_ 0004384g0070、Nitab4.5_ 0011528g0020、Nitab4.5_ 0003282g0020 gene is Serine/threonine -/dual specificity protein kinase, Nitab4.5_ 0002818g0040 gene is NAD dependent emergency/dehydratase, Nitab4.5_ 0009232g0030 gene is Mitochondrial carrier protein ( Figure 8C ). Interestingly, these genes are up-regulated overdominant expression patterns. We believe that these genes may be related to the formation of root heterosis.

RT-qPCR was used to analyze 6 overdominant genes that have been shown to be involved in root development in other crops. In brief, protein sequences of well-characterized genes of each family from Arabidopsis thaliana were used as BLAST query sequences ( Supplementary File 1, Table 15 ), and the following sequences in tobacco were filtered by stringent threshold (E-value <1E-25, % identity>30 and % query coverage>30). In tobacco, protein sequences of root development exhibited a clear topological relationship, with high bootstrap values (Arabidopsis thaliana, Triticum aestivum L.)( Figure 9A ). Actin gene was used as a reference gene to standardize the expression levels of each gene. Supplementary File 1: Table 16 lists the genes and corresponding primers used in the RT-qPCR. The results showed that ( Figure 9B ), the genes of root development also showed a trend of overdominant expression. Although the expression multiples of RT-qPCR data and RNA-seq sequencing results are slightly different, they basically show a consistent expression trend, so RNA-seq data has certain reliability.