Integrated transcriptomic analysis reveals a transcriptional regulation network for the biosynthesis of lignin in Nicotiana tabacum in drought stress response

Integrated transcriptomic analysis reveals a transcriptional regulation network for the biosynthesis of lignin in Nicotiana tabacum in drought stress response | Research Square window.SnipcartSettings = { analytics: { enabled: false } }; (function() { var accessVector = localStorage.getItem('access_vector') || ''; window.dataLayer = window.dataLayer || []; if (accessVector) { window.dataLayer.push({ user: { profile: { profileInfo: { snid: accessVector } } } }); } })(); (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0],j=d.createElement(s),dl=l!='dataLayer'?'&l='+l:'';j.async=true;j.src='https://www.googletagmanager.com/gtm.js?id='+i+dl;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-K279D39R'); Browse Preprints In Review Journals COVID-19 Preprints AJE Video Bytes Research Tools Research Promotion AJE Professional Editing AJE Rubriq About Preprint Platform In Review Editorial Policies Our Team Advisory Board Help Center Sign In Submit a Preprint Cite Share Download PDF Research Article Integrated transcriptomic analysis reveals a transcriptional regulation network for the biosynthesis of lignin in Nicotiana tabacum in drought stress response Maryam Rashidifar, Hossein Askari, Ali Moghadam This is a preprint; it has not been peer reviewed by a journal. https://doi.org/ 10.21203/rs.3.rs-4101335/v1 This work is licensed under a CC BY 4.0 License Status: Posted Version 1 posted You are reading this latest preprint version Abstract Lignin plays a crucial impact on the production of phenolic compounds in tobacco smoke, which have potential health implications associated with tobacco use. The meta-analysis of RNA-seq studies along with high-resolution expression analysis on Nicotiana tabacum clarified a conserved distinctive expression pattern of lignin gene network. According to the results, 67 DEGs associated with lignin biosynthesis network were identified of which 17 genes were introduced by meta-analysis. WGCNA showed 14 clusters for the meta-genes. Various TF families and a number of regulatory factors were identified as the most likely candidate genes associated with the lignin metabolic pathway. Eight major meta-genes were evaluated by using qRT-PCR in two tobacco genotypes with different lignin content under drought stress conditions. Genotype NC100 (high lignin content) and Burly (low lignin content) in response to PEG upregulated CAD 2, ATH12 and CAD 2, CCR , respectively. Despite the accumulation of lignin, the expression levels of CCoAOMT , F5H , COMT , and ODO1 were reduced in both genotypes. The study's exploration into the complex nature of these pathways and meta-analysis highlights the importance of adopting a more comprehensive approach to gene discovery. It suggests that conducting additional individual investigations is crucial for enhancing the reliability and comprehensiveness of gene identification within intricate metabolite pathways. Lignin Meta-analysis RNA-seq Tobacco WGCNA Figures Figure 1 Figure 2 Figure 3 Figure 4 Figure 5 Figure 6 Figure 7 Figure 8 Figure 9 Figure 10 Introduction Meta-analysis allows us to reliably combine comparable high-throughput data to improve the statistical power of studies for detecting data bias and finally to pinpoint high-confidence causal genes regulating a metabolite concentration. Advancements in high-throughput technologies have facilitated the efficient and comprehensive characterization of gene expression profiles responsible for different conditions and attributes. Developing novel data integration approaches help us to integrate relevant data from multiple sources and drive the best decision based on the studies. Lignin is an abundant aromatic biopolymer with a complex mixture of phenolic compounds which is found in plant cell wall (Zeng et al., 2017 ). Thermal degradation of lignin generates some phenolic compounds with toxic and well-known threats to human health (Smith and Hansch, 2000 ). Reduction of toxic phenolic compounds in cigarette smoke is a major public interest in tobacco industry (Dagnon et al., 2011). Molecular breeding of tobacco is a powerful approach to partly control potential reduction of lignin compounds. Content, composition and structure of lignin can be altered by changes in gene regulation (Franke et al., 2000 ; Gu et al., 2020 ; Gui et al., 2011 ; Hu et al., 2009 ; Nakano et al., 2015 ; Sibout et al., 2017 ; Song et al., 2021 ; Van Acker et al., 2017 ; Vanholme et al., 2008 ; Wang et al., 2010 ; Yang et al., 2017 ; Zhao et al., 2010 ). A growing list of candidate regulatory genes with some putative roles in the regulation of monolignol biosynthesis (Dauwe et al., 2007 ; Koutaniemi et al., 2007 ; Sibout et al., 2005 ) demonstrates that lignin deposition is determined by complex and diverse interactions among multiple genes in different times and cell types. Lignin accumulation was consistent with relative expression level of lignin genes and this is dependent on the plant tissue types (Song et al., 2021 ). Lignin deposition in plant cell wall will be induced by drought stress (Bang et al., 2019 ; Barros et al., 2015 ; Cesarino, 2019 ; Choi et al., 2023 ; Lee et al., 2016 ). So that water deficit could significantly alter the expression level of lignin biosynthesis genes and increased lignin content in the both low- and high-lignin-accumulating genotypes (dos Santos et al., 2015 ). As shown, lignin biosynthesis pathway is regulated by key genes such as phenylalanine ammonia-lyase ( PAL ), cinnamate 4-hydroxylase ( C4H ), 4-coumarate: CoA ligase ( 4CL ), hydroxycinnamoyl-CoA shikimate/quinate hydroxycinnamoyl transferase ( HCT ), p-coumarate 3-hydroxylase ( C3H ), caffeoyl-CoA O-methyltransferase ( CCoAOMT ), cinnamoyl-CoA reductase ( CCR ), ferulate 5-hydroxylase ( F5H ), caffeic acid O-methyltransferase ( COMT ), and cinnamyl alcohol dehydrogenase ( CAD ) which are differentially involved in lignin composition and deposition (Boerjan et al., 2003 ; Bonawitz and Chapple, 2010 ). Lignin content, composition and structure can be altered using modification of some key genes (Franke et al., 2000 ; Gu et al., 2020 ; Gui et al., 2011 ; Hu et al., 2009 ; Nakano et al., 2015 ; Sibout et al., 2017 ; Song et al., 2021 ; Van Acker et al., 2017 ; Vanholme et al., 2008 ; Wang et al., 2010 ; Yang et al., 2017 ; Zhao et al., 2010 ). As shown, lignin content is regulated by a huge number of genes which are connected in a gene regulatory network revealing the most complex and extensive interactions between genes. To answer the question which gene(s) is a relevant candidate to be engineered, we assessed integration of transcriptome data to achieve a uniform gene expression network to find the most important determinant lignin gene(s). In the present study, two systems biology approaches of meta-analysis and co-expression gene network analysis were employed to combine multiple RNA-seq datasets from limited studies (four individual datasets) to find and experimentally verify the potential responsive genes involved in lignin biosynthesis induced by drought stress in low- and high-lignin-accumulating tobacco cells. Materials and methods Data collection and processing RNA-seq datasets of drought-exposed N. tabacum were retrieved from Gene Expression Omnibus database. The RNA-seq datasets were selected from studies with similar stress conditions and physiological phases, non-transgenic and non-mutant plants, along with distinct control conditions and biological replications. On the other hand, investigations on different cultivars of N. tabacum which include leaf response to drought stress were applied. Totally, 4 studies with 70 samples were selected for data analysis (Table 1). The qualification of source data was carried out using FastQC. These clean paired-end reads were mapped to the tobacco reference genome ( N. tabacum Ntab-TN90v, https://www.ncbi.nlm.nih. gov/genome/425) using HISAT (v2.0.5) (Kim et al., 2015). HTSeq (ver. 0.11.1) was also implemented to quantify the gene expression levels (Anders et al., 2015). RNA-Seq data was normalized through Read per Kilobase of Exon per Million Fragments Mapped Reads (RPKM) values. After production of the mentioned initial data, we unified databases and estimated the potential non-biological experimental variations (batch effects) derived from combining multiple datasets by Principal Component Analysis (PCA). Finally, the batch effects were removed using SVA v.3.26 ComBat method, which is an empirical Bayes method. Identification of differentially expressed genes The meta-analysis which was conducted using Random Effects Rank Product (Rankprod), as well as MetaDE R package 1.0.5 (Wang et al., 2012) was implemented to detect the robust up/down-regulated genes. By using the effect size combination method, both effect sizes and P-values were achieved (Marot et al., 2009). Genes with a P-value ≤ 0.05 were considered as DEGs. Gene enrichment analysis To illustrate the metabolic category of DEGs, we use Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2016) pathway enrichment analyses (Moriya et al., 2007); moreover, FDR ≤ 0.05 was considered for the statistically significant pathways. DEGs functional classification was conducted through related KEGG pathway enrichment analysis using KOBAS software (http:// kobas. cbi.pku.edu.cn). Co-expression network analysis DEGs with similar expression patterns were determined by applying WGCNA on the matrix of meta-analyzed DEGs normalized expression values (Langfelder and Horvath, 2008). Specifically, the top-ranked DEGs with a 75% coefficient of variation were reserved (Zhao et al., 2020). Optional threshold parameter (b = 14) was selected, and WGCNA R package was utilized to transfer the intergenic correlation coefficient into the adjacent coefficient. Then, the adjacency matrix was converted to a topological overlap (TO) matrix using TOM similarity algorithm. The hierarchical clustering of genes and consequent application of gene dendrogram were carried out with the purpose of module detection through the dynamic tree cut method (minModuleSize = 30). The Maximal Clique Centrality (MCC) algorithm was considered the most efficient method of hub node identification within a co-expression network. It is noteworthy that each node MCC was calculated using CytoHubba which is a plugin in cytoscape (Chin et al., 2014). Subsequently, an investigation was carried out on all hub genes using DAVID (https://david.ncifcrf.gov/) database in order to evaluate the general gene ontology. Identification of miRNAs, transcription factors and promoter motifs To discover potential miRNAs with the capability of lignin genes-associated DEGs targeting, psRNATarget server (http://plantgrn.noble.org/psRNATarget/) was applied. To detect regulators (TFs) in co-expressed modules, the genes were blasted against the plants TF and PK identifier, and iTAK (http://bioinfo.bti.cornell.edu/cgi bin/itak/index.cgi) (Zheng et al., 2016). In addition, 1 kb upstream flanking regions of co-expressed candidate genes were extracted from ensemble plants (https://plants.ensembl.org/index.html). MEME (meme.nbcr.net/meme/intro.html) (Bailey et al., 2006) was implemented to detect the conserved motifs located onto the sequences and predetermined parameters, which do not include the maximum number of motifs (20) with a threshold E-value < 1e-4. Plant materials The experiments were carried out on “OR25”, “NC100”, “Basma”, “K326” and “Burly” tobacco genotypes. The genotypes were obtained from the Tobacco Research Institute, Tirtash, Behshahr, Iran. The plant seeds were kept for germination in a 1/2 strength MS (Murashige and Skoog, 1962). Leaves from 4- to 6-week-old seedlings were used as explants. We selected two genotypes based on lignin content through thioglycolic acid (TGA) method (Graham and Graham, 1991). For raising callus cultures, MS basal medium was fortified with 1.5 mgL −1 each of NAA in combination with 0.5 mgL −1 BA and pH value 5.7, as this combination was found effectual in the preliminary experiments of this laboratory. For raising cell suspension cultures, friable calli were transferred into flasks containing 50 mL of liquid MS medium supplemented with 2 mg/L concentrations of 2,4-D. The flask with cultures were agitated on a horizontal shaker at 100 rpm at 25 ± 2 °C, under dark condition. Stress conditions and sampling To induce osmotic stress, polyethylene glycol (PEG-6000) was added to the basal medium at 25% (w/v), resulting in drought. The samples were harvested after 48 hours of drought induction. Additionally, when significant retarded growth of cells was observed in response to PEG treatments, we took the samples. In addition, the thioglycolic acid (TGA) method was used to estimate the total lignin content in the plant biomass (Graham and Graham, 1991) Selection of candidate genes In this study, we used an integrative data analysis result to select 8 candidate genes for validation in tobacco cells. The genes selected were CCR , CAD , CCoAOMT , COMT , and F5H , which are known to play important roles in lignin biosynthesis. We also focused on lignin-specified transcription factors, ODO1 and HB12 , which are co-regulated in yellow module and are known to be involved in lignin biosynthesis. Although the role of these genes was studied before test to identified the expression pattern of them in lignin biosynthesis in different drought severities and genotypes. Validation of candidate genes Total RNA was extracted by RNA extraction kit (Ana Cell, Iran) according to the manufacturer’s instructions. The RNA samples were treated with deoxyribonuclease I (DNase I; SinaColon, Iran) to remove residual genomic DNA. The treated RNA was rechecked by NanoDrop and agarose gel. Then, 1 μg of DNase-treated RNA was used for first-strand cDNA synthesis using Ana Cell cDNA synthesis kit (Anacell tec, Iran). Gene expression analysis was performed using Real-time PCR (Bio-Rad, Hercules, CA) using specific primers of genes (Table 2). All primers were designed by Allele ID software. The reactions were provided at Parstoos SYBR green master mix Kit (Parstoos, Iran), and real-time PCR system (MIC, Australia). The transcript abundance was calculated from three biological and three technical replicates with 18S rRNA internal control. Three biological and technical replicates were performed for each sample to verify the results of the qRT-PCR, and relative gene expression levels were quantified via the 2-ΔΔCt method (Livak and Schmittgen, 2001) . Statistical analysis Three technical replicates were applied for each sample in qRT-PCR. The data were analyzed using the SAS 9.4 software (SAS Institute Inc., USA) at a P-value of ≤ 0.05. A mean comparison analysis was performed using the Two-way ANOVA and LSD test. The graphs were plotted using GraphPad prism software v8 (GraphPad Software Inc., San Diego, CA). Results 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 assessme Lignin 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. Discussion Lignin is a multifunctional polyphenolic compound which accumulates in response to various physiological and environmental cues. Lignin gene network is an extremely complex and naturally vibrant network. On the other hand, Lignin presence is not essential for plant cell life but its crucial role will be commenced along with developing a secondary cell wall. In spite of essential and important function of lignin in plant physiologic and ecologic life, lignin can play a negative role in some applications of plant materials such as paper industry (Baghel et al., 2020 ) bioethanol production(Broda et al., 2022 ), human nutrition (Tao et al., 2020 ) and pyrolysis during smoking (Liao et al., 2017 ). Depending on the plant production goals, breeding strategies will be designed to increase or decrease lignin content of target plant tissues or organs. Here, we concentrate on the tobacco cells to find highly effective gene(s) regulating and functioning in the lignin biosynthesis pathway by using systems biology approaches and to make a high-resolution assessment of the candidate genes to design a genetic engineering strategy and change lignin accumulation with the least manipulation of the target genome. Our data showed that 67 meta-genes which belong to 13 gene families from phenylpropanoid pathway were significantly regulated in response to drought stress. Comparative analysis of metagenes and individual studies showed that a number of metagenes were higher than all differentially expressed genes of individual studies so that in study one, there was no lignin-related gene and in the other studies with the basis of genotype and stress treatment completely different genes from lignin biosynthetic pathway could be find (Table 6). Of 13 gene families, eight candidate genes which were distinctly major players of lignin biosynthesis pathway were selected and assessed. Genotype NC100 and Burly in response to 25% PEG up-regulated CAD 2, ATH12 and CAD 2, CCR , respectively. Despite significant regulation of some other genes, amount of expression changes was not more than two-fold changes. NC100 as a higher lignin accumulator down-regulated all other functional ( CAD6 , CCoAOMT , COMT , F5H , and CCR ) and regulatory ( ODO1 ) genes in response to drought stress. On the other hand, CCR which plays an essential role in lignin biosynthesis in plants (Cui et al., 2022 ; Goujon et al., 2003 ) showed the highest upregulation in Burly while NC100 significantly decreased expression of CCR . As we know CAD increased the synthesis of cinnamyl alcohols and is thought to be a specific lignification marker (Mitchell et al., 1994 ). Previous studies reported that lignin content had a positive correlation with CAD expression (Georgii et al., 2017 ; Hu et al., 2009 ; Jin et al., 2014 ; Liu et al., 2020 ). As reported in Table 3, meta-analysis demonstrated that all of the lignin-related gene isoforms differentially up- and down-regulated in response to drought stress with the exception of CCR , F5H , and CCoAOMT . This study presents a major conclusion by challenging the prevalent assumption of meta-analysis as an acceptable method for identifying important genes inside complex pathways including the lignin pathway. By exploring into the intricate nature of these pathways and revealing the scarcity of individual studies available for meta-analysis, the study highlights the limits and emphasizes the importance of taking a more complete approach to gene discovery. This innovative concept underlines the importance of undertaking more individual investigations to improve the reliability and comprehensiveness of gene identification within complicated metabolite pathways. Declarations Acknowledgments The authors would like to thank the Plant Biotechnology Department in Shahid Beheshti University and Institute of Biotechnology for supporting this research and the Bioinformatics Research Group in the College of Agriculture (Shiraz University). Funding The authors received no financial support for the research, authorship, and publication of this article. Competing Interests The authors declare no conflict of interest Authors' contributions Maryam Rashidifar performing the work, contributed to data acquisition and manuscript drafting. Hossein Askari design and supervised the research. Ali Moghadam contributed to design and data acquisition. All authors read and approved the final manuscript. References Anders S., Pyl P.T., Huber W. (2015) HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics 31:166-9. DOI: 10.1093/bioinformatics/btu638. Baghel S., Sahariah B.P., Anandkumar J. (2020) Bioremediation of Lignin-Rich Pulp and Paper Industry Effluent, in: M. Shah and A. Banerjee (Eds.), Combined Application of Physico-Chemical & Microbiological Processes for Industrial Effluent Treatment Plant, Springer Singapore, Singapore. pp. 261-278. Bailey T.L., Williams N., Misleh C., Li W.W. (2006) MEME: discovering and analyzing DNA and protein sequence motifs. Nucleic acids research 34:W369-W373. Bang S.W., Choi S., Jin X., Jung S.E., Choi J.W., Seo J.S., Kim J.-K. (2022) Transcriptional activation of rice CINNAMOYL-CoA REDUCTASE 10 by OsNAC5, contributes to drought tolerance by modulating lignin accumulation in roots. Plant Biotechnology Journal 20:736-747. DOI: https://doi.org/10.1111/pbi.13752. Bang S.W., Lee D.-K., Jung H., Chung P.J., Kim Y.S., Choi Y.D., Suh J.-W., Kim J.-K. (2019) Overexpression of OsTF1L, a rice HD-Zip transcription factor, promotes lignin biosynthesis and stomatal closure that improves drought tolerance. Plant Biotechnology Journal 17:118-131. DOI: https://doi.org/10.1111/pbi.12951. Barros J., Serk H., Granlund I., Pesquet E. (2015) The cell biology of lignification in higher plants. Ann Bot 115:1053-74. DOI: 10.1093/aob/mcv046. Boerjan W., Ralph J., Baucher M. (2003) Lignin biosynthesis. Annual review of plant biology 54:519-546. Bonawitz N.D., Chapple C. (2010) The genetics of lignin biosynthesis: connecting genotype to phenotype. Annu Rev Genet 44:337-63. DOI: 10.1146/annurev-genet-102209-163508. Broda M., Yelle D.J., Serwańska K. (2022) Bioethanol Production from Lignocellulosic Biomass-Challenges and Solutions. Molecules 27. DOI: 10.3390/molecules27248717. Cesarino I. (2019) Structural features and regulation of lignin deposited upon biotic and abiotic stresses. Curr Opin Biotechnol 56:209-214. DOI: 10.1016/j.copbio.2018.12.012. Chin C.-H., Chen S.-H., Wu H.-H., Ho C.-W., Ko M.-T., Lin C.-Y. (2014) cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC systems biology 8:1-7. Choi S.J., Lee Z., Kim S., Jeong E., Shim J.S. (2023) Modulation of lignin biosynthesis for drought tolerance in plants. Frontiers in Plant Science 14. DOI: 10.3389/fpls.2023.1116426. Chun H.J., Lim L.H., Cheong M.S., Baek D., Park M.S., Cho H.M., Lee S.H., Jin B.J., No D.H., Cha Y.J., Lee Y.B., Hong J.C., Yun D.J., Kim M.C. (2021) Arabidopsis CCoAOMT1 Plays a Role in Drought Stress Response via ROS- and ABA-Dependent Manners. Plants (Basel) 10. DOI: 10.3390/plants10050831. Cui W., Zhuang Z., Jiang P., Pan J., Zhao G., Xu S., Shen W. (2022) Characterization, Expression Profiling, and Biochemical Analyses of the Cinnamoyl-CoA Reductase Gene Family for Lignin Synthesis in Alfalfa Plants. International Journal of Molecular Sciences 23:7762. Dauwe R., Morreel K., Goeminne G., Gielen B., Rohde A., Van Beeumen J., Ralph J., Boudet A.M., Kopka J., Rochange S.F., Halpin C., Messens E., Boerjan W. (2007) Molecular phenotyping of lignin-modified tobacco reveals associated changes in cell-wall metabolism, primary metabolism, stress metabolism and photorespiration. Plant J 52:263-85. DOI: 10.1111/j.1365-313X.2007.03233.x. Didi V., Jackson P., Hejátko J. (2015) Hormonal regulation of secondary cell wall formation. Journal of Experimental Botany 66:5015-5027. DOI: 10.1093/jxb/erv222. Diehl B.G., Brown N.R. (2014) Lignin cross-links with cysteine- and tyrosine-containing peptides under biomimetic conditions. J Agric Food Chem 62:10312-9. DOI: 10.1021/jf503897n. dos Santos A.B., Bottcher A., Kiyota E., Mayer J.L., Vicentini R., Brito Mdos S., Creste S., Landell M.G., Mazzafera P. (2015) Water stress alters lignin content and related gene expression in two sugarcane genotypes. J Agric Food Chem 63:4708-20. DOI: 10.1021/jf5061858. Franke R., McMichael C.M., Meyer K., Shirley A.M., Cusumano J.C., Chapple C. (2000) Modified lignin in tobacco and poplar plants over‐expressing the Arabidopsis gene encoding ferulate 5‐hydroxylase. The Plant Journal 22:223-234. Georgii E., Jin M., Zhao J., Kanawati B., Schmitt-Kopplin P., Albert A., Winkler J.B., Schäffner A.R. (2017) Relationships between drought, heat and air humidity responses revealed by transcriptome-metabolome co-analysis. BMC Plant Biology 17:120. DOI: 10.1186/s12870-017-1062-y. Goujon T., Ferret V., Mila I., Pollet B., Ruel K., Burlat V., Joseleau J.P., Barrière Y., Lapierre C., Jouanin L. (2003) Down-regulation of the AtCCR1 gene in Arabidopsis thaliana: effects on phenotype, lignins and cell wall degradability. Planta 217:218-28. DOI: 10.1007/s00425-003-0987-6. Graham M.Y., Graham T.L. (1991) Rapid accumulation of anionic peroxidases and phenolic polymers in soybean cotyledon tissues following treatment with Phytophthora megasperma f. sp. glycinea wall glucan. Plant Physiology 97:1445-1455. Gu H., Wang Y., Xie H., Qiu C., Zhang S., Xiao J., Li H., Chen L., Li X., Ding Z. (2020) Drought stress triggers proteomic changes involving lignin, flavonoids and fatty acids in tea plants. Scientific Reports 10:15504. DOI: 10.1038/s41598-020-72596-1. Gui J., Shen J., Li L. (2011) Functional Characterization of Evolutionarily Divergent 4-Coumarate:Coenzyme A Ligases in Rice Plant Physiology 157:574-586. DOI: 10.1104/pp.111.178301. Hoshi T., Heinemann S. (2001) Regulation of cell function by methionine oxidation and reduction. The Journal of physiology 531:1-11. DOI: 10.1111/j.1469-7793.2001.0001j.x. Hu Y., Li W.C., Xu Y., Li G., Liao Y., Fu F.-L. (2009) Differential expression of candidate genes for lignin biosynthesis under drought stress in maize leaves. Journal of applied genetics 50:213-223. Jin Y., Zhang C., Liu W., Qi H., Chen H., Cao S. (2014) The cinnamyl alcohol dehydrogenase gene family in melon (Cucumis melo L.): bioinformatic analysis and expression patterns. PLoS One 9:e101730. DOI: 10.1371/journal.pone.0101730. Kanehisa M., Sato Y., Kawashima M., Furumichi M., Tanabe M. (2016) KEGG as a reference resource for gene and protein annotation. Nucleic acids research 44:D457-D462. Khadr A., Wang Y.H., Zhang R.R., Wang X.R., Xu Z.S., Xiong A.S. (2020) Cytokinin (6-benzylaminopurine) elevates lignification and the expression of genes involved in lignin biosynthesis of carrot. Protoplasma 257:1507-1517. DOI: 10.1007/s00709-020-01527-8. Kim D., Langmead B., Salzberg S.L. (2015) HISAT: a fast spliced aligner with low memory requirements. Nature Methods 12:357-360. DOI: 10.1038/nmeth.3317. Koutaniemi S., Warinowski T., Kärkönen A., Alatalo E., Fossdal C.G., Saranpää P., Laakso T., Fagerstedt K.V., Simola L.K., Paulin L., Rudd S., Teeri T.H. (2007) Expression profiling of the lignin biosynthetic pathway in Norway spruce using EST sequencing and real-time RT-PCR. Plant Mol Biol 65:311-28. DOI: 10.1007/s11103-007-9220-5. Langfelder P., Horvath S. (2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. DOI: 10.1186/1471-2105-9-559. Lee D.K., Jung H., Jang G., Jeong J.S., Kim Y.S., Ha S.H., Do Choi Y., Kim J.K. (2016) Overexpression of the OsERF71 Transcription Factor Alters Rice Root Structure and Drought Resistance. Plant Physiol 172:575-88. DOI: 10.1104/pp.16.00379. Liao J., Lu Z., Hu S., Li Q., Che L., Chen X.D. (2017) Effects of prewash on the pyrolysis kinetics of cut tobacco. Drying Technology 35:1368-1378. DOI: 10.1080/07373937.2017.1320803. Liu C., Yu H., Rao X., Li L., Dixon R.A. (2021) Abscisic acid regulates secondary cell-wall formation and lignin deposition in Arabidopsis thaliana through phosphorylation of NST1. Proceedings of the National Academy of Sciences 118:e2010911118. DOI: 10.1073/pnas.2010911118. Liu W., Jiang Y., Wang C., Zhao L., Jin Y., Xing Q., Li M., Lv T., Qi H. (2020) Lignin synthesized by CmCAD2 and CmCAD3 in oriental melon (Cucumis melo L.) seedlings contributes to drought tolerance. Plant Mol Biol 103:689-704. DOI: 10.1007/s11103-020-01018-7. Livak K.J., Schmittgen T.D. (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2− ΔΔCT method. methods 25:402-408. Mitchell H.J., Hall J.L., Barber M.S. (1994) Elicitor-Induced Cinnamyl Alcohol Dehydrogenase Activity in Lignifying Wheat (Triticum aestivum L.) Leaves. Plant Physiol 104:551-556. DOI: 10.1104/pp.104.2.551. Moriya Y., Itoh M., Okuda S., Yoshizawa A.C., Kanehisa M. (2007) KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res 35:W182-5. DOI: 10.1093/nar/gkm321. Murashige T., Skoog F. (1962) A Revised Medium for Rapid Growth and Bio Assays with Tobacco Tissue Cultures. Physiologia Plantarum 15:473-497. DOI: https://doi.org/10.1111/j.1399-3054.1962.tb08052.x. Nakano Y., Yamaguchi M., Endo H., Rejab N.A., Ohtani M. (2015) NAC-MYB-based transcriptional regulation of secondary cell wall biosynthesis in land plants. Frontiers in Plant Science 6. DOI: 10.3389/fpls.2015.00288. Sibout R., Eudes A., Mouille G., Pollet B., Lapierre C., Jouanin L., Séguin A. (2005) CINNAMYL ALCOHOL DEHYDROGENASE-C and -D are the primary genes involved in lignin biosynthesis in the floral stem of Arabidopsis. Plant Cell 17:2059-76. DOI: 10.1105/tpc.105.030767. Sibout R., Proost S., Hansen B.O., Vaid N., Giorgi F.M., Ho-Yue-Kuang S., Legée F., Cézart L., Bouchabké-Coussa O., Soulhat C., Provart N., Pasha A., Le Bris P., Roujol D., Hofte H., Jamet E., Lapierre C., Persson S., Mutwil M. (2017) Expression atlas and comparative coexpression network analyses reveal important genes involved in the formation of lignified cell wall in Brachypodium distachyon. New Phytol 215:1009-1025. DOI: 10.1111/nph.14635. Smith C.J., Hansch C. (2000) The relative toxicity of compounds in mainstream cigarette smoke condensate. Food Chem Toxicol 38:637-46. DOI: 10.1016/s0278-6915(00)00051-x. Song Z., Wang D., Gao Y., Li C., Jiang H., Zhu X., Zhang H. (2021) Changes of lignin biosynthesis in tobacco leaves during maturation. Functional Plant Biology 48:624-633. DOI: https://doi.org/10.1071/FP20244. Sun S.-C., Xiong X.-P., Zhang X.-L., Feng H.-J., Zhu Q.-H., Sun J., Li Y.-J. (2020) Characterization of the Gh4CL gene family reveals a role of Gh4CL7 in drought tolerance. BMC Plant Biology 20:125. DOI: 10.1186/s12870-020-2329-2. Tao J., Li S., Ye F., Zhou Y., Lei L., Zhao G. (2020) Lignin - An underutilized, renewable and valuable material for food industry. Crit Rev Food Sci Nutr 60:2011-2033. DOI: 10.1080/10408398.2019.1625025. Van Acker R., Déjardin A., Desmet S., Hoengenaert L., Vanholme R., Morreel K., Laurans F., Kim H., Santoro N., Foster C., Goeminne G., Légée F., Lapierre C., Pilate G., Ralph J., Boerjan W. (2017) Different Routes for Conifer- and Sinapaldehyde and Higher Saccharification upon Deficiency in the Dehydrogenase CAD1 Plant Physiology 175:1018-1039. DOI: 10.1104/pp.17.00834. Van Acker R., Vanholme R., Storme V., Mortimer J.C., Dupree P., Boerjan W. (2013) Lignin biosynthesis perturbations affect secondary cell wall composition and saccharification yield in Arabidopsis thaliana. Biotechnology for Biofuels 6:46. DOI: 10.1186/1754-6834-6-46. Vanholme R., Morreel K., Ralph J., Boerjan W. (2008) Lignin engineering. Current Opinion in Plant Biology 11:278-285. DOI: https://doi.org/10.1016/j.pbi.2008.03.005. Wang H., Avci U., Nakashima J., Hahn M.G., Chen F., Dixon R.A. (2010) Mutation of WRKY transcription factors initiates pith secondary wall formation and increases stem biomass in dicotyledonous plants. Proceedings of the National Academy of Sciences 107:22338-22343. Wang X., Kang D.D., Shen K., Song C., Lu S., Chang L.-C., Liao S.G., Huo Z., Tang S., Ding Y., Kaminski N., Sibille E., Lin Y., Li J., Tseng G.C. (2012) An R package suite for microarray meta-analysis in quality control, differentially expressed gene analysis and pathway enrichment detection. Bioinformatics (Oxford, England) 28:2534-2536. DOI: 10.1093/bioinformatics/bts485. Yang Q., He Y., Kabahuma M., Chaya T., Kelly A., Borrego E., Bian Y., El Kasmi F., Yang L., Teixeira P., Kolkman J., Nelson R., Kolomiets M., L Dangl J., Wisser R., Caplan J., Li X., Lauter N., Balint-Kurti P. (2017) A gene encoding maize caffeoyl-CoA O-methyltransferase confers quantitative resistance to multiple pathogens. Nature Genetics 49:1364-1372. DOI: 10.1038/ng.3919. Zeng Y., Himmel M.E., Ding S.-Y. (2017) Visualizing chemical functionality in plant cell walls. Biotechnology for Biofuels 10:263. DOI: 10.1186/s13068-017-0953-3. Zhao D., Luan Y., Shi W., Zhang X., Meng J., Tao J. (2021) A Paeonia ostii caffeoyl-CoA O-methyltransferase confers drought stress tolerance by promoting lignin synthesis and ROS scavenging. Plant Sci 303:110765. DOI: 10.1016/j.plantsci.2020.110765. Zhao J., Zou W., Hu T. (2020) Novel genes associated with folic acid-mediated metabolism in mouse: A bioinformatics study. PloS one 15:e0238940. Zhao Q., Wang H., Yin Y., Xu Y., Chen F., Dixon R.A. (2010) Syringyl lignin biosynthesis is directly regulated by a secondary cell wall master switch. Proceedings of the National Academy of Sciences 107:14496-14501. Zheng Y., Jiao C., Sun H., Rosli H.G., Pombo M.A., Zhang P., Banf M., Dai X., Martin G.B., Giovannoni J.J. (2016) iTAK: a program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Molecular plant 9:1667-1670. Tables Tables 1 to 6 are available in the Supplementary Files section. Supporting information Supporting information is not available with this version. Fig S1 Weighted gene co-expression network analysis (WGCNA) Fig S2 GO classification in 2 selected modules Fig S3 Transcription factor family classifications co-expressed throughout 2 selected modules Table S1 The enriched KEGG pathway of DEGs Table S2 sequences of conserved promoter motifs of up-regulated meta-genes Table S3 sequences of conserved promoter motifs of down-regulated meta-genes Table S4 List of modules in Weighted gene co-expression network analysis (WGCNA) Additional Declarations No competing interests reported. Supplementary Files Table1.xlsx Table2.xlsx table3.xlsx Table4.xlsx Table5.xlsx Table6.xlsx Cite Share Download PDF Status: Posted Version 1 posted You are reading this latest preprint version Research Square lets you share your work early, gain feedback from the community, and start making changes to your manuscript prior to peer review in a journal. As a division of Research Square Company, we’re committed to making research communication faster, fairer, and more useful. We do this by developing innovative software and high quality services for the global research community. Our growing team is made up of researchers and industry professionals working together to solve the most critical problems facing scientific publishing. Also discoverable on Platform About Our Team In Review Editorial Policies Advisory Board Help Center Resources Author Services Accessibility API Access RSS feed Manage Cookie Preferences © Research Square 2026 | ISSN 2693-5015 (online) Privacy Policy Terms of Service Do Not Sell My Personal Information {"props":{"pageProps":{"initialData":{"identity":"rs-4101335","acceptedTermsAndConditions":true,"allowDirectSubmit":true,"archivedVersions":[],"articleType":"Research Article","associatedPublications":[],"authors":[{"id":279774478,"identity":"398e96a3-0f81-4077-b5c0-daa7ac4bd096","order_by":0,"name":"Maryam Rashidifar","email":"","orcid":"","institution":"Shahid Beheshti University","correspondingAuthor":false,"prefix":"","firstName":"Maryam","middleName":"","lastName":"Rashidifar","suffix":""},{"id":279774480,"identity":"3ca8458c-017b-4248-9f8d-97a44118924f","order_by":1,"name":"Hossein Askari","email":"data:image/png;base64,iVBORw0KGgoAAAANSUhEUgAAAZAAAAAyAQMAAABI0h/eAAAABlBMVEX///8AAABVwtN+AAAACXBIWXMAAA7EAAAOxAGVKw4bAAAAzElEQVRIiWNgGAWjYBACAyBmbGBg4OFn4IGIsBGtRbKBVC0MBgd4iHSYOQPz448zau7JGJ8/e0yCocaOgU/6AH4tlg1sZpIbjhXzmN3IS5NgOJbMwMaXQMBhBxjMGB+wJQC18JhJMLAdYGAj5ECDA+yfPz74l8Bj3H8GqOUfUVp4DCQ3tiXwGDDkmEkwthGnpUxyZl8Cj8SNvGSLxL5kHmIctvljz7cEe/7+swdvfPhmJyffQ0ALg/wDJE4CAwOxsTMKRsEoGAWjAB8AAMA7OE7zNFxfAAAAAElFTkSuQmCC","orcid":"","institution":"Shahid Beheshti University","correspondingAuthor":true,"prefix":"","firstName":"Hossein","middleName":"","lastName":"Askari","suffix":""},{"id":279774483,"identity":"db763926-6893-45f1-b6d9-d8e3f8f3b089","order_by":2,"name":"Ali Moghadam","email":"","orcid":"","institution":"Shiraz University","correspondingAuthor":false,"prefix":"","firstName":"Ali","middleName":"","lastName":"Moghadam","suffix":""}],"badges":[],"createdAt":"2024-03-14 13:44:04","currentVersionCode":1,"declarations":"","doi":"10.21203/rs.3.rs-4101335/v1","doiUrl":"https://doi.org/10.21203/rs.3.rs-4101335/v1","draftVersion":[],"editorialEvents":[],"editorialNote":"","failedWorkflow":false,"files":[{"id":53016922,"identity":"17a777ad-e7dd-4351-8a0d-ac00ba19e20b","added_by":"auto","created_at":"2024-03-19 16:07:08","extension":"png","order_by":1,"title":"Figure 1","display":"","copyAsset":false,"role":"figure","size":1172474,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eIdentification and removal of the batch effects\u003c/strong\u003e. Scatter plots showed PCA analysis of normalized gene expression data \u003cstrong\u003eA.\u003c/strong\u003ebefore and\u003cstrong\u003e B\u003c/strong\u003e. after batch effect removal by combat.\u003c/p\u003e","description":"","filename":"1.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/6e59bd7aa0b387e52bb9dc67.png"},{"id":53018144,"identity":"a3eca467-afd0-4b81-8a12-6b728494d182","added_by":"auto","created_at":"2024-03-19 16:15:08","extension":"png","order_by":2,"title":"Figure 2","display":"","copyAsset":false,"role":"figure","size":37015,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eKEGG sub-pathway phenylpropanoid metabolism(meta-genes)\u003c/strong\u003e. The color scale represents genes involved in upregulation(red), down regulation(green) and genes involved in both up and down regulation (red and green).\u003c/p\u003e","description":"","filename":"2.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/0292d2a2ab55cf29438d2521.png"},{"id":53018149,"identity":"c56ce096-36a2-41b2-a961-7845dbf13d52","added_by":"auto","created_at":"2024-03-19 16:15:08","extension":"png","order_by":3,"title":"Figure 3","display":"","copyAsset":false,"role":"figure","size":1916238,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe hierarchical cluster analysis of 67 meta- genes related to lignin biosynthesis pathway\u003c/strong\u003e, which exhibited significantly different expression patterns by comparing the transcriptomes of control and drought treatment samples.The x axis shows samples and the y axis shows meta-gene related to lignin biosynthesis. Samples and gens are grouped using hierarchical clustering approach.\u003c/p\u003e","description":"","filename":"3.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/9f2683aa6736f401f7cbdb03.png"},{"id":53018146,"identity":"6d1025f3-80d2-42c3-aeb3-c485f56d6aab","added_by":"auto","created_at":"2024-03-19 16:15:08","extension":"png","order_by":4,"title":"Figure 4","display":"","copyAsset":false,"role":"figure","size":89535,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eKEGG pathway enrichment analysis of 2 top modules\u003c/strong\u003e. Panels a and b show all functional pathways associated with the mentioned genes in brown and yellow modules through KEGG analysis, respectively (\u003cem\u003eP-value\u003c/em\u003e≤ 0.05). Significantly enriched pathways of the mentioned modules are shown in Y-axis. Rich factor in the X-axis represents the enrichment levels. The larger value of Rich factor represents the higher level of enrichment. The color and size of the dot respectively indicate different \u003cem\u003eP-value\u003c/em\u003e and the number of target genes enriched in the corresponding pathway.\u003c/p\u003e","description":"","filename":"4.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/0f17de2dba4f319fbc0a3852.png"},{"id":53018151,"identity":"d0f0057a-18e5-48c7-8481-8d078797a0d7","added_by":"auto","created_at":"2024-03-19 16:15:08","extension":"png","order_by":5,"title":"Figure 5","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eThe cytoscape representation of co-expressed lignin genes\u003c/strong\u003e in \u003cstrong\u003eA\u003c/strong\u003e. brown and \u003cstrong\u003eB\u003c/strong\u003e. yellow modules. The mapping strategy of using the low number of interactions (neighborhood connectivity) corresponding to smaller size was used for the font of labeling\u003c/p\u003e","description":"","filename":"PlaceholderimageCopy2.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/615293a208b32db3b4e37ec2.png"},{"id":53016931,"identity":"aa355e6a-b0d3-4b52-afe7-a167a17d5c4c","added_by":"auto","created_at":"2024-03-19 16:07:08","extension":"png","order_by":6,"title":"Figure 6","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eFunctional genes and TFs interactions in co-expressed module\u003c/strong\u003e \u003cstrong\u003eA.\u003c/strong\u003e brown and \u003cstrong\u003eB\u003c/strong\u003e. yellow. Node size is in accordance with the connectivity degree (\u003cem\u003ei.e.\u003c/em\u003e, the number of edges connecting each node). The edges with higher weight are shown in red. TFs in the first shell is shown in green, and TFs in second shell are represented in light green.\u003c/p\u003e","description":"","filename":"PlaceholderimageCopy3.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/fdfd0c5111ae301b0c257cf0.png"},{"id":53019484,"identity":"674df839-cb7d-467e-b032-8b96ea971b67","added_by":"auto","created_at":"2024-03-19 16:23:08","extension":"png","order_by":7,"title":"Figure 7","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eAssessment of different Tobacco genotypes for leaf lignin content\u003c/strong\u003e. accumulation and OR25, Basma, and K326 do not demonstrate a significant difference of lignin content. Data are presented as mean ± SD of three replicates. (* \u003cem\u003eP-value\u003c/em\u003e\u0026lt; 0.05).\u003c/p\u003e","description":"","filename":"PlaceholderimageCopy4.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/9cac32105af7567203d39d52.png"},{"id":53019485,"identity":"aa8026a9-bf59-48c8-bb88-3fcc0c0406bc","added_by":"auto","created_at":"2024-03-19 16:23:08","extension":"png","order_by":8,"title":"Figure 8","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eCell suspension growth retarded drought stress levels.\u003c/strong\u003eError bars represent standard deviation among three independent replicates. Data are presented as mean ± SD of three replicates. Asterisks represent statistically significant difference between control condition and stress conditions (* \u003cem\u003eP-value\u003c/em\u003e \u0026lt; 0.05).\u003c/p\u003e","description":"","filename":"PlaceholderimageCopy5.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/9fdab7cb6f4cb0e1bc0f93de.png"},{"id":53016927,"identity":"691bcd59-3af5-4af3-b823-e3902ee86dc2","added_by":"auto","created_at":"2024-03-19 16:07:08","extension":"png","order_by":9,"title":"Figure 9","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eDetermination of lignin content in ‘NC100’ and ‘Burly’ cell suspension culture at drought stress levels.\u003c/strong\u003e Error bars represent standard deviation among three independent replicates. Data are presented as mean ± SD of three replicates. Asterisks represent statistically significant difference between control condition and stress conditions (* \u003cem\u003eP-value\u003c/em\u003e\u0026lt; 0.05).\u003c/p\u003e","description":"","filename":"PlaceholderimageCopy.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/efa6e2f27f04c324b13f493d.png"},{"id":53016933,"identity":"5d3923ca-d6ed-43f7-935b-2a50164426b6","added_by":"auto","created_at":"2024-03-19 16:07:09","extension":"png","order_by":10,"title":"Figure 10","display":"","copyAsset":false,"role":"figure","size":5713,"visible":true,"origin":"","legend":"\u003cp\u003e\u003cstrong\u003eQuantitative RT-PCR analysis of putative lignin biosynthesis genes in tobacco under control and drought stress conditions\u003c/strong\u003e. Relative expression was normalized based on average of control expression=1. Error bars represent SE (n=3). Asterisks represent statistically significant difference between control and stress condition (Two-way ANOVA; * p\u0026lt;0.05)\u003c/p\u003e","description":"","filename":"Placeholderimage.png","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/cf339de60f5bf0595b6c820a.png"},{"id":54021651,"identity":"7f129ad5-e25e-4e86-bc40-435c97bb145e","added_by":"auto","created_at":"2024-04-03 13:37:55","extension":"pdf","order_by":0,"title":"","display":"","copyAsset":false,"role":"manuscript-pdf","size":1502792,"visible":true,"origin":"","legend":"","description":"","filename":"manuscript.pdf","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/3e67c468-d658-4dad-997b-2e0e9ce62028.pdf"},{"id":53019482,"identity":"f0cec29b-204f-482c-a292-f4d851a4a43b","added_by":"auto","created_at":"2024-03-19 16:23:08","extension":"xlsx","order_by":1,"title":"","display":"","copyAsset":false,"role":"supplement","size":10420,"visible":true,"origin":"","legend":"","description":"","filename":"Table1.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/6d4cfe8401c30573325ffdeb.xlsx"},{"id":53016920,"identity":"2ce98585-da0a-48cf-bff8-e03e96dd75fb","added_by":"auto","created_at":"2024-03-19 16:07:08","extension":"xlsx","order_by":2,"title":"","display":"","copyAsset":false,"role":"supplement","size":9824,"visible":true,"origin":"","legend":"","description":"","filename":"Table2.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/c619d19c1fb345e6c784b982.xlsx"},{"id":53016918,"identity":"cfff206e-2b60-4260-a61c-d6529a4fcb74","added_by":"auto","created_at":"2024-03-19 16:07:08","extension":"xlsx","order_by":3,"title":"","display":"","copyAsset":false,"role":"supplement","size":10071,"visible":true,"origin":"","legend":"","description":"","filename":"table3.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/7dc41061984ce20eb03a992b.xlsx"},{"id":53019483,"identity":"062dac9f-e13d-4f6c-8a98-b14505d9432d","added_by":"auto","created_at":"2024-03-19 16:23:08","extension":"xlsx","order_by":4,"title":"","display":"","copyAsset":false,"role":"supplement","size":9887,"visible":true,"origin":"","legend":"","description":"","filename":"Table4.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/302521c2eeebc0c64f017029.xlsx"},{"id":53016928,"identity":"3cce2e50-cb61-404c-964b-f13660869ffe","added_by":"auto","created_at":"2024-03-19 16:07:08","extension":"xlsx","order_by":5,"title":"","display":"","copyAsset":false,"role":"supplement","size":9222,"visible":true,"origin":"","legend":"","description":"","filename":"Table5.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/76dd6b3d4be2bed385940eca.xlsx"},{"id":53016934,"identity":"b33cb4e6-b2b6-426f-b0e1-d119d9ac45c7","added_by":"auto","created_at":"2024-03-19 16:07:09","extension":"xlsx","order_by":6,"title":"","display":"","copyAsset":false,"role":"supplement","size":9814,"visible":true,"origin":"","legend":"","description":"","filename":"Table6.xlsx","url":"https://assets-eu.researchsquare.com/files/rs-4101335/v1/0b05c1eecbef3473a1ac9754.xlsx"}],"financialInterests":"No competing interests reported.","formattedTitle":"Integrated transcriptomic analysis reveals a transcriptional regulation network for the biosynthesis of lignin in Nicotiana tabacum in drought stress response","fulltext":[{"header":"Introduction","content":"\u003cp\u003eMeta-analysis allows us to reliably combine comparable high-throughput data to improve the statistical power of studies for detecting data bias and finally to pinpoint high-confidence causal genes regulating a metabolite concentration. Advancements in high-throughput technologies have facilitated the efficient and comprehensive characterization of gene expression profiles responsible for different conditions and attributes. Developing novel data integration approaches help us to integrate relevant data from multiple sources and drive the best decision based on the studies.\u003c/p\u003e \u003cp\u003eLignin is an abundant aromatic biopolymer with a complex mixture of phenolic compounds which is found in plant cell wall (Zeng et al., \u003cspan citationid=\"CR54\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Thermal degradation of lignin generates some phenolic compounds with toxic and well-known threats to human health (Smith and Hansch, \u003cspan citationid=\"CR44\" class=\"CitationRef\"\u003e2000\u003c/span\u003e). Reduction of toxic phenolic compounds in cigarette smoke is a major public interest in tobacco industry (Dagnon et al., 2011).\u003c/p\u003e \u003cp\u003eMolecular breeding of tobacco is a powerful approach to partly control potential reduction of lignin compounds. Content, composition and structure of lignin can be altered by changes in gene regulation (Franke et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Gu et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Gui et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Hu et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Nakano et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Sibout et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Song et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Van Acker et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Vanholme et al., \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Yang et al., \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Zhao et al., \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). A growing list of candidate regulatory genes with some putative roles in the regulation of monolignol biosynthesis (Dauwe et al., \u003cspan citationid=\"CR15\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Koutaniemi et al., \u003cspan citationid=\"CR31\" class=\"CitationRef\"\u003e2007\u003c/span\u003e; Sibout et al., \u003cspan citationid=\"CR42\" class=\"CitationRef\"\u003e2005\u003c/span\u003e) demonstrates that lignin deposition is determined by complex and diverse interactions among multiple genes in different times and cell types. Lignin accumulation was consistent with relative expression level of lignin genes and this is dependent on the plant tissue types (Song et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2021\u003c/span\u003e). Lignin deposition in plant cell wall will be induced by drought stress (Bang et al., \u003cspan citationid=\"CR5\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Barros et al., \u003cspan citationid=\"CR6\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Cesarino, \u003cspan citationid=\"CR10\" class=\"CitationRef\"\u003e2019\u003c/span\u003e; Choi et al., \u003cspan citationid=\"CR12\" class=\"CitationRef\"\u003e2023\u003c/span\u003e; Lee et al., \u003cspan citationid=\"CR33\" class=\"CitationRef\"\u003e2016\u003c/span\u003e). So that water deficit could significantly alter the expression level of lignin biosynthesis genes and increased lignin content in the both low- and high-lignin-accumulating genotypes (dos Santos et al., \u003cspan citationid=\"CR18\" class=\"CitationRef\"\u003e2015\u003c/span\u003e). As shown, lignin biosynthesis pathway is regulated by key genes such as phenylalanine ammonia-lyase (\u003cem\u003ePAL\u003c/em\u003e), cinnamate 4-hydroxylase (\u003cem\u003eC4H\u003c/em\u003e), 4-coumarate: CoA ligase (\u003cem\u003e4CL\u003c/em\u003e), hydroxycinnamoyl-CoA shikimate/quinate hydroxycinnamoyl transferase (\u003cem\u003eHCT\u003c/em\u003e), p-coumarate 3-hydroxylase (\u003cem\u003eC3H\u003c/em\u003e), caffeoyl-CoA O-methyltransferase (\u003cem\u003eCCoAOMT\u003c/em\u003e), cinnamoyl-CoA reductase (\u003cem\u003eCCR\u003c/em\u003e), ferulate 5-hydroxylase (\u003cem\u003eF5H\u003c/em\u003e), caffeic acid O-methyltransferase (\u003cem\u003eCOMT\u003c/em\u003e), and cinnamyl alcohol dehydrogenase (\u003cem\u003eCAD\u003c/em\u003e) which are differentially involved in lignin composition and deposition (Boerjan et al., \u003cspan citationid=\"CR7\" class=\"CitationRef\"\u003e2003\u003c/span\u003e; Bonawitz and Chapple, \u003cspan citationid=\"CR8\" class=\"CitationRef\"\u003e2010\u003c/span\u003e). Lignin content, composition and structure can be altered using modification of some key genes (Franke et al., \u003cspan citationid=\"CR19\" class=\"CitationRef\"\u003e2000\u003c/span\u003e; Gu et al., \u003cspan citationid=\"CR23\" class=\"CitationRef\"\u003e2020\u003c/span\u003e; Gui et al., \u003cspan citationid=\"CR24\" class=\"CitationRef\"\u003e2011\u003c/span\u003e; Hu et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Nakano et al., \u003cspan citationid=\"CR41\" class=\"CitationRef\"\u003e2015\u003c/span\u003e; Sibout et al., \u003cspan citationid=\"CR43\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Song et al., \u003cspan citationid=\"CR45\" class=\"CitationRef\"\u003e2021\u003c/span\u003e; Van Acker et al., \u003cspan citationid=\"CR48\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Vanholme et al., \u003cspan citationid=\"CR50\" class=\"CitationRef\"\u003e2008\u003c/span\u003e; Wang et al., \u003cspan citationid=\"CR51\" class=\"CitationRef\"\u003e2010\u003c/span\u003e; Yang et al., \u003cspan citationid=\"CR53\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Zhao et al., \u003cspan citationid=\"CR57\" class=\"CitationRef\"\u003e2010\u003c/span\u003e).\u003c/p\u003e \u003cp\u003eAs shown, lignin content is regulated by a huge number of genes which are connected in a gene regulatory network revealing the most complex and extensive interactions between genes. To answer the question which gene(s) is a relevant candidate to be engineered, we assessed integration of transcriptome data to achieve a uniform gene expression network to find the most important determinant lignin gene(s). In the present study, two systems biology approaches of meta-analysis and co-expression gene network analysis were employed to combine multiple RNA-seq datasets from limited studies (four individual datasets) to find and experimentally verify the potential responsive genes involved in lignin biosynthesis induced by drought stress in low- and high-lignin-accumulating tobacco cells.\u003c/p\u003e"},{"header":"Materials and methods","content":"\u003cp\u003e\u003cstrong\u003eData collection and\u0026nbsp;\u003c/strong\u003e\u003cstrong\u003eprocessing\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eRNA-seq datasets of drought-exposed \u003cem\u003eN. tabacum\u003c/em\u003e were retrieved from Gene Expression Omnibus database. The RNA-seq datasets were selected from studies with similar stress conditions and physiological phases, non-transgenic and non-mutant plants, along with distinct control conditions and biological replications. On the other hand, investigations on different cultivars of \u003cem\u003eN. tabacum\u003c/em\u003e which include leaf response to drought stress were applied. Totally, 4 studies with 70 samples were selected for data analysis (Table 1). The qualification of source data was carried out using FastQC. These clean paired-end reads were mapped to the tobacco reference genome (\u003cem\u003eN. tabacum\u003c/em\u003e Ntab-TN90v, \u0026nbsp;https://www.ncbi.nlm.nih. gov/genome/425) using HISAT (v2.0.5) (Kim et al., 2015). HTSeq (ver. 0.11.1) was also implemented to quantify the gene expression levels (Anders et al., 2015). RNA-Seq data was normalized through Read per Kilobase of Exon per Million Fragments Mapped Reads (RPKM) values. After production of the mentioned initial data, we unified databases and estimated the potential non-biological experimental variations (batch effects) derived from combining multiple datasets by Principal Component Analysis (PCA). Finally, the batch effects were removed using SVA v.3.26 ComBat method, which is an empirical Bayes method.\u003cem\u003e\u0026nbsp;\u003c/em\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIdentification of differentially expressed genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe meta-analysis which was conducted using Random Effects Rank Product (Rankprod), as well as MetaDE R package 1.0.5 (Wang et al., 2012) was implemented to detect the robust up/down-regulated genes. By using the effect size combination method, both effect sizes and \u003cem\u003eP-values\u003c/em\u003e were achieved (Marot et al., 2009). Genes with a \u003cem\u003eP-value\u003c/em\u003e \u0026le; 0.05 were considered as DEGs.\u003cstrong\u003e\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eGene enrichment analysis\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo illustrate the metabolic category of DEGs, we use Kyoto Encyclopedia of Genes and Genomes (KEGG) (Kanehisa et al., 2016) pathway enrichment analyses (Moriya et al., 2007); moreover, FDR \u0026le; 0.05 was considered for the statistically significant pathways. DEGs functional classification was conducted through related KEGG pathway enrichment analysis using KOBAS software (http:// kobas. cbi.pku.edu.cn).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCo-expression network analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eDEGs with similar expression patterns were determined by applying WGCNA on the matrix of meta-analyzed DEGs normalized expression values (Langfelder and Horvath, 2008). Specifically, the top-ranked DEGs with a 75% coefficient of variation were reserved (Zhao et al., 2020). Optional threshold parameter (b = 14) was selected, and WGCNA R package was utilized to transfer the intergenic correlation coefficient into the adjacent coefficient. Then, the adjacency matrix was converted to a topological overlap (TO) matrix using TOM similarity algorithm. The hierarchical clustering of genes and consequent application of gene dendrogram were carried out with the purpose of module detection through the dynamic tree cut method (minModuleSize = 30). The Maximal Clique Centrality (MCC) algorithm was considered the most efficient method of hub node identification within a co-expression network. It is noteworthy that each node MCC was calculated using CytoHubba which is a plugin in cytoscape (Chin et al., 2014). Subsequently, an investigation was carried out on all hub genes using DAVID (https://david.ncifcrf.gov/) database in order to evaluate the general gene ontology.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIdentification of miRNAs, transcription factors and promoter motifs\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo discover potential miRNAs with the capability of lignin genes-associated DEGs targeting, psRNATarget server (http://plantgrn.noble.org/psRNATarget/) was applied. To detect regulators (TFs) in co-expressed modules, the genes were blasted against the plants TF and PK identifier, and iTAK (http://bioinfo.bti.cornell.edu/cgi bin/itak/index.cgi) (Zheng et al., 2016).\u003cspan dir=\"RTL\"\u003e\u0026nbsp;\u003c/span\u003eIn addition, 1 kb upstream flanking regions of co-expressed candidate genes were extracted from ensemble plants (https://plants.ensembl.org/index.html). MEME (meme.nbcr.net/meme/intro.html) (Bailey et al., 2006) was implemented to detect the conserved motifs located onto the sequences and predetermined parameters, which do not include the maximum number of motifs (20) with a threshold \u003cem\u003eE-value\u003c/em\u003e \u0026lt; 1e-4.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePlant materials\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe experiments were carried out on \u0026ldquo;OR25\u0026rdquo;, \u0026ldquo;NC100\u0026rdquo;, \u0026ldquo;Basma\u0026rdquo;, \u0026ldquo;K326\u0026rdquo; and \u0026ldquo;Burly\u0026rdquo; tobacco genotypes. The genotypes were obtained from the Tobacco Research Institute, Tirtash, Behshahr, Iran. The plant seeds were kept for germination in a 1/2 strength MS (Murashige and Skoog, 1962). Leaves from 4- to 6-week-old seedlings were used as explants. We selected two genotypes based on lignin content through thioglycolic acid (TGA) method (Graham and Graham, 1991). For raising callus cultures, MS basal medium was fortified with 1.5 mgL\u003csup\u003e\u0026minus;1\u003c/sup\u003e each of NAA in combination with 0.5 mgL\u003csup\u003e\u0026minus;1\u003c/sup\u003e BA and pH value 5.7, as this combination was found effectual in the preliminary experiments of this laboratory. For raising cell suspension cultures, friable calli were transferred into flasks containing 50 mL of liquid MS medium supplemented with 2 mg/L \u0026nbsp; concentrations of 2,4-D. The flask with cultures were agitated on a horizontal shaker at 100 rpm at 25 \u0026plusmn; 2 \u0026deg;C, under dark condition.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eStress conditions and sampling\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo induce osmotic stress, polyethylene glycol (PEG-6000) was added to the basal medium at 25% (w/v), resulting in drought. The samples were harvested after 48 hours of drought induction. Additionally, when significant retarded growth of cells was observed in response to PEG treatments, we took the samples. In addition, the thioglycolic acid (TGA) method was used to estimate the total lignin content in the plant biomass (Graham and Graham, 1991)\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eSelection of candidate genes\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eIn this study, we used an integrative data analysis result to select 8 candidate genes for validation in tobacco cells. The genes selected were \u003cem\u003eCCR\u003c/em\u003e, \u003cem\u003eCAD\u003c/em\u003e, \u003cem\u003eCCoAOMT\u003c/em\u003e, \u003cem\u003eCOMT\u003c/em\u003e, and \u003cem\u003eF5H\u003c/em\u003e, which are known to play important roles in lignin biosynthesis. We also focused on lignin-specified transcription factors, \u003cem\u003eODO1\u003c/em\u003e and \u003cem\u003eHB12\u003c/em\u003e, which are co-regulated in yellow module and are known to be involved in lignin biosynthesis. Although the role of these genes was studied before test to identified the expression pattern of them in lignin biosynthesis in different drought severities and genotypes.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eValidation of candidate genes\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTotal RNA was extracted by RNA extraction kit (Ana Cell, Iran) according to the manufacturer\u0026rsquo;s instructions. The RNA samples were treated with deoxyribonuclease I (DNase I; SinaColon, Iran) to remove residual genomic DNA. The treated RNA was rechecked by NanoDrop and agarose gel. Then, 1 \u0026mu;g of DNase-treated RNA was used for first-strand cDNA synthesis using Ana Cell cDNA synthesis kit (Anacell tec, Iran). Gene expression analysis was performed using Real-time PCR (Bio-Rad, Hercules, CA) using specific primers of genes (Table 2). All primers were designed by Allele ID software. The reactions were provided at Parstoos SYBR green master mix Kit (Parstoos, Iran), and real-time PCR system (MIC, Australia). The transcript abundance was calculated from three biological and three technical replicates with \u003cem\u003e18S rRNA\u0026nbsp;\u003c/em\u003einternal control. Three biological and technical replicates were performed for each sample to verify the results of the qRT-PCR, and relative gene expression levels were quantified via the 2-\u0026Delta;\u0026Delta;Ct method\u0026nbsp;(Livak and Schmittgen, 2001)\u003cem\u003e.\u003c/em\u003e\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eStatistical analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThree technical replicates were applied for each sample in qRT-PCR. The data were analyzed using the SAS 9.4 software (SAS Institute Inc., USA) at a \u003cem\u003eP-value\u0026nbsp;\u003c/em\u003eof \u0026le; 0.05. A mean comparison analysis was performed using the \u0026nbsp;Two-way ANOVA and LSD test. The graphs were plotted using GraphPad prism software v8 (GraphPad Software Inc., San Diego, CA).\u003c/p\u003e"},{"header":"Results","content":"\u003cp\u003e\u003cstrong\u003ePreprocessing and quality control\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFiltration of low-quality sequences with less than 25 scores led to the achievement of 105 Gb of cleaned data, which indicated \u0026nbsp;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.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eKEGG enrichment analysis of meta DEGs\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eFour studies with 70 samples were introduced to the process of meta-analysis.\u003cspan dir=\"RTL\"\u003e\u0026nbsp;\u003c/span\u003eNon-biological heterogenicity was removed from studies by applying batch effect correction to gene expression data. PCA results showed that normalization and standardization\u0026nbsp;declined\u0026nbsp;batch effects and\u0026nbsp;noises by which\u0026nbsp;the direct merging of the datasets was facilitated (Fig. 1).\u003c/p\u003e\n\u003cp\u003eRankprod method extracted 7897 DEGs (\u003cem\u003eFDR\u003c/em\u003e \u0026le; 0.05) of which 3393 and 4504 were associated with up- and down-regulation, respectively.\u003cspan dir=\"RTL\"\u003e\u0026nbsp;\u003c/span\u003eKEGG pathway analysis showed that among 7897 DEGs, 67 were\u003cspan dir=\"RTL\"\u003e\u0026nbsp;\u003c/span\u003edirectly 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.\u003cu\u003e\u0026nbsp;\u003c/u\u003e\u003c/p\u003e\n\u003cp\u003eIn the current study, 67 differential genes were identified from 13 gene families located on the lignin biosynthesis pathway throughout the \u003cem\u003eN. tabacum\u003c/em\u003e genome. Maximum numbers of isoforms were identified for \u003cem\u003eHCT\u003c/em\u003e (Shikimate O-Hydroxy-Cinnamoyl transferase), \u003cem\u003eCCoAoMT\u003c/em\u003e (Cinnamyl Alcohol Dehydrogenase), and \u003cem\u003ePRX\u003c/em\u003e (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 \u003cem\u003e4CL\u003c/em\u003e, \u003cem\u003eCCR\u003c/em\u003e, \u003cem\u003eF5H\u003c/em\u003e, \u003cem\u003eF6H\u003c/em\u003e, \u003cem\u003eCCoAoMT\u003c/em\u003e were up-regulated, while all DEGs isomers of \u003cem\u003ePAL\u003c/em\u003e and \u003cem\u003eCOMT\u003c/em\u003e and \u003cem\u003eUGT72E\u003c/em\u003e were down-regulated (Table 3, Fig. 2).\u003c/p\u003e\n\u003cp\u003ehierarchical 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).\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eInsilico analysis of lignin associated Meta-gene relationships\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePromoter analysis\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eAmong 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.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eLignin biosynthesis-related co‑expressed modules\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo achieve a better understanding of lignin biosynthesis-involved gene regulation network within \u003cem\u003eN. tabacum\u003c/em\u003e, 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).\u0026nbsp;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\u0026nbsp;lignin\u0026nbsp;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\u0026apos;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).\u003cem\u003e\u0026nbsp;\u003c/em\u003eA 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). \u003cstrong\u003e\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eIdentification of TFs\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTranscription 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.\u0026nbsp;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).\u003c/p\u003e\n\u003cp\u003eCo-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,\u003cspan dir=\"RTL\"\u003e\u0026nbsp;\u003c/span\u003ebZIP61. MYB35, and bZIPs uncharacterized TFs.\u0026nbsp;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.\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eEmpirical Assessment of Candidate Meta-genes\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003ePhenotype assessme\u003c/strong\u003eLignin contents were measured among leaves of \u0026ldquo;OR25\u0026rdquo;, \u0026ldquo;NC100\u0026rdquo;, \u0026ldquo;Basma\u0026rdquo;, \u0026ldquo;K326\u0026rdquo;, and \u0026ldquo;Burly\u0026rdquo; tobacco genotypes. The results showed that \u0026lsquo;Burly\u0026rsquo; genotype contains the lowest amount of lignin accumulation while \u0026lsquo;OR25\u0026rsquo;, \u0026lsquo;Basma\u0026rsquo;, and \u0026lsquo;K326\u0026rsquo; genotypes do not demonstrate a significant difference. However, lignin content in \u0026lsquo;NC100\u0026rsquo; leaves and cells from cell suspension culture was approximately 1.5 times higher than \u0026lsquo;Burly\u0026rsquo; (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 \u0026lsquo;NC100\u0026rsquo; cells was relatively higher than that in \u0026apos;Burly\u0026apos; cells in both treatment levels. \u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eExpression profiles of genes involved in lignin biosynthesis\u0026nbsp;\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eTo 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.\u0026nbsp;\u003c/p\u003e"},{"header":"Discussion","content":"\u003cp\u003eLignin is a multifunctional polyphenolic compound which accumulates in response to various physiological and environmental cues. Lignin gene network is an extremely complex and naturally vibrant network. On the other hand, Lignin presence is not essential for plant cell life but its crucial role will be commenced along with developing a secondary cell wall. In spite of essential and important function of lignin in plant physiologic and ecologic life, lignin can play a negative role in some applications of plant materials such as paper industry (Baghel et al., \u003cspan citationid=\"CR2\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) bioethanol production(Broda et al., \u003cspan citationid=\"CR9\" class=\"CitationRef\"\u003e2022\u003c/span\u003e), human nutrition (Tao et al., \u003cspan citationid=\"CR47\" class=\"CitationRef\"\u003e2020\u003c/span\u003e) and pyrolysis during smoking (Liao et al., \u003cspan citationid=\"CR34\" class=\"CitationRef\"\u003e2017\u003c/span\u003e). Depending on the plant production goals, breeding strategies will be designed to increase or decrease lignin content of target plant tissues or organs. Here, we concentrate on the tobacco cells to find highly effective gene(s) regulating and functioning in the lignin biosynthesis pathway by using systems biology approaches and to make a high-resolution assessment of the candidate genes to design a genetic engineering strategy and change lignin accumulation with the least manipulation of the target genome. Our data showed that 67 meta-genes which belong to 13 gene families from phenylpropanoid pathway were significantly regulated in response to drought stress. Comparative analysis of metagenes and individual studies showed that a number of metagenes were higher than all differentially expressed genes of individual studies so that in study one, there was no lignin-related gene and in the other studies with the basis of genotype and stress treatment completely different genes from lignin biosynthetic pathway could be find (Table\u0026nbsp;6). Of 13 gene families, eight candidate genes which were distinctly major players of lignin biosynthesis pathway were selected and assessed.\u003c/p\u003e \u003cp\u003eGenotype NC100 and Burly in response to 25% PEG up-regulated \u003cem\u003eCAD\u003c/em\u003e2, \u003cem\u003eATH12\u003c/em\u003e and \u003cem\u003eCAD\u003c/em\u003e2, \u003cem\u003eCCR\u003c/em\u003e, respectively. Despite significant regulation of some other genes, amount of expression changes was not more than two-fold changes. NC100 as a higher lignin accumulator down-regulated all other functional (\u003cem\u003eCAD6\u003c/em\u003e, \u003cem\u003eCCoAOMT\u003c/em\u003e, \u003cem\u003eCOMT\u003c/em\u003e, \u003cem\u003eF5H\u003c/em\u003e, and \u003cem\u003eCCR\u003c/em\u003e) and regulatory (\u003cem\u003eODO1\u003c/em\u003e) genes in response to drought stress. On the other hand, \u003cem\u003eCCR\u003c/em\u003e which plays an essential role in lignin biosynthesis in plants (Cui et al., \u003cspan citationid=\"CR14\" class=\"CitationRef\"\u003e2022\u003c/span\u003e; Goujon et al., \u003cspan citationid=\"CR21\" class=\"CitationRef\"\u003e2003\u003c/span\u003e) showed the highest upregulation in Burly while NC100 significantly decreased expression of \u003cem\u003eCCR\u003c/em\u003e. As we know \u003cem\u003eCAD\u003c/em\u003e increased the synthesis of cinnamyl alcohols and is thought to be a specific lignification marker (Mitchell et al., \u003cspan citationid=\"CR38\" class=\"CitationRef\"\u003e1994\u003c/span\u003e). Previous studies reported that lignin content had a positive correlation with \u003cem\u003eCAD\u003c/em\u003e expression (Georgii et al., \u003cspan citationid=\"CR20\" class=\"CitationRef\"\u003e2017\u003c/span\u003e; Hu et al., \u003cspan citationid=\"CR26\" class=\"CitationRef\"\u003e2009\u003c/span\u003e; Jin et al., \u003cspan citationid=\"CR27\" class=\"CitationRef\"\u003e2014\u003c/span\u003e; Liu et al., \u003cspan citationid=\"CR36\" class=\"CitationRef\"\u003e2020\u003c/span\u003e). As reported in Table\u0026nbsp;3, meta-analysis demonstrated that all of the lignin-related gene isoforms differentially up- and down-regulated in response to drought stress with the exception of \u003cem\u003eCCR\u003c/em\u003e, \u003cem\u003eF5H\u003c/em\u003e, and \u003cem\u003eCCoAOMT\u003c/em\u003e.\u003c/p\u003e \u003cp\u003eThis study presents a major conclusion by challenging the prevalent assumption of meta-analysis as an acceptable method for identifying important genes inside complex pathways including the lignin pathway. By exploring into the intricate nature of these pathways and revealing the scarcity of individual studies available for meta-analysis, the study highlights the limits and emphasizes the importance of taking a more complete approach to gene discovery. This innovative concept underlines the importance of undertaking more individual investigations to improve the reliability and comprehensiveness of gene identification within complicated metabolite pathways.\u003c/p\u003e"},{"header":"Declarations","content":"\u003cp\u003e\u003cstrong\u003eAcknowledgments\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors would like to thank the Plant Biotechnology Department in Shahid Beheshti University and Institute of Biotechnology for supporting this research and the Bioinformatics Research Group in the College of Agriculture (Shiraz University).\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFunding\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors received no financial support for the research, authorship, and publication of this article.\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eCompeting Interests\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eThe authors declare no conflict of interest\u0026nbsp;\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eAuthors\u0026apos; contributions\u003c/strong\u003e\u003c/p\u003e\n\u003cp\u003eMaryam Rashidifar performing the work, contributed to data acquisition and manuscript drafting. Hossein Askari design and supervised the research. Ali Moghadam contributed to design and data acquisition. All authors read and approved the final manuscript.\u003c/p\u003e"},{"header":"References","content":"\u003col\u003e\n\u003cli\u003eAnders S., Pyl P.T., Huber W. (2015) HTSeq--a Python framework to work with high-throughput sequencing data. Bioinformatics 31:166-9. DOI: 10.1093/bioinformatics/btu638.\u003c/li\u003e\n\u003cli\u003eBaghel S., Sahariah B.P., Anandkumar J. (2020) Bioremediation of Lignin-Rich Pulp and Paper Industry Effluent, in: M. Shah and A. Banerjee (Eds.), Combined Application of Physico-Chemical \u0026amp; Microbiological Processes for Industrial Effluent Treatment Plant, Springer Singapore, Singapore. pp. 261-278.\u003c/li\u003e\n\u003cli\u003eBailey T.L., Williams N., Misleh C., Li W.W. (2006) MEME: discovering and analyzing DNA and protein sequence motifs. Nucleic acids research 34:W369-W373.\u003c/li\u003e\n\u003cli\u003eBang S.W., Choi S., Jin X., Jung S.E., Choi J.W., Seo J.S., Kim J.-K. (2022) Transcriptional activation of rice CINNAMOYL-CoA REDUCTASE 10 by OsNAC5, contributes to drought tolerance by modulating lignin accumulation in roots. Plant Biotechnology Journal 20:736-747. DOI: https://doi.org/10.1111/pbi.13752.\u003c/li\u003e\n\u003cli\u003eBang S.W., Lee D.-K., Jung H., Chung P.J., Kim Y.S., Choi Y.D., Suh J.-W., Kim J.-K. (2019) Overexpression of OsTF1L, a rice HD-Zip transcription factor, promotes lignin biosynthesis and stomatal closure that improves drought tolerance. Plant Biotechnology Journal 17:118-131. DOI: https://doi.org/10.1111/pbi.12951.\u003c/li\u003e\n\u003cli\u003eBarros J., Serk H., Granlund I., Pesquet E. (2015) The cell biology of lignification in higher plants. Ann Bot 115:1053-74. DOI: 10.1093/aob/mcv046.\u003c/li\u003e\n\u003cli\u003eBoerjan W., Ralph J., Baucher M. (2003) Lignin biosynthesis. Annual review of plant biology 54:519-546.\u003c/li\u003e\n\u003cli\u003eBonawitz N.D., Chapple C. (2010) The genetics of lignin biosynthesis: connecting genotype to phenotype. Annu Rev Genet 44:337-63. DOI: 10.1146/annurev-genet-102209-163508.\u003c/li\u003e\n\u003cli\u003eBroda M., Yelle D.J., Serwańska K. (2022) Bioethanol Production from Lignocellulosic Biomass-Challenges and Solutions. Molecules 27. DOI: 10.3390/molecules27248717.\u003c/li\u003e\n\u003cli\u003eCesarino I. (2019) Structural features and regulation of lignin deposited upon biotic and abiotic stresses. Curr Opin Biotechnol 56:209-214. DOI: 10.1016/j.copbio.2018.12.012.\u003c/li\u003e\n\u003cli\u003eChin C.-H., Chen S.-H., Wu H.-H., Ho C.-W., Ko M.-T., Lin C.-Y. (2014) cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC systems biology 8:1-7.\u003c/li\u003e\n\u003cli\u003eChoi S.J., Lee Z., Kim S., Jeong E., Shim J.S. (2023) Modulation of lignin biosynthesis for drought tolerance in plants. Frontiers in Plant Science 14. DOI: 10.3389/fpls.2023.1116426.\u003c/li\u003e\n\u003cli\u003eChun H.J., Lim L.H., Cheong M.S., Baek D., Park M.S., Cho H.M., Lee S.H., Jin B.J., No D.H., Cha Y.J., Lee Y.B., Hong J.C., Yun D.J., Kim M.C. (2021) Arabidopsis CCoAOMT1 Plays a Role in Drought Stress Response via ROS- and ABA-Dependent Manners. Plants (Basel) 10. DOI: 10.3390/plants10050831.\u003c/li\u003e\n\u003cli\u003eCui W., Zhuang Z., Jiang P., Pan J., Zhao G., Xu S., Shen W. (2022) Characterization, Expression Profiling, and Biochemical Analyses of the Cinnamoyl-CoA Reductase Gene Family for Lignin Synthesis in Alfalfa Plants. International Journal of Molecular Sciences 23:7762.\u003c/li\u003e\n\u003cli\u003eDauwe R., Morreel K., Goeminne G., Gielen B., Rohde A., Van Beeumen J., Ralph J., Boudet A.M., Kopka J., Rochange S.F., Halpin C., Messens E., Boerjan W. (2007) Molecular phenotyping of lignin-modified tobacco reveals associated changes in cell-wall metabolism, primary metabolism, stress metabolism and photorespiration. Plant J 52:263-85. DOI: 10.1111/j.1365-313X.2007.03233.x.\u003c/li\u003e\n\u003cli\u003eDidi V., Jackson P., Hej\u0026aacute;tko J. (2015) Hormonal regulation of secondary cell wall formation. Journal of Experimental Botany 66:5015-5027. DOI: 10.1093/jxb/erv222.\u003c/li\u003e\n\u003cli\u003eDiehl B.G., Brown N.R. (2014) Lignin cross-links with cysteine- and tyrosine-containing peptides under biomimetic conditions. J Agric Food Chem 62:10312-9. DOI: 10.1021/jf503897n.\u003c/li\u003e\n\u003cli\u003edos Santos A.B., Bottcher A., Kiyota E., Mayer J.L., Vicentini R., Brito Mdos S., Creste S., Landell M.G., Mazzafera P. (2015) Water stress alters lignin content and related gene expression in two sugarcane genotypes. J Agric Food Chem 63:4708-20. DOI: 10.1021/jf5061858.\u003c/li\u003e\n\u003cli\u003eFranke R., McMichael C.M., Meyer K., Shirley A.M., Cusumano J.C., Chapple C. (2000) Modified lignin in tobacco and poplar plants over‐expressing the Arabidopsis gene encoding ferulate 5‐hydroxylase. The Plant Journal 22:223-234.\u003c/li\u003e\n\u003cli\u003eGeorgii E., Jin M., Zhao J., Kanawati B., Schmitt-Kopplin P., Albert A., Winkler J.B., Sch\u0026auml;ffner A.R. (2017) Relationships between drought, heat and air humidity responses revealed by transcriptome-metabolome co-analysis. BMC Plant Biology 17:120. DOI: 10.1186/s12870-017-1062-y.\u003c/li\u003e\n\u003cli\u003eGoujon T., Ferret V., Mila I., Pollet B., Ruel K., Burlat V., Joseleau J.P., Barri\u0026egrave;re Y., Lapierre C., Jouanin L. (2003) Down-regulation of the AtCCR1 gene in Arabidopsis thaliana: effects on phenotype, lignins and cell wall degradability. Planta 217:218-28. DOI: 10.1007/s00425-003-0987-6.\u003c/li\u003e\n\u003cli\u003eGraham M.Y., Graham T.L. (1991) Rapid accumulation of anionic peroxidases and phenolic polymers in soybean cotyledon tissues following treatment with Phytophthora megasperma f. sp. glycinea wall glucan. Plant Physiology 97:1445-1455.\u003c/li\u003e\n\u003cli\u003eGu H., Wang Y., Xie H., Qiu C., Zhang S., Xiao J., Li H., Chen L., Li X., Ding Z. (2020) Drought stress triggers proteomic changes involving lignin, flavonoids and fatty acids in tea plants. Scientific Reports 10:15504. DOI: 10.1038/s41598-020-72596-1.\u003c/li\u003e\n\u003cli\u003eGui J., Shen J., Li L. (2011) Functional Characterization of Evolutionarily Divergent 4-Coumarate:Coenzyme A Ligases in Rice Plant Physiology 157:574-586. DOI: 10.1104/pp.111.178301.\u003c/li\u003e\n\u003cli\u003eHoshi T., Heinemann S. (2001) Regulation of cell function by methionine oxidation and reduction. The Journal of physiology 531:1-11. DOI: 10.1111/j.1469-7793.2001.0001j.x.\u003c/li\u003e\n\u003cli\u003eHu Y., Li W.C., Xu Y., Li G., Liao Y., Fu F.-L. (2009) Differential expression of candidate genes for lignin biosynthesis under drought stress in maize leaves. Journal of applied genetics 50:213-223.\u003c/li\u003e\n\u003cli\u003eJin Y., Zhang C., Liu W., Qi H., Chen H., Cao S. (2014) The cinnamyl alcohol dehydrogenase gene family in melon (Cucumis melo L.): bioinformatic analysis and expression patterns. PLoS One 9:e101730. DOI: 10.1371/journal.pone.0101730.\u003c/li\u003e\n\u003cli\u003eKanehisa M., Sato Y., Kawashima M., Furumichi M., Tanabe M. (2016) KEGG as a reference resource for gene and protein annotation. Nucleic acids research 44:D457-D462.\u003c/li\u003e\n\u003cli\u003eKhadr A., Wang Y.H., Zhang R.R., Wang X.R., Xu Z.S., Xiong A.S. (2020) Cytokinin (6-benzylaminopurine) elevates lignification and the expression of genes involved in lignin biosynthesis of carrot. Protoplasma 257:1507-1517. DOI: 10.1007/s00709-020-01527-8.\u003c/li\u003e\n\u003cli\u003eKim D., Langmead B., Salzberg S.L. (2015) HISAT: a fast spliced aligner with low memory requirements. Nature Methods 12:357-360. DOI: 10.1038/nmeth.3317.\u003c/li\u003e\n\u003cli\u003eKoutaniemi S., Warinowski T., K\u0026auml;rk\u0026ouml;nen A., Alatalo E., Fossdal C.G., Saranp\u0026auml;\u0026auml; P., Laakso T., Fagerstedt K.V., Simola L.K., Paulin L., Rudd S., Teeri T.H. (2007) Expression profiling of the lignin biosynthetic pathway in Norway spruce using EST sequencing and real-time RT-PCR. Plant Mol Biol 65:311-28. DOI: 10.1007/s11103-007-9220-5.\u003c/li\u003e\n\u003cli\u003eLangfelder P., Horvath S. (2008) WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics 9:559. DOI: 10.1186/1471-2105-9-559.\u003c/li\u003e\n\u003cli\u003eLee D.K., Jung H., Jang G., Jeong J.S., Kim Y.S., Ha S.H., Do Choi Y., Kim J.K. (2016) Overexpression of the OsERF71 Transcription Factor Alters Rice Root Structure and Drought Resistance. Plant Physiol 172:575-88. DOI: 10.1104/pp.16.00379.\u003c/li\u003e\n\u003cli\u003eLiao J., Lu Z., Hu S., Li Q., Che L., Chen X.D. (2017) Effects of prewash on the pyrolysis kinetics of cut tobacco. Drying Technology 35:1368-1378. DOI: 10.1080/07373937.2017.1320803.\u003c/li\u003e\n\u003cli\u003eLiu C., Yu H., Rao X., Li L., Dixon R.A. (2021) Abscisic acid regulates secondary cell-wall formation and lignin deposition in \u0026lt;em\u0026gt;Arabidopsis thaliana\u0026lt;/em\u0026gt; through phosphorylation of NST1. Proceedings of the National Academy of Sciences 118:e2010911118. DOI: 10.1073/pnas.2010911118.\u003c/li\u003e\n\u003cli\u003eLiu W., Jiang Y., Wang C., Zhao L., Jin Y., Xing Q., Li M., Lv T., Qi H. (2020) Lignin synthesized by CmCAD2 and CmCAD3 in oriental melon (Cucumis melo L.) seedlings contributes to drought tolerance. Plant Mol Biol 103:689-704. DOI: 10.1007/s11103-020-01018-7.\u003c/li\u003e\n\u003cli\u003eLivak K.J., Schmittgen T.D. (2001) Analysis of relative gene expression data using real-time quantitative PCR and the 2\u0026minus; \u0026Delta;\u0026Delta;CT method. methods 25:402-408.\u003c/li\u003e\n\u003cli\u003eMitchell H.J., Hall J.L., Barber M.S. (1994) Elicitor-Induced Cinnamyl Alcohol Dehydrogenase Activity in Lignifying Wheat (Triticum aestivum L.) Leaves. Plant Physiol 104:551-556. DOI: 10.1104/pp.104.2.551.\u003c/li\u003e\n\u003cli\u003eMoriya Y., Itoh M., Okuda S., Yoshizawa A.C., Kanehisa M. (2007) KAAS: an automatic genome annotation and pathway reconstruction server. Nucleic Acids Res 35:W182-5. DOI: 10.1093/nar/gkm321.\u003c/li\u003e\n\u003cli\u003eMurashige T., Skoog F. (1962) A Revised Medium for Rapid Growth and Bio Assays with Tobacco Tissue Cultures. Physiologia Plantarum 15:473-497. DOI: https://doi.org/10.1111/j.1399-3054.1962.tb08052.x.\u003c/li\u003e\n\u003cli\u003eNakano Y., Yamaguchi M., Endo H., Rejab N.A., Ohtani M. (2015) NAC-MYB-based transcriptional regulation of secondary cell wall biosynthesis in land plants. Frontiers in Plant Science 6. DOI: 10.3389/fpls.2015.00288.\u003c/li\u003e\n\u003cli\u003eSibout R., Eudes A., Mouille G., Pollet B., Lapierre C., Jouanin L., S\u0026eacute;guin A. (2005) CINNAMYL ALCOHOL DEHYDROGENASE-C and -D are the primary genes involved in lignin biosynthesis in the floral stem of Arabidopsis. Plant Cell 17:2059-76. DOI: 10.1105/tpc.105.030767.\u003c/li\u003e\n\u003cli\u003eSibout R., Proost S., Hansen B.O., Vaid N., Giorgi F.M., Ho-Yue-Kuang S., Leg\u0026eacute;e F., C\u0026eacute;zart L., Bouchabk\u0026eacute;-Coussa O., Soulhat C., Provart N., Pasha A., Le Bris P., Roujol D., Hofte H., Jamet E., Lapierre C., Persson S., Mutwil M. (2017) Expression atlas and comparative coexpression network analyses reveal important genes involved in the formation of lignified cell wall in Brachypodium distachyon. New Phytol 215:1009-1025. DOI: 10.1111/nph.14635.\u003c/li\u003e\n\u003cli\u003eSmith C.J., Hansch C. (2000) The relative toxicity of compounds in mainstream cigarette smoke condensate. Food Chem Toxicol 38:637-46. DOI: 10.1016/s0278-6915(00)00051-x.\u003c/li\u003e\n\u003cli\u003eSong Z., Wang D., Gao Y., Li C., Jiang H., Zhu X., Zhang H. (2021) Changes of lignin biosynthesis in tobacco leaves during maturation. Functional Plant Biology 48:624-633. DOI: https://doi.org/10.1071/FP20244.\u003c/li\u003e\n\u003cli\u003eSun S.-C., Xiong X.-P., Zhang X.-L., Feng H.-J., Zhu Q.-H., Sun J., Li Y.-J. (2020) Characterization of the Gh4CL gene family reveals a role of Gh4CL7 in drought tolerance. BMC Plant Biology 20:125. DOI: 10.1186/s12870-020-2329-2.\u003c/li\u003e\n\u003cli\u003eTao J., Li S., Ye F., Zhou Y., Lei L., Zhao G. (2020) Lignin - An underutilized, renewable and valuable material for food industry. Crit Rev Food Sci Nutr 60:2011-2033. DOI: 10.1080/10408398.2019.1625025.\u003c/li\u003e\n\u003cli\u003eVan Acker R., D\u0026eacute;jardin A., Desmet S., Hoengenaert L., Vanholme R., Morreel K., Laurans F., Kim H., Santoro N., Foster C., Goeminne G., L\u0026eacute;g\u0026eacute;e F., Lapierre C., Pilate G., Ralph J., Boerjan W. (2017) Different Routes for Conifer- and Sinapaldehyde and Higher Saccharification upon Deficiency in the Dehydrogenase CAD1 Plant Physiology 175:1018-1039. DOI: 10.1104/pp.17.00834.\u003c/li\u003e\n\u003cli\u003eVan Acker R., Vanholme R., Storme V., Mortimer J.C., Dupree P., Boerjan W. (2013) Lignin biosynthesis perturbations affect secondary cell wall composition and saccharification yield in Arabidopsis thaliana. Biotechnology for Biofuels 6:46. DOI: 10.1186/1754-6834-6-46.\u003c/li\u003e\n\u003cli\u003eVanholme R., Morreel K., Ralph J., Boerjan W. (2008) Lignin engineering. Current Opinion in Plant Biology 11:278-285. DOI: https://doi.org/10.1016/j.pbi.2008.03.005.\u003c/li\u003e\n\u003cli\u003eWang H., Avci U., Nakashima J., Hahn M.G., Chen F., Dixon R.A. (2010) Mutation of WRKY transcription factors initiates pith secondary wall formation and increases stem biomass in dicotyledonous plants. Proceedings of the National Academy of Sciences 107:22338-22343.\u003c/li\u003e\n\u003cli\u003eWang X., Kang D.D., Shen K., Song C., Lu S., Chang L.-C., Liao S.G., Huo Z., Tang S., Ding Y., Kaminski N., Sibille E., Lin Y., Li J., Tseng G.C. (2012) An R package suite for microarray meta-analysis in quality control, differentially expressed gene analysis and pathway enrichment detection. Bioinformatics (Oxford, England) 28:2534-2536. DOI: 10.1093/bioinformatics/bts485.\u003c/li\u003e\n\u003cli\u003eYang Q., He Y., Kabahuma M., Chaya T., Kelly A., Borrego E., Bian Y., El Kasmi F., Yang L., Teixeira P., Kolkman J., Nelson R., Kolomiets M., L Dangl J., Wisser R., Caplan J., Li X., Lauter N., Balint-Kurti P. (2017) A gene encoding maize caffeoyl-CoA O-methyltransferase confers quantitative resistance to multiple pathogens. Nature Genetics 49:1364-1372. DOI: 10.1038/ng.3919.\u003c/li\u003e\n\u003cli\u003eZeng Y., Himmel M.E., Ding S.-Y. (2017) Visualizing chemical functionality in plant cell walls. Biotechnology for Biofuels 10:263. DOI: 10.1186/s13068-017-0953-3.\u003c/li\u003e\n\u003cli\u003eZhao D., Luan Y., Shi W., Zhang X., Meng J., Tao J. (2021) A Paeonia ostii caffeoyl-CoA O-methyltransferase confers drought stress tolerance by promoting lignin synthesis and ROS scavenging. Plant Sci 303:110765. DOI: 10.1016/j.plantsci.2020.110765.\u003c/li\u003e\n\u003cli\u003eZhao J., Zou W., Hu T. (2020) Novel genes associated with folic acid-mediated metabolism in mouse: A bioinformatics study. PloS one 15:e0238940.\u003c/li\u003e\n\u003cli\u003eZhao Q., Wang H., Yin Y., Xu Y., Chen F., Dixon R.A. (2010) Syringyl lignin biosynthesis is directly regulated by a secondary cell wall master switch. Proceedings of the National Academy of Sciences 107:14496-14501.\u003c/li\u003e\n\u003cli\u003eZheng Y., Jiao C., Sun H., Rosli H.G., Pombo M.A., Zhang P., Banf M., Dai X., Martin G.B., Giovannoni J.J. (2016) iTAK: a program for genome-wide prediction and classification of plant transcription factors, transcriptional regulators, and protein kinases. Molecular plant 9:1667-1670.\u003c/li\u003e\n\u003c/ol\u003e"},{"header":"Tables","content":"\u003cp\u003eTables 1 to 6 are available in the Supplementary Files section.\u003c/p\u003e"},{"header":"Supporting information","content":"\u003cp\u003eSupporting information is not available with this version. \u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFig S1\u003c/strong\u003e Weighted gene co-expression network analysis (WGCNA)\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFig S2\u003c/strong\u003e GO classification in 2 selected modules\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eFig S3\u003c/strong\u003e Transcription factor family classifications co-expressed throughout 2 selected modules\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTable S1\u003c/strong\u003e The enriched KEGG pathway of DEGs\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTable S2\u003c/strong\u003e sequences of conserved promoter motifs of up-regulated meta-genes\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTable S3\u003c/strong\u003e sequences of conserved promoter motifs of down-regulated meta-genes\u003c/p\u003e\n\u003cp\u003e\u003cstrong\u003eTable S4\u003c/strong\u003e List of modules in Weighted gene co-expression network analysis (WGCNA)\u003c/p\u003e"}],"fulltextSource":"","fullText":"","funders":[],"hasAdminPriorityOnWorkflow":false,"hasManuscriptDocX":true,"hasOptedInToPreprint":true,"hasPassedJournalQc":"","hasAnyPriority":false,"hideJournal":true,"highlight":"","institution":"","isAcceptedByJournal":false,"isAuthorSuppliedPdf":false,"isDeskRejected":"","isHiddenFromSearch":false,"isInQc":false,"isInWorkflow":false,"isPdf":false,"isPdfUpToDate":true,"isWithdrawnOrRetracted":false,"journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true},"keywords":"Lignin, Meta-analysis, RNA-seq, Tobacco, WGCNA","lastPublishedDoi":"10.21203/rs.3.rs-4101335/v1","lastPublishedDoiUrl":"https://doi.org/10.21203/rs.3.rs-4101335/v1","license":{"name":"CC BY 4.0","url":"https://creativecommons.org/licenses/by/4.0/"},"manuscriptAbstract":"\u003cp\u003eLignin plays a crucial impact on the production of phenolic compounds in tobacco smoke, which have potential health implications associated with tobacco use. The meta-analysis of RNA-seq studies along with high-resolution expression analysis on \u003cem\u003eNicotiana tabacum\u003c/em\u003e clarified a conserved distinctive expression pattern of lignin gene network. According to the results, 67 DEGs associated with lignin biosynthesis network were identified of which 17 genes were introduced by meta-analysis. WGCNA showed 14 clusters for the meta-genes. Various TF families and a number of regulatory factors were identified as the most likely candidate genes associated with the lignin metabolic pathway. Eight major meta-genes were evaluated by using qRT-PCR in two tobacco genotypes with different lignin content under drought stress conditions. Genotype NC100 (high lignin content) and Burly (low lignin content) in response to PEG upregulated \u003cem\u003eCAD\u003c/em\u003e2, \u003cem\u003eATH12\u003c/em\u003e and \u003cem\u003eCAD\u003c/em\u003e2, \u003cem\u003eCCR\u003c/em\u003e, respectively. Despite the accumulation of lignin, the expression levels of \u003cem\u003eCCoAOMT\u003c/em\u003e, \u003cem\u003eF5H\u003c/em\u003e, \u003cem\u003eCOMT\u003c/em\u003e, and \u003cem\u003eODO1\u003c/em\u003e were reduced in both genotypes. The study's exploration into the complex nature of these pathways and meta-analysis highlights the importance of adopting a more comprehensive approach to gene discovery. It suggests that conducting additional individual investigations is crucial for enhancing the reliability and comprehensiveness of gene identification within intricate metabolite pathways.\u003c/p\u003e","manuscriptTitle":"Integrated transcriptomic analysis reveals a transcriptional regulation network for the biosynthesis of lignin in Nicotiana tabacum in drought stress response","msid":"","msnumber":"","nonDraftVersions":[{"code":1,"date":"2024-03-19 16:07:03","doi":"10.21203/rs.3.rs-4101335/v1","editorialEvents":[{"type":"communityComments","content":0}],"status":"published","journal":{"display":true,"email":"[email protected]","identity":"researchsquare","isNatureJournal":false,"hasQc":true,"allowDirectSubmit":true,"externalIdentity":"","sideBox":"","snPcode":"","submissionUrl":"/submission","title":"Research Square","twitterHandle":"researchsquare","acdcEnabled":true,"dfaEnabled":false,"editorialSystem":"","reportingPortfolio":"","inReviewEnabled":false,"inReviewRevisionsEnabled":true}}],"origin":"","ownerIdentity":"e60c3449-cc42-4f2e-ad49-ffa366dca58b","owner":[],"postedDate":"March 19th, 2024","published":true,"recentEditorialEvents":[],"rejectedJournal":[],"revision":"","amendment":"","status":"posted","subjectAreas":[],"tags":[],"updatedAt":"2024-04-03T13:29:49+00:00","versionOfRecord":[],"versionCreatedAt":"2024-03-19 16:07:03","video":"","vorDoi":"","vorDoiUrl":"","workflowStages":[]},"version":"v1","identity":"rs-4101335","journalConfig":"researchsquare"},"__N_SSP":true},"page":"/article/[identity]/[[...version]]","query":{"redirect":"/article/rs-4101335","identity":"rs-4101335","version":["v1"]},"buildId":"qtupq5eGEP_6zYnWcrvyt","isFallback":false,"isExperimentalCompile":false,"dynamicIds":[84888],"gssp":true,"scriptLoader":[]}