NAC and MYB Families and Lignin Biosynthesis-Related Members Identification and Expression Analysis in Melilotus albus

Melilotus albus is an annual or biennial legume species that adapts to extreme environments via its high stress tolerance. NAC and MYB transcription factors (TFs) are involved in the regulation of lignin biosynthesis, which has not been studied in M. albus. A total of 101 MaNAC and 299 MaMYB members were identified based on M. albus genome. Chromosome distribution and synteny analysis indicated that some genes underwent tandem duplication. Ka/Ks analysis suggested that MaNACs and MaMYBs underwent strong purifying selection. Stress-, hormone- and development-related cis-elements and MYB-binding sites were identified in the promoter regions of MaNACs and MaMYBs. Five MaNACs, two MaMYBs and ten lignin biosynthesis genes were identified as presenting coexpression relationships according to weighted gene coexpression network analysis (WGCNA). Eleven and thirteen candidate MaNAC and MaMYB genes related to lignin biosynthesis were identified, respectively, and a network comprising these genes was constructed which further confirmed the MaNAC and MaMYB relationship. These candidate genes had conserved gene structures and motifs and were highly expressed in the stems and roots, and qRT-PCR further verified the expression patterns. Overall, our results provide a reference for determining the precise role of NAC and MYB genes in M. albus and may facilitate efforts to breed low-lignin-content forage cultivars in the future.


Introduction
Transcription factors (TFs) play an important role in plant growth and development. NAC (no apical meristem (NAM), Arabidopsis transcription activation factor (ATAF1/2), and cup-shaped cotyledon (CUC2)) family members are plant-specific TFs. NAC TFs have a highly conserved N-terminal region and a variable C-terminal region; the N-terminal regions share a conserved DNA-binding region, whereas the C-terminal regions have various transcriptional regulatory domains [1][2][3]. In addition, the C-terminal region of some NAC proteins has α-helical transmembrane motifs that are related to stress responses [4,5]. These features help these proteins repress or activate various downstream genes to regulate multiple molecular or cellular processes. NAC proteins are involved in the regulation of various biological functions. Overexpression of STRESS-RESPONSIVE NAC 1 significantly enhanced the drought resistance of transgenic rice in the field, with no phenotypic changes or yield penalty [6]. ATAF1, which encodes an NAC TF, is a stimulus-dependent attenuator of abscisic acid signaling involved in mediating penetration resistance [7].
The MYB gene family is one of the largest TF families. MYB proteins are characterized by a highly conserved DNA-binding domain that generally consists of up to four imperfect amino acid sequence repeats (R) of 50-53 amino acids, each encoding three α-helices [8]. The second and third helices of each repeat form a helix-turn-helix structure with regularly interspersed triplet tryptophan residues forming a hydrophobic core [9]. The functions of

Genome Distribution and Gene Synteny Analysis
Among all the NAC and MYB family members, 100 MaNACs and 287 MaMYBs were mapped onto eight chromosomes ( Figure 1). The distribution of MaNACs and MaMYBs on the chromosomes was uneven. Chromosome 5 contained the lowest density, with eight MaNACs and 24 MaMYBs (Figure 1). Chromosome 1 contained the largest number of genes, with 19 MaNACs and 42 MaMYBs. In addition, the distribution of NAC and MYB genes at chromosomal loci revealed that 30 (10.03%) and 15 (14.85%), respectively, were the result of tandem duplications (Figure 1). Eight tandem repeats of MaNAC genes were present on chromosome 2, and 11 tandem repeats of MaMYB genes were present on chromosome 1. Five tandem duplication genes of MaNAC were found on chromosome 2. Three tandem duplication genes of MaNAC were found on chromosome 8. Three and three MaMYB tandem duplication genes were found on chromosomes 1 and 7, respectively ( Figure 1). The synteny analysis showed a total of 48 NAC duplicated gene pairs and 136 MYB duplicated gene pairs in M. albus (Figure 2a, Table S3). Furthermore, some of the synteny gene pairs distributed in the same chromosome ( Figure 2a). The tandem duplication genes and synteny gene pairs distributed in the same chromosome may be caused by segmental duplication [34]. The Ks could be used to evaluate the gene duplication evolution that had occurred after divergence from their ancestors. To estimate the selection pressure on the MaNACs and MaMYBs during the M. albus evolutionary, we calculated the Ka and Ks values and Ka/Ks ratio. The median Ks and Ka/Ks values of MYB was slightly larger than that of NAC (Figure 2b,c, Table S4). The results indicate that the gene duplication events of MYB occurred earlier than that of NAC. In addition, the Ka/Ks values of all gene pairs were less than 1, which indicated that the NAC and MYB genes of M. albus experienced strong purifying selection (Figure 2c).

