Cavendish banana (Musa acuminata L. AAA group cv. Cavendish) fruits were harvested from the banana plantation of the Institute of Bioscience and Biotechnology (ITBB) of Chinese Academy of Tropical Agricultural Sciences (CATAS) (Wenchang, Hainan, 19.61N, 110.75E). The pulp of Cavendish banana fruit was harvested at 0, 20, and 80 DAF (day after flowering, DAF), which represented fruit developmental stages of budding, cutting flower, and harvest stages, respectively. During ripening stage, pulp of Cavendish banana fruit was harvested at 0, 8 and 14 DPH (day post-harvest, DPH) representing green, yellowish green, and yellow based on color of the fruit, respectively.

One month old Cavendish banana plantlets were grown under greenhouse conditions with 26°C and 80% humidity. Cavendish banana plantlets were irrigated with 300 mM NaCL and 200 mM mannitol for 7 days for salt and osmotic treatments, respectively. Cavendish banana plants were maintained at 4°C for 22 h for cold treatment. The leaves without main vein were harvested for analysis.

Fusarium Wilt of Banana (FWB) is prevalent in main banana producing areas worldwide, which is seriously destructive for banana industry. It is a typical soil borne fungus disease caused by Fusarium oxysporum f. sp. cubense (Foc). Cavendish banana is highly susceptible to Foc TR4. Giant Cavendish Tissue Culture Variants (GCTCV, GC) have acquired resistance to TR4 through somaclonal variation (Hwang and Ko, 2004). GCTCV-119 is the best Foc TR4-resistant alternative cultivar for the Cavendish (Ploetz, 2015). Variety Pahang (PH) is a wild banana belonging to subspecies Musa acuminata ssp. malaccensis, and highly resistant to Foc TR4 (Zhang et al., 2018). The fungus was grown on potato dextrose agar medium and incubated for 5-7 days at 28°C. The roots of 2-month old susceptible cultivar Cavendish banana, resistant cultivar GCTCV-119 and variety Pahang were dipped in Foc TR4 spore suspension of 1.0×10⁶ condia/mL, and the roots were harvested at 0 and 2 days post-infection (DPI) (Wang et al., 2012). All of the above samples were immediately frozen in liquid nitrogen and stored at -80°C.

