Preprocessing and quality control
Filtration of low-quality sequences with less than 25 scores led to the achievement of 105 Gb of cleaned data, which indicated 6 Gb per sample. The cleaned sequence GC content varied from 42.1 to 42.7%. The mapping rates for the cleaned sample reads were opposed to the reference genome sequence ranging from 91.6 to 97.8%. The sequencing quality and gene expression levels were generally consistent across the sequenced samples.
KEGG enrichment analysis of meta DEGs
Four studies with 70 samples were introduced to the process of meta-analysis. Non-biological heterogenicity was removed from studies by applying batch effect correction to gene expression data. PCA results showed that normalization and standardization declined batch effects and noises by which the direct merging of the datasets was facilitated (Fig. 1).
Rankprod method extracted 7897 DEGs (FDR ≤ 0.05) of which 3393 and 4504 were associated with up- and down-regulation, respectively. KEGG pathway analysis showed that among 7897 DEGs, 67 were directly and indirectly involved in lignin biosynthesis pathway. The mentioned pathway was categorized in the biosynthesis of secondary metabolites class (nta01110) (Table S1). The secondary metabolite pathway, which included phenylpropanoid group, was significantly enriched during the classification of KEGG. Generally, the meta-analysis results were consistent with lignin biosynthesis under drought stress (Bang et al., 2022; Choi et al., 2023; Chun et al., 2021; Sun et al., 2020; Zhao et al., 2021). Therefore, they were considered as a basis for lignin biosynthesis-related candidate genes detection.
In the current study, 67 differential genes were identified from 13 gene families located on the lignin biosynthesis pathway throughout the N. tabacum genome. Maximum numbers of isoforms were identified for HCT (Shikimate O-Hydroxy-Cinnamoyl transferase), CCoAoMT (Cinnamyl Alcohol Dehydrogenase), and PRX (Peroxidase) that contained the maximum homologous isoforms (21) amongst all monolignol genes in the current study. Five genes encoding the same catalytic enzymes showed an obvious difference in the regulation. However, all DEG isomers of 4CL, CCR, F5H, F6H, CCoAoMT were up-regulated, while all DEGs isomers of PAL and COMT and UGT72E were down-regulated (Table 3, Fig. 2).
hierarchical cluster analysis of the 67 meta-genes related to lignin biosynthesis pathway showed that these genes could be divided into two main clusters, based on their up-regulation or down-regulation in response to drought treatment. Further analysis showed that some of the lignin-related genes were co-regulated in response to drought treatment, while others exhibited more variable expression patterns (Fig. 3).
Insilico analysis of lignin associated Meta-gene relationships
Promoter analysis
Among promoter region of 61 lignin meta-genes, twenty conserved motifs were discovered by MEME algorithm (Table S2 and Table S3) with lengths ranging from 20 to 50 bp. Identified motifs are distributed on both positive and negative strands. Motif No. 3 was a frequent motif shared by the majority (83.6%) of the promoter regions of lignin meta-genes. However, motifs 6, 13 for down- and 14 for up-regulated genes were specified motifs. Motif 6 and 13 were located in positive strand, while motif 14 located in both strands. The result suggested that local DNA-sequence elements and their positional context may play a crucial role in transcriptional regulation of lignin genes.
Lignin biosynthesis-related co‑expressed modules
To achieve a better understanding of lignin biosynthesis-involved gene regulation network within N. tabacum, we conducted WGCNA for Meta-genes using the dynamic tree-cutting algorithm. Totally, 13591 meta-genes were identified from the datasets and 5789 genes were transferred to WGCNA after genes screening with low coefficient of variation (CV). As it could be observed in Fig. S1b, the main branches of dendrogram identified 14 different modules, which were represented by branches of different colors ranging from 32 to 1760 genes per module. The hierarchical topological overlap matrix (TOM)-derived clustering of the meta-genes was represented in Fig. S1a. Identification of genes with high correlation in various phases of lignin biosynthesis was carried out through WGCNA. Considering the functional enrichment of the identified modules, two lignin- enriched significant modules were selected (Table S4). The brown and yellow modules contained 19 and 15 lignin genes, respectively, which were reported in previous studies (Liu et al., 2021). Results of KEGG enrichment analysis carried out on 2 selected modules showed that there was a considerable number of the genes expressed in the metabolic pathway; moreover, a dramatic change of plants metabolism was observed in the secondary metabolic pathway that was resulted from stress adaptation. According to GO survey (Fig. S2a,b) carried out on a brown module, the enrichment in cinnamic acid biosynthesis, negative regulation of cytokinin-activated signaling pathway, and regulation of monopolar cell growth were observed. Also, the co-expressed gene in yellow module was enriched in sulfate reduction, intracellular auxin transport, and lignin biosynthetic process. Therefore, it was concluded that functionally-characterized lignin biosynthesis genes were distributed in 2 modules. The above-mentioned modules included 42 of all 66 genes that were preferentially expressed in the leaf; also, they were active in the lignin biosynthesis co-expressed with vital biosynthetic pathways, such as photosynthesis and amino acid metabolism, as well as a number of secondary metabolite pathways including isoquinoline, stilbenoid, gingerol, and flavonoid biosynthesis correlations (Fig. 4a,b). David's database detected several uncharacterized genes correlated with lignin pathway. Therefore, we visualized the co-expression networks of this module using Cytoscape. According to this visualization, each node represented a gene and connecting lines (edges) between genes indicated the coexpression (Fig. 5a,b). A direct interaction was observed between lignin biosynthesis genes and flavonoid genes, plant hormone, starch and pentose, as well as glucuronate pathway genes. the yellow module showed the co-expressions of cytokinin regulation, flavonoid biosynthesis, stilbenoid, diarylheptanoid, gingerol biosynthesis, Galactose metabolism, sulfur-containing amino acid, as well as cysteine- / tyrosine-rich proteins (Didi et al., 2015; Diehl and Brown, 2014; Hoshi and Heinemann, 2001; Khadr et al., 2020; Van Acker et al., 2013).
Identification of TFs
Transcription factors (TFs) play important roles in plant biological regulatory networks. Moreover, they are initially evaluated through comparative co-expression analyses using structural genes and target genes. There were also 72 additional regulatory TFs observed in the brown module that belonged to 21 families and consisted 69 down-regulated and 3 up-regulated TFs. In the mentioned module, most of the TFs were related to C2C2 and bHLH families (Fig. S3). Moreover, MYB306, MYB1R1, and bZIP61 showed the maximum connection to other genes (Fig. 6a). Their weights were also between 0.1 to 0.36. Considering the co-expression analysis for the yellow module, there were 93 up-regulated TF targeted genes derived from 15 families. The maximum TFs quantities was related to MYB (11) and AP2 (8) family; moreover, MYB308, ATHB12, ODORANAT, and MYB306 showed the maximum connection to other genes. Their weights were between 0.1 to 0.36 (Fig. 6b).
Co-expression network analysis revealed that the most of the monolignol biosynthesis-connected TF transcripts were derived from C2C2 and MYB classes. A direct interaction was observed between the genes of lignin biosynthesis and GATA8, bZIP61. MYB35, and bZIPs uncharacterized TFs. The results showed that the regulation of lignin biosynthesis was carried out through MYB (MYB83, MYB46), and downstream TF that included MYB58, MYB63, MYB85, MYB4, MYB32, and MYB7. By contrast, MYB4, MYB32, and MYB7 negatively regulated the expression of lignin biosynthetic genes.
Empirical Assessment of Candidate Meta-genes
Phenotype assessmeLignin contents were measured among leaves of “OR25”, “NC100”, “Basma”, “K326”, and “Burly” tobacco genotypes. The results showed that ‘Burly’ genotype contains the lowest amount of lignin accumulation while ‘OR25’, ‘Basma’, and ‘K326’ genotypes do not demonstrate a significant difference. However, lignin content in ‘NC100’ leaves and cells from cell suspension culture was approximately 1.5 times higher than ‘Burly’ (Fig. 7; Fig. 9). Under Drought conditions, cell growth was significantly affected by PEG 25% treatment. Under this condition, the growth rate of NC100 and Burly was decreased by 53% and 67%, respectively (Fig. 8). Cells of NC100 and Burly accumulated lignin content up to 3.34 and 2.35 times under PEG 25% in comparison with control condition, respectively (Fig.9). However, the lignin content of ‘NC100’ cells was relatively higher than that in 'Burly' cells in both treatment levels.
Expression profiles of genes involved in lignin biosynthesis
To confirm the relationship between 8 meta-genes correlated to lignin biosynthesis and lignin content, we performed qRT-PCR for tobacco cell suspensions affected by drought stress in two genotypes (Fig. 10,11). The results showed that expression level of CAD2 and ATH12 was induced in NC100 and expression of CAD2 and CCR was upregulated in Burly. On the other hand, ODO1, CAD6, CCoAoMT, COMT, F5H, and CCR were downregulated in NC100 while in Burly genotype ATH12, CoMT and F5H were down-regulated and the other genes did not change in response to drought stress treatment.