Phylogenetic Analysis of NAC and MYB TF Genes
Phylogenetic comparative analysis of NACs and MYBs from different plant species revealed the conservation and diversification of the NAC and MYB TF family members. To examine the evolutionary relationships among the NAC or MYB gene family members from M. albus, M. truncatula, Arabidopsis and rice, a neighbor-joining (NJ) phylogenetic tree was constructed (Figure 3a,b). Most of the subgroups contained the NACs from M. albus, M. truncatula, Arabidopsis and rice, which suggested a common ancestor before the divergence of these species (Figure 3a). Some clades included only M. albus and M. truncatula members. One clade included only M. albus NACs, indicating that these NAC genes may have been acquired in M. albus after its divergence from the last common ancestor or that those genes had been lost either in rice, M. truncatula or Arabidopsis. The gain or loss of species-specific NACs led to functional divergence. Compared with the rice (155), M. truncatula (185) and Arabidopsis (197) MYB family, the M. albus (299) MYB family is substantially larger. M. albus, rice, M. truncatula and Arabidopsis had independent clades, and clades that included only M. albus MYB members were more abundant than those that included the other three species (Figure 3b). These results indicate that the NAC family was more conserved than the MYB family and that the MYB family members evolved faster. In the MaNAC protein phylogenetic tree, all proteins were divided into six clades-group I to group VI (Figure 4a). The phylogenetic tree of the MaMYB proteins showed that these proteins were clustered into two clades: one clade included mainly R2R3-MYB proteins, and the other clade included mainly R1-MYB proteins (Figure 3c).

Phylogenetic Analysis of NAC and MYB TF Genes
Phylogenetic comparative analysis of NACs and MYBs from different plant species revealed the conservation and diversification of the NAC and MYB TF family members. To examine the evolutionary relationships among the NAC or MYB gene family members

Gene Structure and Protein Motif Analysis
To better understand the similarity and diversity of gene structure among MaNACs or MaMYBs, their intron/exon structure pattern was investigated by aligning their gene sequences and CDSs. With respect to the MaNAC TF family, 71 MaNAC genes had three exons, and only two MaNAC genes had one exon (Figure 4b). In addition, genes in the same subgroups generally contained similar intron and motif patterns (Figure 4b). For example, most of the genes in subgroups I, II and V had three exons. Thirteen motifs were identified in the MaNACs. With the exceptions of MaNAC1, MaNAC92, and MaNAC30, nearly all the MaNACs contained motif 2, and most MaNACs contained motifs 1, 2, 3, 4 and 5. For the MaMYB TF family, 11 genes had one exon, but the number of exons varied from 1 to 12 ( Figure S1). The motif structures were similar in one clade ( Figure S2). The length of MYB protein motifs varied from 11 to 51 amino acids, and the number of motifs in each MYB varied from 2 to 10 ( Figure S2).

Gene Promoter Region and Functional Analysis
The identification of cis-regulatory elements within promoter regions is beneficial for understanding the expression patterns of genes. The cis-elements in the promoter regions of MaNACs and MaMYBs were similar and comprised stress-, hormone-and development-related cis-elements (Figure 5a,b). ARE cis-elements are essential for anaerobic induction [35]. ARE cis-elements were found in the promoter regions of 83.17% and 85.28% of the MaNAC and MaMYB genes, respectively ( Figure 5a, Table S5). Lowtemperature responsiveness (LTR) cis-elements were found in the promoter regions of 38.61% and 37.46% of the MaNAC and MaMYB genes, respectively ( Figure 5a, Table S5). MBS and MBSI cis-elements are MYB-binding sites involved in drought inducibility and the regulation of flavonoid biosynthesis genes, respectively. In total, 46.53% of MaNAC and 42.14% of the MaMYB promoter regions had MBS cis-elements. MBSI cis-elements were found in 8.91% and 9.70% of the MaNAC and MaMYB gene promoter regions, respectively. Most development-related cis-elements are related to the light response, and the others were related to the endosperm, roots, seeds and the circadian rhythm ( Figure 5b). nearly all the MaNACs contained motif 2, and most MaNACs contained motifs 1, 2, 3, 4 and 5. For the MaMYB TF family, 11 genes had one exon, but the number of exons varied from 1 to 12 ( Figure S1). The motif structures were similar in one clade ( Figure S2). The length of MYB protein motifs varied from 11 to 51 amino acids, and the number of motifs in each MYB varied from 2 to 10 ( Figure S2).