To analyze the genome-wide lignin biosynthesis genes, we first downloaded all gene coding protein sequences of Musa acuminata genome (DH Pahang) from the Banana Genome Hub (https://banana-genome-hub.southgreen.fr) (Martin et al., 2016). We compared genes related to lignin biosynthesis genes in the Musa acuminata genome with the genes annotated in genomes of rice (http://rice.uga.edu/index.shtml) and Arabidopsis (http://www.arabidopsis.org/). In the PAL and HCT gene families, we also selected representative genes to draw the phylogenetic tree, including EsPAL (ATD50344.1) from Eucalyptus saligna, JcPAL (XP_012077001.1) from Jatropha curcas, HbPAL (XP_021684949.1) from Hevea brasiliensis, PtHCT1 (XP_024452069.1) from Populus trichocarpa, PtHCT2 (XP_006368492.1), NtHCT (NP_001312552.1) and NtHQT(NP_001312079.1) from Nicotiana tabacum, and CcHCT (ABO47805.1) from Coffea canephora. We retrieved protein sequences of these gene families for homologue-based searches with the criteria: similarity > 80% and coverage > 80%. Then, we confirmed the presence of the conserved domain within all the protein sequences in CDD (http://www.ncbi.nlm.nih.gov/cdd/) databases in NCBI and removed members without complete domain. The full-length lignin biosynthesis genes protein sequences from Musa acuminata, rice, Arabidopsis and other woody plants were aligned using ClustalW. Relationships were assessed using a Maximum Likelihood evolutionary tree with 1,000 bootstrap replicates and were created using MEGA 6.0 software (Tamura et al., 2013). The accession number of identified lignin biosynthesis genes was listed in Supplementary Table 1 . The molecular weight and isoelectric points of the lignin biosynthesis genes were predicted from the ExPASy database (http://expasy.org/). The sequence logo for lignin biosynthesis genes domain was created by WebLogo server (http://weblogo.berkeley.edu/logo.cgi).

To determine the physical locations of lignin biosynthesis genes, the starting and ending positions of all lignin biosynthesis genes on each chromosome were obtained from the Musa acuminata genome database. MapInspect software was used to draw the images of the locations of the lignin biosynthesis genes (http://mapinspect.software.informer.com/). Tandem and segmental duplications were also identified according to the plant genome duplication database (Lee et al., 2012). Based on physical chromosomal location, homologous lignin biosynthesis genes on a single chromosome within 30 kb of each other were characterized as tandem duplication (Shiu and Bleecker, 2003; Du et al., 2013). Syntenic blocks were detected using MCSCAN (parameters: -a -e 1e-5 -s 5) (Wang et al., 2019), and all lignin biosynthesis genes located in the syntenic blocks were extracted. Circos (0.63) software was used to draw the images of the locations and synteny of the lignin biosynthesis genes (http://circos.ca/).

Total RNAs were extracted using Eastep Super Total RNA Extraction Kit (Promega, Beijing, China). Five μg of total RNA from each sample was converted to cDNA using GoScript Reverse Transcription Mix (Promega, Beijing, China), and constructed the cDNA libraries using TruSeq RNA Library Preparation Kit v2, and then were sequenced on an IIIumina HiSeq 2000 platform (San Diego, CA, USA) using the IIIumina RNA-seq protocol. Two biological replicates were used for each sample. Gene expression levels were calculated as Reads Per Kilobases per Million reads (RPKM) (Mortazavi et al., 2008). Differentially expressed genes were identified with the mean of read count from two replicates for each gene (fold change ≥ 2) (Audic and Claverie, 1997). MeV 4.9 and Java Treeview software constructed the heat-map of genes according to the manufacturer’s protocol. We analyzed the expression patterns of lignin biosynthesis genes during fruit development and postharvest ripening (0DAF, 20DAF, 80DAF-0DPH, 8DPH and 14DPH), and the response to abiotic (4°C for 22h, 300 mM NaCl for 7 days and 200 mM mannitol for 7 days) and biological (Roots of Cavendish banana, GCTCV-119 and Pahang varieties inoculated Foc TR4 after 2 days) stressors. All the data of RNA-seq had been uploaded to deposit in the CNSA (https://db.cngb.org/cnsa/) of CNGBdb with accession number CNP0000292. The accession numbers of all samples were listed in Supplementary Table 2 .

Gene expression patterns for all identified TFs and the lignin biosynthesis genes from fruits, leaves and roots tissues were used to construct the co-expression network using WGCNA (version 1.47) within the R platform (version 3.2.2) (Langfelder and Horvath, 2008). Genes whose expression had not been detected in all tissues were removed in advance. Soft thresholds were set according to the scale-free topology criterion adopted by Zhang and Horvath (2005). An adjacency matrix was developed using the squared Euclidean distance values, and Pearson method was used to calculate the topological overlap matrix of unsigned network detection. Then, the co-expression coefficients greater than 0.50 between the target genes were selected. Finally, we extracted the co-expression network of lignin biosynthesis genes, and the network connections were visualized using cytoscape (Shannon et al., 2003). We enriched analysis of target genes using the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.

Yeast one-hybrid screening was performed using MatchmakerTM Gold Yeast One-Hybrid Library Screening System (Clontech, Dalian, China). The bait fragment (the 1481 bp fragment of MaC3H2 promoter) was cloned into the pAbAi vector. The MaC3H2-AbAi and p53-AbAi were linearized and transformed into Y1HGold to make a bait-reporter strain, respectively. Transformants were initially screened on plates containing SD medium without Ura (SD/-Ura) and added with 0-1000 ng/ml aureobasidin A (AbA) for auto-activation analysis. Full-length coding sequences of MaMYB1 (Ma02_t00290.1), MaMYB2 (Ma06_t05680.1), MaMYB3 (Ma11_t11940.1), MaNAC1 (Ma07_t06080.1) and MabHLH1 (Ma07_t28500.1) were cloned into the pGADT7 (AD) prey vector and transferred into the bait-reporter yeast strain, respectively. Transformed Y1HGold were incubated on SD medium with 300 ng/ml AbA and without leucine (SD-Leu+AbA300) at 28°C for 3 days to test the interaction. pGADT7-Rec (AD-Rec-P53) and p53-promoter fragments were co-transformed into Y1HGold as positive controls, while the AD-empty and MaC3H2-AbAi were used as negative controls. The primers used for yeast one-hybrid assay were listed in Supplementary Table 3 .

Full-length coding sequences of the MaMYB1, MaMYB2, MaMYB3, MaNAC1 and MabHLH1 were inserted into the pGreen II 62-SK vector (SK), respectively, and the fragment of MaC3H2 promoter was inserted into the pGreen II 0800-LUC vector. All the above constructs were transformed into Agrobacterium tumefaciens GV3101 using freeze-thaw method (Holsters et al., 1978). The dual luciferase assays were performed with Nicotiana benthamiana grown in green house for forty days. Agrobacterium cultures were prepared with infiltration buffer (10 mM MgCl₂, 10 mM MES and 150 mM acetosyringone, pH5.6) to an OD₆₀₀ of 0.8. Agrobacterium culture mixtures of TF genes (5 ml) and promoter (95 μL) were infiltrated into the abaxial side of Nicotiana benthamiana leaves. Leaves were collected after two days infiltration for luciferase (LUC) and Renilla luciferase (REN) activities analyses using Dual-Luciferase Reporter Assay System (Promega, USA, E1910) with Modulus Luminometers (Promega). Luciferase activity was analyzed in three independent experiments with six replications for each assay.

Statistical analysis was performed using Student’s t-test. The experimental results obtained were expressed as the means ± standard deviation (SD). P values<0.05 were considered statistically significant (*), and P values<0.01 were considered highly statistically significant (**).

In total, 104 genes belonging to ten key enzyme gene families, including eight MaPALs, five MaC4Hs, eighteen Ma4CLs, three MaHCTs, two MaC3Hs, seven MaCCoAOMTs, twenty-four MaCCRs, four MaF5Hs, twenty-four MaCOMTs and ten MaCADs were identified in the Musa acuminata genome. In addition, we found 31 genes identified as ‘like’ genes (candidate genes excluded from bona fide clade), including 9 from Ma4CLs, 21 from MaCCRs and 1 from MaCOMTs ( Supplementary Table 1 ). Gene names were based on the location of chromosomes. The gene ID, position, and protein physicochemical properties of all genes were shown in Supplementary Table 1 .

In conservative domain analysis, we mainly compared the lignin biosynthesis genes from banana with genes from Arabidopsis. We found all genes contained the core conserved domains belonging to the features of gene family in the protein sequence. MaPALs contained the conserved motifs (GTITASGDLUPLSYIAG) ( Figure 1A ). MaC3Hs, MaC4Hs and MaF5Hs had cytochrome P450 conserved domains, such as proline rich region (PPGPKGLP), PERF conserved motif and heme-binding motif of CYP450 proteins (PFGSGRRSCP) ( Figure 1B ). In Ma4CLs, all of them contained conservative Box I (SSGTTGLPKGV) and Box II (GEICVRS) in the protein sequence ( Figure 1C ). We found MaHCTs genes contained the conserved active motif (HXXXDG) and conservative motif (DFGWGR) of acetyltransferase gene family ( Figure 1D ). Generally, the amino acid sequence of CCoAOMT gene family contains 8 conserved motifs, in which motif D, motif E, motif F, motif G and motif H are unique to the CCoAOMT family and are as the tag sequences (Zhao et al., 2004). We found all MaCCoAOMTs contained motif D, motif E, motif F, motif G and motif H ( Figure 1E ). In these sequences of COMT proteins, there were five amino acid conserved elements in proper order from the N-terminal to the C-terminal of the amino acid sequence (conservative motif 1: LVDGGGxG, conservative motif II: GINFDLPHV, conservative motif III: EHVGGDMF and conservative motif VI: NGKVI) (Ibrahim et al., 1998). In Musa acuminata genome, we found that 23 COMT genes had the above conserved domains ( Figure 1F ). Three bona fide MaCCRs (MaCCR1, MaCCR2 and MaCCR3) were identified to contain the conserved sequence (KNWYCYGK) belonging to CCR gene family ( Figure 1G ). Ten MaCAD proteins contained binding domain for the structurally important zinc (CX2CX2CX7C), NADPH (GxGGVG) and catalytic Zn²⁺ binding motif (GHExxGxxxxxGxxV), which was the signature sequence of CAD genes ( Figure 1H ). The above results show that all 104 lignin biosynthesis genes have typical conserved domains of their respective gene families.

Analyses of the evolutionary relationships of lignin synthetase genes from Musa acuminata, Arabidopsis, rice, and other plants have been functionally verified. The phylogenetic tree was analyzed by MEGA6 using the method with 1000 bootstrap replicates and the evolutionary distances were calculated using the Poisson correction method.

MaPAL family can be clearly divided into two groups. MaPAL1, MaPAL2, MaPAL3, MaPAL4, MaPAL5, MaPAL7 and MaPAL8 were clustered together with all PALs from Arabidopsis and rice. MaPAL6 was placed near the PALs from woody plants, such as EsPAL, JcPAL and HbPAL, indicating that the evolution of MaPALs is more conservative, MaPAL6 may have different origins from the other seven MaPALs ( Figure 2A ). In plants, 4CL can be divided into three clades, clade I is thought to be related to lignin synthesis, and clade II 4CL is thought to be related to flavonoid synthesis (Wei and Wang, 2004). The Ma4CLs can be divided into three clades, Ma4CL1, Ma4CL3, Ma4CL4, Ma4CL5, Ma4CL6, Ma4CL7 were located in clade I, Ma4CL2 and Ma4CL8 were located in clade II, and others were located in 4CL-like clade. Moreover, we also found that the Ma4CLs were often closer to Os4CLs and far away from At4CLs ( Figure 2B ). It is suggested that Ma4CL11, Ma4CL14, Ma4CL6 and Ma4CL7 may be related to lignin synthesis and formed after the differentiation of monocotyledons and dicotyledons. MaC4Hs can be divided into two clades: MaC4H, MaC4H3, MaC4H4, MaC4H5, MaC4H7, AtC4H, OsC4H1 and OsC4H2 which were located in same clade, while MaC4H1, MaC4H2, MaC4H6, OsC4H3 and OsC4H4 were located in the same clade ( Figure 2C ). MaC3Hs were more closely related to OsC3H than AtC3H1 and no member was close to AtC3H2 and AtC3H3 ( Figure 2C ). Those results indicate that the evolution of MaC4Hs and MaC3Hs may be formed after the differentiation of monocotyledons and dicotyledons. The members of MaF5Hs were clustered with the AtF5Hs, and far from OsF5Hs ( Figure 2C ), indicating that the differentiation of F5H gene is formed before the differentiation of monocotyledons and dicotyledons. Three MaHCTs were distributed in the same clade with OsHCT1 and OsHCT2, but far away from the HCTs from Arabidopsis thaliana, Populus trichocarpa (PtHCT1, PtHCT2), Nicotiana tabacum (NtHCT and NtHQT) and Coffea canephora (CcHCT), indicating that MaHCTs may be formed after differentiation of monocotyledons and dicotyledons ( Figure 2C ). CCoAOMTs can be mainly divided into two clades (Rakoczy et al., 2018). All seven MaCCoAOMTs were distributed in the same clade with OsCCoAOMT and AtCCoAOMT1 ( Figure 2D ). Similarly, all the MaCOMTs located with AtOMT1 and OsCOMT1 belonged to the Clade I, and only MaCOMT21 close to AtOMT-like11 and OsCOMT-like5 belonged to COMT-like ( Figure 2D ). Bona fide CCR clade and CCR-like clade were contained in CCR gene family. In particular, MaCCR1 and MaCCR2 were closely related to TaCCR1, ZmCCR1, AtCCR1 and AtCCR2. MaCCR3 was grouped with ZmCCR2, TaCCR2, and PvCCR2 ( Figure 2E ). Thirty-one CADs were clearly divided into three classes (class I, class II and class III). MaCAD1 and MaCAD2 were located in class I (bona fide CAD clade) with AtCAD4, AtCAD5 and OsCAD2. MaCAD5, MaCAD6, MaCAD7, MaCAD8, MaCAD9 and MaCAD10 were distributed in class II. MaCAD3 and MaCAD4 were located in class III with AtCAD1, OsCAD1 and OsCAD4 ( Figure 2F ).

In order to determine the distribution of lignin biosynthesis genes on the chromosome, the relative positions of 101 genes on the chromosome were identified according to their physical positions in the Musa acuminata genome database. Ma00_t01670.1, Ma00_t02580.1 and Ma00_t03560.1 were not anchored on the chromosome and were named MaCOMT22, MaCOMT23 and MaCAD10, respectively. As shown in Figure 3 , although each of the 11 banana chromosomes contained lignin biosynthesis genes, the distribution appeared to be uneven. There were 17 members on chromosome 4 from 7 different gene families (MaHCTs, Ma4CLs, MaCADs, MaC3Hs, MaCCRs, MaF5Hs and MaCOMTs). There were 14 members on chromosome 9 from six families (Ma4CLs, MaC4Hs, MaCADs, MaCCoAOMTs, MaCOMTs and MaPALs), while only three genes on chromosome 3 were from three families (Ma4CL6, MaCAD6 and MaC4H2). In the multi member gene families, eighteen Ma4CLs were widely distributed on the 10 chromosomes except for Chr10. Twenty-four MaCCRs were distributed on 6 chromosomes as follows Chr3, Chr4, Chr5, Chr6, Chr8 and Chr9. Most of MaCOMTs were only distributed on Chr2 and Chr9.

Segmental duplication, tandem duplication and retrotransposition are known to be key factors driving gene family expansion (Hurles, 2004; Freeling, 2009). In our study, 11 colinear genes pairs from nine gene families except for MaC3Hs were identified to be located in the synthetic blocks and belonged to segmental duplication ( Figure 3 and Supplementary Table 4 ). It is worth noting that there were four MaCADs (MaCAD1, MaCAD2, MaCAD3 and MaCAD4) with segmental duplication located in Block4_chr04, Block8_chr01 and Block8_chr04, indicating that the segmental duplication plays important roles in the expansion of MaCADs. We also found four tandem duplication gene clusters on Chr2, Chr8, Chr9 and Chr10. The number of gene clusters on chromosome 2 was the largest, containing 11 genes, and all of them were from MaCOMT gene family ( Figure 3 ). We did not find the member belonging to retrotransposition among the 104 lignin biosynthesis genes. Based on our analysis, we propose that the expansion of lignin biosynthesis genes may be mainly via segment and tandem duplication during the Musa acuminata genome evolutionary processes.

To analyze the expression profiles of lignin biosynthesis genes, RNA-Seq data were derived from the fruit development and ripening stages ( Supplementary Table 2 ). We deleted 56 lignin biosynthesis genes with RPMK values less than 5 in 0DAF, 20DAF, 80DAF_0DPH, 8DPH and 14DPH. There were 25 and 26 lignin biosynthesis genes highly expressed at 0 DAF and 20 DAF (RPKM > 20), contributing about 24.04% and 25.00% of the total, respectively. In 0DAF and 20 DAF, the expression levels of MaC3H2, MaC4H2, MaC4H5, MaCAD7, MaCCR-like17, MaCCR-like19, MaCCR-like2, MaCCR-like3, MaCCR-like4, MaCCR-like6, MaCCR-like7 and MaPAL3 were more than 50 (RPKM>50). However, in 80 DAF-0 DPH, 8 DPH and 14 DPH, only 14, 16 and 10 genes were highly expressed (RPKM>20), respectively ( Figure 4 ). The number of highly expressed genes during postharvest ripening was significantly less than that during fruit development. These results indicate that the lignin biosynthesis genes are expressed actively in the fruit development, which is closely related to the growth and development of banana fruit. It is worth noting that MaC3H gene family contained two copies (MaC3H1 and MaC3H2) ( Figure 3 ). During fruit development and ripening, we found that MaC3H1 was not expressed in 11 transcriptome data and checked by RT-PCR (results not shown), while MaC3H2 was constitutively expressed with high RPKM value (RPKM>40) ( Supplementary Table 5 ), showing that only MaC3H2 may have biological functions. We also found that the expression levels of Ma4CL8, MaPAL4, MaC4H5, Ma4CL-like6, MaCCoAOMT2, MaCAD3, MaHCT2 and MaCAD2 were increased gradually in 80DAF-0DPH, 8DPH and 14DPH, and some genes such as Ma4CL-like9, MaHCT2, MaF5H4 and MaC4H5 reached the highest level in 14DPH (RPKM>190) ( Figure 4 ). The above results suggest that these lignin biosynthesis genes may play roles in the lignin synthesis during the postharvest process of banana fruit and participate in the quality formation.

To analyze the expression profiles of lignin biosynthesis genes, RNA-Seq data were derived from Cavendish banana plantlets in response to osmotic, salt and cold treatments. To present the differentially expressed genes visually and exactly, we filtered out the genes with RPKM values less than 5 in the control or treatments. Cold treatment (4°C) caused the fifteen lignin biosynthesis genes were differentially expressed (log2 Ratio Cold/Control>1), among which eleven were up-regulated and four down-regulated ( Figure 5 and Supplementary Table 6 ). Treatment with drought (200 mM mannitol) caused the eighteen genes were found to be differentially expressed (log2 Ratio Osmotic/Control > 1), among which nine were up-regulated and nine down-regulated ( Figure 5 and Supplementary Table 6 ). Treatment with salt (300 mM NaCl) caused fourteen members were differentially expressed (log2 Ratio Salt/Control > 1), among which thirteen were up-regulated and one was down regulated ( Figure 5 and Supplementary Table 6 ). It is worth noting that MaPAL1, MaPAL3, MaPAL4, MaPAL7, MaC4H2 and MaCCR9 were up-regulated during the drought, salt, and cold stresses. We found that forty-one genes were significantly differentially expressed under cold, drought and salt stresses, but nine genes were only involved in one of the abiotic stresses ( Figure 5 ). Therefore, most of the genes were not differentially expressed in the responses to above stresses. For example, we did not find that bona fide MaCCRs was involved in the process of the responses to above abiotic stresses. Although several MaCCR-likes were found to be up-regulated under low temperature, drought and salt stresses, the function of MaCCR-like was not clear. Similarly, we also found that only one member of MaCADs was up-regulated under salt stress, other eight MaCADs were not differentially expressed during the responses to abiotic stresses. Moreover, in phylogenetic tree, the differentially expressed MaCAD4 was closely related to AtCAD1, which played a partial role in the synthesis of lignin in Arabidopsis (Eudes et al., 2006). MaCAD2 and MaCAD3, belonging to the group I (the bona fide CAD clade) ( Figure 2 ), were not differentially expressed under the drought, salt and cold stresses ( Figure 5 ). These results indicate that most of lignin synthesis genes are not involved in banana response to stresses, which may be one of the possible reasons for banana being more sensitive to abiotic stresses.

The expression pattern of a gene is generally correlated with its function; hence, we analyzed the expression patterns of the lignin biosynthesis genes in susceptible and resistant varieties responding to Foc TR4 infection. To present the differentially expressed genes visually and exactly, we filtered out the genes with RPKM values less than 5 in the 0 DPI or 2 DPI. In variety Cavendish, there were only 11 differentially expressed genes (marked with black stars), of which only one was up-regulated and 10 down-regulated ( Figure 6 and Supplementary Table 7 ). The results showed that lignin synthesis genes were not activated in susceptible varieties. In variety Pahang, there were 12 differentially expressed genes from seven gene families, all of which were up-regulated ( Figure 6 and Supplementary Table 7 ), indicating that lignin biosynthesis genes were activated in banana responding to Foc TR4 infection. In variety GCTCV-119, we found 20 differentially expressed genes, 17 of which were up-regulated and 3 down-regulated ( Figure 6 and Supplementary Table 7 ), indicating that most of the differentially expressed genes were activated. It is noteworthy that we found that the differentially expressed genes between Pahang and GCTCV-119 were not consistent, indicating that ploidy still has a great impact on the banana resistance to Foc TR4 infection. In general, our results indicate that these differentially expressed lignin biosynthesis genes are involved in banana resistance to Foc TR4 infection.

Based on RNA-seq data, 8 differentially expressed lignin biosynthesis genes were selected for RT-qPCR. The results showed that the expression patterns of the selected genes had the same trend and consistent results between RNA-seq data and RT-qPCR data ( Figure 7 ).

To explore the functions, 104 lignin biosynthesis genes were selected as ‘guide genes’ to seek co-expressed genes using an RNA-Seq dataset from 11 different transcriptomes, including fruit development and ripening stages, banana plantlets response to osmotic, salt, and cold treatment, and banana roots inoculated with Foc TR4 ( Supplementary Table 2 ). A total of 725 TFs were as the ‘target nodes’, whose expression patterns were closely correlated with 51 lignin biosynthesis genes from ten gene families, which were identified with weighted values larger than 0.5 ( Supplementary Table 6 ). Following visualization using Cytoscape (Otasek et al., 2019), the co-expression network of lignin biosynthesis genes was divided into four modules. Module I-IV (purple, red, light blue and yellow, respectively) contained 37, 3, 4, 1 and 1 lignin biosynthesis genes, respectively ( Figure 8 ). The correlated TFs were in 68 different type TF classes, such as MYB, ERF, basic helix-loop-helix (bHLH), WRKY, C2H2, NAC, and LOB ( Supplementary Table 8 ).

The -2000bp upstream sequence of MaC3H2 was obtained by RT-PCR based on the sequence of Musa acuminata genome, the 1481bp promoter sequence was obtained by transcriptional active site prediction, and the main regulatory elements of MaC3H2 promoter were predicted online by plantcare website. We found that there was conserved AC element (AC-I element, -1378bp) and MBS and MRE elements (-32 and -1066bp) bounding to MYB TF in the MaC3H2 promoter sequence ( Supplementary Table 3 ). The results show that MaC3H2 may be closely related to lignin synthesis in banana.

We extracted the co-expression network of MaC3H2 and found 36 transcription factors from 13 gene families were contained in the network, and MaMYBs had the largest number (10 MaMYBs) ( Figure 9A ). We further identified whether the MaC3H2 promoter had a relationship with these transcription factors. Five transcription factors including MaMYB1 (Ma03_t25780.1), MaMYB2 (Ma09_t20280.1), MaMYB3 (Ma11_t11940.1), MaNAC1 (Ma07_t06080.1) and MabHLH1 (Ma10_t14490.1) were verified by yeast one hybridization. Using the pGADT7-MaMYB1, pGADT7-MaMYB2, pGADT7-MaMYB3, pGADT7-MaNAC1, pGADT7-MabHLH1 vectors as effectors and using the pAbAi vector carrying MaC3H2 promoter as a reporter, the Y1HGold yeast co-transformed with MaC3H2 promoter-AbAi and pGADT7-MaMYB1, pGADT7-MaMYB2, pGADT7-MaMYB3, pGADT7-MaNAC1, pGADT7-MabHLH1 could grow normally on the SD/-Leu medium containing 300 ng/mL AbA, respectively, whereas the yeast carrying MaC3H2 promoter-AbAi and pGADT7 could not grow ( Figure 9B ), which proved that MaMYB1, MaMYB2, MaMYB3, MaNAC1 and MabHLH1 can definitely bind with the sequence of the MaC3H2 promoter, respectively.

In the dual luciferase assay, the MaC3H2 promoter-driven luciferase reporter and CaMV35S-driven MaMYB1, MaMYB2, MaMYB3, MaNAC1 and MabHLH1 effector were separately constructed and then co-infected tobacco leaf epidermal cells by Agrobacterium-mediated transformation ( Figure 9C ). MaMYB1, MaMYB2, MaMYB3, MaNAC1 and MabHLH1 were found to activate the activity of the MaC3H2 promoter with a LUC/REN ratio higher than 2.12, 4.57, 3.81, 6.41 and 3.04-fold of empty vector (pGreen II 62-SK/pGreen II 0800-LUC), respectively, indicating that MaMYB1, MaMYB2, MaMYB3, MaNAC1 and MabHLH1 are activative effectors of MaC3H2 gene expression ( Figure 9D ).