Gene Promoter Region and Functional Analysis
The identification of cis-regulatory elements within promoter regions is beneficial for understanding the expression patterns of genes. The cis-elements in the promoter regions of MaNACs and MaMYBs were similar and comprised stress-, hormone-and develop ment-related cis-elements (Figure 5a,b). ARE cis-elements are essential for anaerobic in duction [35]. ARE cis-elements were found in the promoter regions of 83.17% and 85.28% of the MaNAC and MaMYB genes, respectively ( Figure 5a, Table S5). Low-temperature responsiveness (LTR) cis-elements were found in the promoter regions of 38.61% and 37.46% of the MaNAC and MaMYB genes, respectively ( Figure 5a, Table S5). MBS and MBSI cis-elements are MYB-binding sites involved in drought inducibility and the regu lation of flavonoid biosynthesis genes, respectively. In total, 46.53% of MaNAC and 42.14% of the MaMYB promoter regions had MBS cis-elements. MBSI cis-elements were found in 8.91% and 9.70% of the MaNAC and MaMYB gene promoter regions, respec tively. Most development-related cis-elements are related to the light response, and the others were related to the endosperm, roots, seeds and the circadian rhythm ( Figure 5b). Given that NAC and MYB TFs play important roles in development and stress responses, Gene Ontology (GO) classification of MaNAC and MaMYB members was performed ( Figure 5c,d, Table S6). The NAC family members were enriched mainly in biological regulation, developmental processes, and multicellular organismal processes (biological processes); cells, cell parts and organelles (cellular components); and binding (molecular functions) (Figure 5c). The MYB family members were enriched mainly in biological regu-lation, signaling and rhythmic process (biological process); cells, cell parts and organelles (cellular components); and binding, nucleic acid-binding TFs and TF activity, and protein binding (molecular functions) ( Figure 5d).

Lignin Biosynthesis Genes Regulated by NAC Genes
In the production and application of forage, lignin content significantly affects the quality of forage. NAC and MYB are known to be involved in regulating lignin biosynthesis. Using the Arabidopsis lignin-related NAC (AtVND1-7, AtNST1-3 and AtSND1-2) CDSs and protein sequences as queries, we identified six VND genes, one NST gene and four SND genes in M. albus (Figure 6a). We subsequently constructed a phylogenetic tree of M. albus and Arabidopsis NAC proteins according to the maximum likelihood method of MEGA 7.0 ( Figure 6a). The 24 NAC members of M. albus and Arabidopsis were subdivided into three subgroups: the VND subgroup, NST subgroup and SND subgroup ( Figure 6a). The similarity of amino acids was higher within the same subgroup than across subgroups ( Figure S3). With the exceptions of AtVND3 and AtVND6, all these lignin biosynthesisrelated genes had similar structures, with three exons (Figure 6b). Proteins in the VND and NST clades had the seven motifs, and proteins in the SND clade had four motifs ( Figure 6c). All proteins had NAM domains that spanned motifs 1, 2, 3 and 5 in the VND and NST clades and motifs 1, 5, 6 and 7 in the SND clade ( Figure 6c). The NAM domain with 153 amino acids was visualized by Gene Structures Display Server (GSDS) (Figure 6d). All amino acid sequences of Arabidopsis and M. albus were aligned by DNAMAN software; the results are shown in Figure S3.

NAC and MYB Gene Expression and Coexpression Analysis
To elucidate the possible roles of NACs and MYBs in the lignin biosynthesis, growth and development of M. albus, we analyzed RNA-seq data from different tissues. A total of 63 MaNAC genes and 185 MaMYBs were expressed in the roots, stems, flowers and seeds

NAC and MYB Gene Expression and Coexpression Analysis
To elucidate the possible roles of NACs and MYBs in the lignin biosynthesis, growth and development of M. albus, we analyzed RNA-seq data from different tissues. A total of 63 MaNAC genes and 185 MaMYBs were expressed in the roots, stems, flowers and seeds (Figure 8a,b, Tables S1 and S2). Some of the MaNACs and MaMYBs were not expressed in any tissue, which indicated that they were functionally redundant. Thirty-two MaNACs and 104 MaMYBs were expressed in all tissues. A total of 44 MaNAC genes and 108 MaMYBs presented tissue-specific expression (Figure 8c,d), among which 6 MaNACs and 14 MaMYBs were only expressed in the roots, one MaNAC and one MaMYB were only expressed in the leaf, one MaMYB was only expressed in the stem, three MaNACs and one MaMYB were only expressed in the seed, and three MaNACs and 20 MaMYBs were only expressed in the flower. A total of 45 MaNAC genes and 128 MaMYB genes were expressed in five leaf samples ( Figure S5). The expression level analysis suggested that tandem repeat genes exhibited similar expression patterns. For example, the expression levels of five MaNAC tandem repeat genes were notably low, and the expression patterns of three MaNAC tandem repeat genes were similar (Figure 8, Figure S5).
To determine the relationship of MaNAC and MaMYB TFs further, we performed a weighted gene coexpression network analysis (WGCNA) based on RNA-seq data of the roots, stems, flowers, and seeds as well as five leaf samples, with three replications. Thirteen MaNACs, 34 MaMYBs, and 53 lignin biosynthesis-related genes were identified (Table S7)

Network of Lignin Biosynthesis-Related MaNACs and MaMYBs
The first layer of master switch genes (SND1/NST3, NST1, VND6, and VND7) and the second layer of master switch genes (MYB46 and MYB83) regulate downstream genes to control lignin biosynthesis during secondary wall formation [19]. MYB46 and MYB83 are direct targets of NST3/SND1, VND6 and VND7 [20,[40][41][42]. MYB103 is a direct target of NST1, which is involved in syringyl lignin biosynthesis [43] and may be a direct target of VND6 and VND7 [41,42]. AtMYB58, AtMYB63, and AtMYB103 are direct targets of AtMYB46/83 [22]. AtMYB20, AtMYB42, AtMYB43, and AtMYB85 activate AtMYB4, AtMYB7, and AtMYB32 and directly activate lignin biosynthesis genes during secondary wall formation in Arabidopsis [22]. AtMYB58, AtMYB63, and AtMYB103 regulate lignin biosynthesis genes [19,21,44], and At-MYB4, AtMYB32, and AtMYB7 negatively regulate lignin biosynthesis genes [19]. On the basis of previous studies, we assumed a MaNAC and MaMYB coregulation network (Figure 9b). The expression pattern of MaMYB234 was opposite that of MaNAC43 and MaNAC8 (VND) and MaMYB29. The expression pattern of MaMYB234 was similar to that of MaNAC78 and MaNAC32, MaMYB118, MaMYB235 and MaMYB249 and MaMYB266 (Figure 9b). In the different examined tissues, most of the genes (16/24) were highly expressed in the stems. However, MaMYB118, MaMYB266, MaNAC8 and MaNAC65 (VND) were highly expressed in the roots. To further validate the reliability of gene expression data obtained by RNA-Seq analysis, four lignin biosynthesis-related NAC and MYB genes each carried out qRT-PCR analysis in five tissues ( Figure 10). NAC and MYB genes were highly expressed in stems and roots, which were consistent with the RNA-seq analysis, with the exception of three genes higher expressed in leaf than in stem or root. MaNAC tandem repeat genes were notably low, and the expression patterns of three MaNAC tandem repeat genes were similar (Figure 8, Figure S5).

Network of Lignin Biosynthesis-Related MaNACs and MaMYBs
The first layer of master switch genes (SND1/NST3, NST1, VND6, and VND7) and the second layer of master switch genes (MYB46 and MYB83) regulate downstream genes to control lignin biosynthesis during secondary wall formation [19]. MYB46 and MYB83 are direct targets of NST3/SND1, VND6 and VND7 [20,[40][41][42]. MYB103 is a direct target of NST1, which is involved in syringyl lignin biosynthesis [43] and may be a direct target of VND6 and VND7 [41,42]. AtMYB58, AtMYB63, and AtMYB103 are direct targets of AtMYB46/83 [22]. AtMYB20, AtMYB42, AtMYB43, and AtMYB85 activate AtMYB4, AtMYB7, and AtMYB32 and directly activate lignin biosynthesis genes during secondary wall formation in Arabidopsis [22]. AtMYB58, AtMYB63, and AtMYB103 regulate lignin biosynthesis genes [19,21,44], and AtMYB4, AtMYB32, and AtMYB7 negatively regulate lignin biosynthesis genes [19]. On the basis of previous studies, we assumed a MaNAC and MaMYB coregulation network (Figure 9b). The expression pattern of MaMYB234 was opposite that of MaNAC43 and MaNAC8 (VND) and MaMYB29. The expression pattern of MaMYB234 was similar to that of MaNAC78 and MaNAC32, MaMYB118, MaMYB235 and MaMYB249 and MaMYB266 (Figure 9b). In the different examined tissues, most of the genes (16/24) were highly expressed in the stems. However, MaMYB118, MaMYB266, MaNAC8 and MaNAC65 (VND) were highly expressed in the roots. To further validate the reliability of gene expression data obtained by RNA-Seq analysis, four lignin biosynthesisrelated NAC and MYB genes each carried out qRT-PCR analysis in five tissues ( Figure 10). NAC and MYB genes were highly expressed in stems and roots, which were consistent with the RNA-seq analysis, with the exception of three genes higher expressed in leaf than in stem or root.

Discussion
NACs and MYBs compose two large TF families. As an increasing number of genomes of legume species are sequenced, additional NAC and MYB TF family members are being identified. For instance, there are, respectively, 269 and 438 NACs and MYBs in Glycine max, 123 and 185 in M. truncatula, 116 and 126 in Lotus japonicus, and 96 and 166 in Cicer arietinum (data from PlantTFDB 5.0, http://planttfdb.cbi.pku.edu.cn/index.php?sp=Mtr, accessed on 28 February 2020). In the present study, based on whole-genome data of the M. albus, 101 MaNAC and 299 MaMYB genes were identified. Because G. max underwent independent whole-genome duplication, there are substantially more of these family members in this species than other legume species. The number of NACs is similar between M. albus and both L. japonicus and C. arietinum, but there are substantially more MYBs in M. albus than in M. truncatula, L. japonicus and C. arietinum. M. albus can adapt to drought, cold and high-salinity environments [45], and the expansive MaMYB family may contribute to the increased stress tolerance of M. albus. The phylogenetic tree M. albus, M. truncatula, rice and Arabidopsis showed that most clades contained NAC or MYB members of all three species, but some clades contained members of only one species. In addition, the intron/exon and motif structures of NAC and MYB genes were similar in one clade. Genes in the same clade may share the same ancestor before species divergence, and genes in independent clades with genes from only one species may have formed after species divergence.
Previous studies have shown that NAC TFs regulate cell division [46] and flowering in response to salt stress [47] and are involved in the regulation of the drought response [1] and in AtNAC1-mediated auxin signaling to promote lateral root development [48]. In addition, the various functions of MYB TFs in plants include lignin biosynthesis [49,50], hormone signal transduction [51,52], stress tolerance [53,54], and organ development [55,56]. The promoter regions of MaNACs and MaMYBs contain stress-, development-and hormoneresponsive cis-elements. Some of the gene promoter regions contain ARE cis-elements (>80%) and LTR cis-elements (>37%), which are involved in stress responses [35]. CAT-box cis-elements are related to meristem expression, GCN4_motif cis-elements are involved in endosperm expression, CGTCA-motif cis-elements are involved in methyl jasmonate (MeJA) responsiveness, and TCA elements are involved in salicylic acid responsiveness [35]. The various cis-elements within the promoter regions of MaNACs and MaMYBs are strongly associated with gene function. Furthermore, GO annotation of the MaNACs and MaMYBs indicated that they are associated with stimulus response and development, and tissuespecific expression of the MaNAC and MaMYB genes indicated that these genes may regulate organ development. Furthermore, 23 MaMYB genes were expressed only in the flowers, suggesting that they may regulate M. albus floral development. Therefore, MaNAC and MaMYB members may be involved in stress, development and hormone responses.
Some NAC members target MYB members that are involved in lignin biosynthesis [19]. Analysis of the cis-elements within the promoter regions revealed that the promoters of some MaNACs and MaMYBs contained MBS and MBSI cis-elements. Moreover, WGCNA indicated that MaNAC and MaMYB members presented coexpression relationships. Our results are consistent with those of previous studies, which showed that NAC members regulate MYB members. On the basis of the relationship between MYBs and NACs, and MYBs in Arabidopsis, a network comprising MaNACs and MaMYBs with respect to lignin biosynthesis in M. albus was constructed which further confirmed the MaNAC and MaMYB relationship. Stems and roots are highly lignified, and our sequencing results are consistent with this fact. Comparing gene expression levels in different tissues, we identified genes that are highly expressed in the stems and roots, as well as flower-specific genes. Furthermore, the results of qRT-PCR also showed higher relative expression level in stems and roots. VND and NST subfamily members may function redundantly [19], compared with seven AtVNDs and three AtNSTs in Arabidopsis, only six VNDs and one NST were identified in M. albus. Two SND2s, two SND3s, two MYB83s, two MYB20s, and two MYB103s were identified that had expanded in M. albus compared with Arabidopsis.
Except for MaMYB284, all these genes were expressed in different tissues. These expanded genes may increase M. albus lignin biosynthesis, as mature M. albus plants have a greater degree of lignification than mature Arabidopsis plants. Sequences of the Arabidopsis lignin biosynthesis-related NAC and MYB genes were retrieved from previous studies and used as queries to identify the corresponding M. albus genes. Local BLAST searches were performed: an e-value cut-off of 1e −100 and 1e −70 was applied to NAC and MYB homologous sequence recognition, respectively (Table S9).

Chromosomal Location, Motif Analysis, and Gene Structure of the NAC and MYB Genes in M. albus
Each MaNAC and MaMYB gene was mapped to the M. albus genome to determine its position on the chromosome via TBtools with default parameters [57]. The protein sequences of the NAC and MYB family members were submitted to MEME (http:// meme-suite.org/tools/meme, accessed on 25 March 2020) (with the default parameters) to identify conserved motifs. To examine the conserved domains of M. albus in detail, multiple amino acid sequence alignment was performed via DNAMAN 6.0 software. The multiple alignment sequences were submitted to WebLogo (http://weblogo.berkeley.edu/ logo.cgi, accessed on 3 April 2020) to create sequence logos of NAM and R2 and R3 Myb repeats. Furthermore, the CDSs of NAC and MYB were compared with their corresponding genomic sequences to determine all gene structures, which were visualized using the Gene Structure Display Server (http://gsds.cbi.pku.edu.cn/, accessed on 10 April 2020) with default parameters.

Phylogenetic Analysis
The NAC and MYB protein sequences of Arabidopsis, O. sativa and M. albus were aligned via the Muscle method, which is the MEGA 7.0 built-in method. Phylogenetic trees of NAC and MYB family members were then constructed by MEGA 7.0 [58] with the neighbor-joining (NJ) method (bootstraps with 1000 replicates). Lignin biosynthesisrelated NAC and MYB phylogenetic trees were subsequently constructed according to the maximum likelihood method.

Analysis of the Promoter Regions of MaNAC and MaMYB Genes
The 2000-bp upstream of the transcriptional start site of MaNAC and MaMYB genes was extracted from the M. albus genome to identify the cis-elements within the promoter regions. The plant Cis-Acting Regulatory Element (http://bioinformatics.psb.ugent.be/ webtools/plantcare/html/, accessed on 9 May 2020) database was used to identify ciselements within the promoter regions of each gene (with the default parameters), and the function of each cis-element was annotated.

Synteny Analysis of MYB and NAC Genes
Genes within families were aligned with one another by the use of BLAST (e-value < 1 × 10 −10 ), and the best matched gene pairs were defined as homologous pairs. The synonymous substitution rates (Ks) and nonsynonymous substitution rates (Ka) of the gene pairs were calculated by TBtools [57].
The syntenic relationships of M. albus genes were determined via MCScanX [59], with the default parameters. The MYB and NAC gene pairs were then selected from the MCScanX results. Circos [60] was subsequently used to visualize the syntenic relationships.

Expression Analysis and Gene Annotation
Seeds were grown in turfy soil in a glasshouse, which was maintained at 28/24 • C (day/night) and provided a 16 h light/8 h dark photoperiod. During the inflorescence stage, tissue samples of fresh stems, flowers and roots were collected, immediately frozen in liquid nitrogen and then stored at −80 • C. Seed samples were mixed together, imbibed and then allowed to germinate for 2 days. Each sample was a composite of independent biological replicates.
Total RNA from each sample was isolated by the use of TRIzol reagent (Invitrogen) according to the manufacturer's instructions. Messenger RNAs (mRNAs) were separated from the total RNA by oligo (dT) primers and were cleaved into short fragments at random. Sequencing libraries were generated via rRNA-depleted RNA with a NEBNext ® Ultra™ Directional RNA Library Prep Kit for Illumina ® (NEBNext, Ipswich, MA, USA) according to the manufacturer's recommendations. The qualified cDNA libraries were constructed by PCR enrichment and sequenced on a HiSeq 2500 platform, with a read length of 125 bp. Three independent biological replicates were sequenced per sample. The RNA sequencing (RNA-seq) data of leaf samples from 5 different genotypes from a previous study were also included [32].
Clean reads were obtained by removing reads containing adapters, reads containing poly-N sequences and reads of low quality from the raw reads. The clean reads were mapped to the M. albus genome by HISAT2. StringTie 1.3.1 [61] was used to calculate fragments per kilobase of transcript per million mapped reads (FPKMs) of mRNAs in each sample. The FPKMs of the genes were computed by summing the FPKMs of the transcripts in each gene group. Genes with FPKM < 0.5 were defined as non-expressed genes.
All gene expression levels were used to calculate the weight of gene pairs via WGCNA. The relationships of the MYB and NAC gene pairs were visualized by Cytoscape 3.5.1, and gene functions were annotated via the Gene Ontology (GO) database.
PrimeScript ® RT reagent Kit (TaKaRa, Kyoto, Japan) was used to reverse transcribe RNA into first-stand cDNA. The product was used as a template for qRT-PCR with specific primers (Table S10) which were designed by PerlPrimer v1.1.21. MaTublin was used as the reference gene. qRT-PCR was performed using SYBR Premix Ex TaqTM (TaKaRa). Each experiment was performed with three replicates. The relative transcription levels were calculated by the comparative delta-delta Ct method.

Conclusions
In summary, a total of 101 MaNAC and 299 MaMYB family members were identified in M. albus. Phylogenetic and evolutionary relationships, gene structure analysis and motif analysis indicated the conservation of MaNAC and MaMYB families and that some members diverged from their ancestors. MaNAC and MaMYB members experienced purifying selection and some of the members were functionally redundant. Promoter region analysis and GO functional annotation of MaNAC and MaMYB TFs suggested that they may be involved in stress, development and hormone responses. WGCNA revealed that MaNAC, MaMYB and lignin biosynthesis-related genes presented coexpression relationships, and a regulatory network comprising some of these genes in M. albus further confirmed their relationship. Gene expression levels in different tissues were compared, some of the MaNAC and MaMYB members presented tissue-specific expression, lignin biosynthesis-related MaNACs and MaMYBs were highly expressed in stem and root, and the relative expression of several members was also high in stem and root, which further confirmed that these candidate genes may regulated lignin biosynthesis. These results provide a useful framework for future research to understand the evolution and function of the MaNAC and MaMYB gene families and to understand the potential physiological role of the MaNAC and MaMYB genes during lignin biosynthesis. Since no MaNAC and MaMYB genes have been functionally characterized to date, our results may facilitate future efforts to breed forage cultivars with low amounts of lignin.

Conflicts of Interest:
The authors declare no conflict of interest.