Comparative transcriptome analyses reveal that the MsNST1 gene affects lignin synthesis in alfalfa (Medicago sativa L.)

Qing Zhou, Pei Mo, Dong Luo, Xutin Chi, Ho Deng, Qingen Fng, Longf Fng, Zhiio Nn,Jingqi Wen, Zhipeng Liu,

a State Key Laboratory of Grassland Agro-ecosystems, Key Laboratory of Grassland Livestock Industry Innovation, Ministry of Agriculture and Rural Affairs, Engineering Research Center of Grassland Industry, Ministry of Education, College of Pastoral Agriculture Science and Technology, Lanzhou University, Lanzhou 730020, Gansu, China

b Xinjiang Key Laboratory of Desert Plant Roots Ecology and Vegetation Restoration, Xinjiang Institute of Ecology and Geography, Chinese Academy of Sciences, Urumqi 830011,Xinjiang, China

c College of Grassland Science/Key Laboratory of Grassland Ecosystem of the Ministry of Education, Gansu Agricultural University, Lanzhou 730020, Gansu, China

d Noble Research Institute, LLC, Ardmore, OK 73401, USA

Keywords:Medicago sativa L.Lignin synthesis Transcriptome WGCNA MsNST1

A B S T R A C T As the second most abundant natural polymer,accounting for approximately 30%of the organic carbon in the biosphere,lignin plays an essential role in plant development.However,a high lignin content affects the nutritional quality of alfalfa(Medicago sativa L.),the most widely cultivated perennial legume forage crop.Histological analysis indicated that G-lignin and S-lignin were present in the stem,leaf,and petiole of alfalfa, and the deposition of lignin increased gradually in descending internodes. Neutral detergent fiber (NDF), acid detergent fiber (ADF), and acid detergent lignin (ADL) contents continually increased from the top to the bottom of the stem,and ADL content showed a similar trend in leaves.Alfalfa leaves and stems from five different nodes(1,2,4,6,and 8)were used as materials to investigate molecular regulatory mechanisms in lignin synthesis by RNA sequencing. Respectively 8074 and 7752 differentially expressed genes (DEGs) were identified in leaves and stems, and 1694 DEGs were common to the two tissues. ‘‘Phenylpropanoid biosynthesis” was the most enriched pathway in both leaves and stems, and 134 key regulatory genes in lignin synthesis were identified by a weighted gene co-expression network analysis.The NAC family transcription factor MsNST1 gene was highly expressed in old leaf and stem tissues.The deposition pattern of G-and S-lignin differed among M.truncatula wild-type,nst1 mutants,and overexpression lines,and the transcription levels of lignin synthesis genes such as HCT,F5H,and COMT in these three materials also differed. These results suggest that MsNST1 affects lignin synthesis in alfalfa.These findings provide a genetic basis and abundant gene resources for further study of the molecular mechanisms of lignin synthesis, laying a foundation for low-lignin alfalfa breeding research.

1. Introduction

Lignin is a natural aromatic polymer in plant cell walls that contributes cell wall integrity, stem strength, water transport, and pathogen resistance[1],and strongly affects the utilization of sustainable lignocellulosic biomass [2]. As the second most plentiful natural polymer on earth after cellulose,lignin accounts for around 30% of organic carbon [3]. Generally, lignin is an underutilized feedstock in the growing pulping and biofuel industries and is used as a fuel source instead. However, in recent decades, advances in lignin applications include its use in the synthesis of polymers,dyes,adhesives,and fertilizers[4,5].Lignin can be used for wound dressings, pharmaceuticals, and drug delivery materials in the medical sector and also for electrochemical energy materials and 3D printing of lignin-plastic composite materials [6-8].

Lignin is a phenolic polymer thought to be formed by oxidative polymerization of the five primary hydroxycinamyl alcohols(monolignols)[9],including p-coumaryl,caffeyl,coniferyl,sinapyl,and 5-hydroxylcinnamyl alcohols, which give rise to phydroxylphenyl (H), catechyl (C), guaiacyl (G), syringyl (S), and 5-hydroxyguaiacyl (5H/5-OH-G) lignin. Lignin composition varies among plant species and tissue types.Lignins are composed mainly of G and S units in vascular tissues, whereas seed lignins vary in their monolignol composition. In gymnosperm, lignins contain mostly G units and very low levels of H units. Lignins from dicots are composed of both G and S units with only traces of H units,whereas in monocot grasses they contain G and S units along with higher levels of H units (but still ＜5%) [10]. Lignin structure and the monomer ratio influence the use of wood fiber in biomass energy production, pulping, and papermaking, and the properties of lignin are influenced by its monomer ratio [11]. A higher S/G value contributes to delignification and accelerates the decomposition of lignocellulose in these industries [12].

Following the biosynthesis of monolignols,lignin monomers are translocated to the secondary plant cell wall, where they are oxidized for polymerization [2]. A series of enzymes are involved in the lignin biosynthetic pathway,including phenylalanine ammonia lyase (PAL), cinnamate 4-hydroxylase (C4H), 4-coumarate CoA ligase (4CL), shikimate/quinate hydroxy cinnamoyl transferase(HCT), p-coumarate 3-hydroxylase (C3H), caffeoyl CoA 3-Omethyltransferase (CCoAOMT), cinnamoyl CoA reductase (CCR),caffeic acid O-methyltransferase (COMT), cinnamyl alcohol dehydrogenase (CAD), and ferulate 5-hydroxylase (F5H) [13]. These related genes have different functions in the lignin synthesis pathway.COMTis responsible for O-methylation at the C5 position of the phenolic ring, while the function ofCCoAOMTis in Omethylation of the C3 position, and the methylation affects the efficiency of lignin synthesis and lignin composition [14]. The S unit, H unit, and S/G ratio inC3H-down-regulated transgenic poplar are increased, and its G unit is decreased; but the S unit and S/G ratio inHCT-down-regulated transgenic poplar are decreased while its G unit is increased [12]. Transcription factors are also involved in the lignin biosynthesis pathway;NAC SECONDARY WALL THICKENING PROMOTING FACTOR(NST), a NAC-family transcription factor, affects the entire biosynthetic pathways for cellulose, xylan, and lignin inArabidopsisandMedicago truncatula[15-17]. One effective method to improve plants for biofuel production is to alter the expression of these genes in the lignin biosynthesis pathway.

High-throughput sequencing has dramatically improved the efficiency of transcriptome data collection. It offers clear advantages over microarray methods and has been used in gene discovery and regulatory network studies [18], such as investigating the response of plants to abiotic stress[19,20],molecular mechanisms of seed mucilage formation[21,22],and flavor volatiles biosynthesis [23]. The method has also been used in studying the lignin biosynthetic pathway, including identification of key genes. Transcriptome profiles ofPaphiopedilum armeniacumseed at several development stages were compared to identify the molecular clues for non-methylated lignin synthesis, and many genes associated with phenylpropanoid biosynthesis and phenylalanine metabolism during seed maturation were differentially expressed[9].The IAAand ABA-responsive transcription factorCgMYB58was identified by transcriptomic analyses using the juice sacs of threeCitrus grandiscultivars at 55, 85, 115, 145, and 205 days post-anthesis, and overexpression ofCgMYB58promoted lignin accumulation in the pummelo fruit mesocarp and up-regulated 19 lignin biosynthetic genes[13].Establishing the lignin synthesis pathway and identifying the mechanism of critical genes will contribute to improving the efficiency of lignin utilization.

Alfalfa(Medicago sativaL.)is the most widely cultivated perennial legume forage crop, with more than 40 million hectares planted worldwide [24]. In addition to being high-quality forage,alfalfa offers potential as a feedstock for the production of ethanol and other bioproducts[25].However,alfalfa contains a high lignin content compared to grass, hindering the degradation of other components [26]. The cultivation of low-lignin cultivars has become the main direction of alfalfa breeding research, and reduced-lignin alfalfa cultivars are now available and have the potential to increase the digestibility of alfalfa forage [27]. The Monsanto Company (St. Louis, MO, USA) developed a low-lignin alfalfa cultivar KK179 by down-regulation ofCCoAOMT, which is one of the lignin biosynthesis genes, and this cultivar has been approved in many countries [28]. Some genes have been shown to regulate lignin synthesis in alfalfa, includingMsASMT1[29],COMT[30,31],4CL[32],TT8, andHB12[28]. The identification of these lignin synthesis genes will promote the development of low-lignin alfalfa cultivars, but the lignin biosynthesis pathway in alfalfa is still unknown. The objective of the present study was to characterize the lignin synthesis pathway of alfalfa by transcriptome analysis of tissue samples collected from multiple plant organs.

2. Materials and methods

2.1. Plant materials and sample collection

The alfalfa cultivar Longzhong was grown at the Yuzhong Experimental Station of Lanzhou University (N35°57′, E104°09′,1720 m above sea level),Lanzhou,China.The mean annual precipitation at the site is 350 mm and the mean annual temperature is 6.7 °C. Ten sample types were collected in the early bud stage,comprising leaves(including the entire petiole)and stems(excluding the node part)from five nodes(1,2,4,6,and 8,counting from the top of the stems).These samples were named L1,L2,L4,L6,L8,S1, S2, S4, S6, and S8 (Figs. 1A, S1). All samples were harvested from 60 individual plants,and three biological replicates per sample type were prepared. These materials were dried to constant weight at 60°C before the determination of neutral detergent fiber(NDF), acid detergent fiber (ADF), and acid detergent lignin (ADL).Fresh materials were selected to characterize the distribution of lignin in the sample types.For RNA isolation,six plants(three biological replicates containing two plants per replicate) were selected to collect the 10 sample types.Samples were immediately frozen in liquid nitrogen and stored at -80 °C.

2.2. Determination of NDF, ADF, and ADL

All samples were ground into powder with an SM2000 mill(Retsch, Haan, Germany) and passed through a 0.5 mm screen.The contents of NDF, ADF, and ADL were determined as described previously [33]. All data were subjected to one-way analysis of variance (ANOVA) for a completely randomized design using SPSS 19.0(IBM,Armonk,New York,NY,USA).Differences among means were evaluated by Duncan’s multiple range test atP＜0.05.

2.3. Histological analysis

To characterize changes in lignin in the plant samples, the distribution of lignin was investigated as in previous studies [34,35].Fresh samples were cut into 20-25 μm thin slices along a cross section. For guaiacyl (G) lignin, the samples were stained for 3 min with 2% m-phloroglucinol ethanol and 1 min with 3 mol L-1HCl,and for syringyl (S) lignin, samples were stained for 2 min with 0.5%KMnO4,1 min with 3 mol L-1HCl,and 1 min with 10%ammonia solution.After staining,the distribution of lignin was observed under a Nikon Microphot FXA microscope (Nikon Corporation,Tokyo, Japan).

2.4.RNA extraction and library construction for transcriptome analysis

Fig.1. Alfalfa accumulated high lignin content in their old leaves and stems.A total of 10 sample types were collected in the early bud stage of alfalfa,including the leaves and stems from five different nodes: 1, 2, 4, 6, and 8 (A); and NDF, ADF, and ADL contents were determined (B). Different letters above bars indicate significant differences at P ＜0.05 by Duncan’s multiple range test. Scale bars, 1 cm.

Total RNA of each sample was isolated using the Trizol method(Sangon Biotech,Shanghai,China)according to the manufacturer’s instructions.RNA concentration was determined using a NanoDrop ND8000 spectrophotometer (Thermo Scientific, Waltham, MA,USA), and RNA integrity was checked with an RNA Nano 6000 Assay Kit with an Agilent Bioanalyzer 2100 system (Agilent Technologies,Santa Clara,CA,USA).Thirty cDNA libraries(three biological replicates for each sample type) were constructed using an mRNA-Seq Sample Preparation Kit, according to the manufacturer’s instructions (MGIEasy mRNA library preparation kit, MGI,Shenzhen, Guangdong, China). Poly (A)-containing mRNA was enriched by oligo magnetic adsorption, and the enriched mRNA was fragmented. Double-stranded cDNA (dscDNA) was then synthesized by reverse transcription with an N6 random primer.Purification, adapter ligation, and PCR amplification were performed according to the RNA-Seq protocol. These libraries were sequenced on a BGISEQ-500 platform at the Beijing Genomics Institute (https://www.genomics.org.cn, BGI, Shenzhen, Guangdong, China).

2.5. Sequence filtering, de novo assembly, and functional annotation

Sequencing read filtering was performed before downstream analyses because raw reads contain low-quality, adapter-polluted reads and also have a high content of unknown-base (N) reads.Trinity (Broad Institute, Cambridge, MA, USA) was then used to perform de novo assembly with clean reads, and Tgicl was used to cluster transcripts into unigenes [36]. The unigenes were divided into two types: clusters, for which the prefix is ‘‘CL” with the cluster ID following it (one cluster contains several unigenes,and the similarity between them is more than 70%)and singletons,for which the prefix is ‘‘Unigene”. The final unigenes were annotated using five public databases, including the NCBI Nucleotide Sequence Database (Nt), NCBI Non-redundant Protein Sequences(Nr), Gene Ontology (GO), euKaryotic Orthologous Groups (KOG),and Kyoto Encyclopedia of Genes and Genomes(KEGG)databases.Unigenes were aligned to the Nt, Nr, KOG, and KEGG databases using Blastn and Blastx (NCBI, Rockville Pike, Bethesda, MD, USA)with default parameters, and Blast2GO (BioBam, Valencia, Spain)was used to perform GO annotation.

2.6. Quantitative real-time PCR

To verify the reliability of RNA-Seq, eight unigenes were selected for a quantitative real-time PCR (qRT-PCR) test. qRT-PCR was performed using 2×SG Fast qPCR Master Mix(Sangon Biotech,Shanghai, China) on a CFX96 Touch Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA) and the original data was analyzed with CFX Manager software (Bio-Rad). qRT-PCR was programmed as follows: 95 °C for 3 min and 39 cycles of 95 °C for 10 s and 55 °C for 30 s. The specific primers were designed using Primer6 software (Premier Biosoft International, Palo Alto,CA, USA) as shown in Table S1. Three technical replicates were completed for each sample, and relative expression levels were normalized to the expression of theMedicago actingene and calculated using the 2-ΔΔCTmethod [37].

2.7. Differential expression of unigenes

Transcript expression levels were identified with the RSEM software package (https://deweylab.github.io/RSEM/) and expressed as fragments per kilobase of transcript per million mapped transcript (FPKM) method for each sample. Based on the mean FPKM values in each sample, differential expression among the leaf and stem samples was assessed with NOISeq(https://www.bioconductor.org/). Differentially expressed genes (DEGs) were assigned by the following criteria: absolute value of fold change ≥4 and adjustedP-value ≤0.05.

2.8. DEGs analysis

Cluster analysis and expression pattern assessment of DEGs were performed with the Pheatmap R package (https://cran.r-project.org/web/packages/). For annotation of the biological functions of the DEGs,GO and KEGG pathway enrichment analyses were performed with agriGO 2.0 (https://systemsbiology.cau.edu.cn/agri-GOv2/) and KOBAS 3.0 (https://kobas.cbi.pku.edu.cn/),respectively.To predict transcription factors(TFs)in DEGs,the open reading frame (ORF) of each unigene was examined with Getorf[38] and aligned to TF domains using Hmmsearch [39]. TFs were identified according to the rules described in the PlnTFDB database(https://plntfdb.bio.uni-potsdam.de/v3.0/). Weighted gene coexpression network analysis(WGCNA)was performed using the R package WGCNA 4.0.3 to further identify key genes in alfalfa lignin synthesis [40]. An applicable soft-thresholding power based on a scale-free topology criterion was employed to transform the correlation matrix to a signed weighted adjacency matrix by calculating the Pearson correlations between FPKM of all genes. The resulting adjacency matrix was used to calculate a topological overlap matrix(TOM).All genes were then hierarchically clustered based on similarity using a dynamic tree cut to form a module with a threshold of minimum module size of 20 genes[22].

2.9. Phenotypic analysis of transgenic plants

The full-length coding DNA sequence(CDS)ofMsNST1was isolated from the stem of the alfalfa cultivar Longzhong and cloned into a pBI121 overexpression vector. The vector was then introduced into theAgrobacterium tumefaciensstrain EHA105 for transformation of calli ofM. truncatula. Callus transformation and growth conditions were established using the stable transformation protocol described in the previous report[41].TwoM.truncatula nst1mutants (NF12301,nst1-1; NF15718,nst1-2) were obtained from the Noble Research Institute. Thus, a wild type(WT), two T-DNA insertion lines (nst1-1andnst1-2), and two OEMsNST1lines (OE7 and OE14) were used. AllM. truncatulaplants were of the R108 genotype. Plants were grown in peat soil at 22 °C, with a 16-h day/8-h night photoperiod, 70%-80% relative humidity, and 150 μmol m-2s-1light intensity.

The expression ofMsNST1in leaves and stems from five nodes was measured by qRT-PCR. To quantify the relative expression ofNST1and genes involved in lignin synthesis (4CL,C3H,CCoAOMT,CCR,COMT,F5H, andHCT), qRT-PCR was performed with RNA isolated from stems (internodes 1-4, counting from the top of stems) of WT,nst1-1,nst1-2, OE7, and OE14. Gene expression was calculated by the 2-ΔΔCTmethod andMtactinwas used as a reference gene.The specific primers are listed in Table S1.The fourth internode of these five 12-week-oldM. truncatulamaterials was used to characterize changes in lignin by histological analysis.Stem length, stem diameter, and number of branches of 15-week-old material were measured. The lignin content of the stem material (all internodes) was determined according to the manufacturer’s instructions (MZS-1-G, Komin, Suzhou, Jiangsu, China),and lignin composition was determined by nitrobenzene oxidation[42]. Shoot tissue of 15-week-oldM. truncatulamaterials was collected, and all the samples were ground to powder and sieved as described above.Then the NDF,ADF,crude protein,crude fat,crude ash, calcium ion, potassium ion, phosphorus ion contents,and relative feed value (RFV) of all samples were measured with a near infrared spectrometer (Unity SpectraStar 1400XTR, Unity, USA),according to the constructed prediction model of alfalfa.

Fig.3. Scatter plot of enriched KEGG pathways for DEGs in leaves(A)and stems(B).Only the 20 most strongly represented pathways are displayed.The enrichment factor is the ratio of the total of annotated genes to the DEG number in a given pathway. The colors of dots represent the range of the - log10 (Q-value).

3. Results

3.1. NDF, ADF, and ADL contents and histological analysis

NDF,ADF,and ADL contents were measured to characterize the lignification level and nutritional value of leaves and stems from multiple alfalfa nodes. As shown in Fig. 1B, NDF and ADF contents increased continually in S1-S8, from 37% to 63% and 16% to 44%,respectively, but there was no significant difference in L1-L8.ADL content also increased continually in S1-S8 and increased(P＜0.05) in L1-L2 and L6-L8. NDF, ADF, and ADL contents were higher in stems, except for S1 and S2, than in leaves.

Fig. 2. Histochemical staining of stems from five nodes in alfalfa. G-lignin was stained with m-phloroglucinol ethanol and HCl (A-J),and S-lignin was stained with KMnO4,HCl, and ammonia solution (K-T). Figures F-J and P-T (Scale bars, 200 μm) are enlarged images of A-E and K-O (Scale bars, 500 μm), respectively. Co, collenchyma; Pf,primary phloem fiber;Pi,pith;Pp,primary phloem;Px,primary xylem;Sx,secondary xylem;Xv,primary xylem vessel.S1 to S8 indicate stems from five different nodes(1,2,4, 6, and 8, counting from the top of the stems).

To further investigate the dynamics of lignin changes, the distribution of lignin in the sample types was observed by histological analysis. Both G-lignin and S-lignin were present in the stem,leaf,and petiole of alfalfa,and lignin deposition increased gradually with lower internodes (Figs. 2, S2). S-lignin occurred later than G-lignin in the three alfalfa tissues, and there were differences in the distribution of G-lignin and S-lignin among these three tissues.

3.2. De novo transcriptome assembly and functional annotation of unigenes

The 30 cDNA libraries yielded 3579.59 Mb raw reads(Table S2).Read cleaning yielded 330.62 Gb clean bases,and each library was longer than 10 Gb.Finally,199,583 unigenes were assembled,with mean, N50, N70 and N90 lengths of 1332, 2038, 1435, and 663 nt,respectively(Table S2).Additional detail of the length distribution of all unigenes is provided in Fig. S3.

Of these 199,583 unigenes,129,197(64.73%),163,325(81.83%),79,106 (39.64%), 93,566 (46.88%), and 94,737 (47.47%) unigenes were successfully annotated in the Nr, Nt, GO, KOG, and KEGG databases, respectively (Fig. S4). A total of 56,975 (28.55%) unigenes were annotated in all five databases,and 173,387(86.87%)unigenes were annotated in at least one database.

3.3. qRT-PCR verification

To verify the reliability and reproducibility of the RNA-Seq analysis,eight unigenes were randomly selected for qRT-PCR validation(Table S1). The expression of these selected unigenes in our transcriptome data was generally consistent with that measured by qRT-PCR, indicating that our RNA-Seq data were reliable (Fig. S5).

3.4. Identification and expression pattern analysis of DEGs

The FPKM values from the different sample libraries were collected and analyzed to investigate the changes in gene expression and to identify critical genes involved in lignin synthesis in alfalfa.Respectively 8074 and 7752 unigenes were differentially expressed in the leaf and stem samples(Fig.S6).These DEGs were divided into two categories according to their expression pattern:up-regulated and down-regulated transcripts (Fig. S7). In the leaf samples, the up-regulated exceeded the down-regulated DEG,but the opposite phenomenon was observed in the stem samples.Respectively 6380 (79.02%) and 6058 (78.15%) DEGs were found to be leaf-and stem-specific(Fig.S6C).A set of 1694 unigenes were differentially expressed in both leaf and stem samples, and their expression pattern in the two sample types was similar (Fig. S8).

3.5. GO functional analysis

To further investigate whether similar mechanisms are involved in lignin synthesis in alfalfa leaf and stem, the GO database was used to classify DEG functions. In total, 60 functional groups were distributed into three categories (Fig. S9). Under the cellular component (CC), molecular function (MF), and biological process (BP) categories, ‘‘cell”, ‘‘catalytic activity”, and ‘‘metabolic process”were the largest groups in leaves and stems,respectively.Of the DEGs in leaves and stems, respectively 5084 (62.97%) and 5014 (64.68%) were successfully annotated in the GO database;1081 DEGs were detected in both tissues,and the remaining DEGs were tissue-specific, suggesting the presence of novel transcripts with tissue-specific functions.

GO enrichment (top 20) analysis of the DEGs in leaves, stems,and parts shared by the two tissues were performed with aPscore cut-off of 0.05 to identify significantly enriched GO terms among the DEGs.The results of 6 BP,11 MF,and 3 CC categories were considered significantly enriched among the DEGs in leaves(Fig. S10A). Seven BP, 9 MF, and 4 CC categories were significantly enriched in stems (Fig. S10B). Nine BP and 11 MF categories were significantly enriched in both tissues, but GO terms of CC did not appear (Fig. S10C). ‘‘Catalytic activity” was the most significantly enriched category both in leaves and stems, followed by ‘‘integral component of membrane” and ‘‘intrinsic component of membrane”. The DEGs most significantly enriched in the leaves were‘‘metabolic process”, followed by ‘‘phosphorylation” and‘‘phosphate-containing compound metabolic process” in the BP category, whereas ‘‘oxidation-reduction process”, ‘‘phosphorylation” and ‘‘carbohydrate metabolic process” were the most abundant categories in stems. Among the 1081 DEGs shared in both tissues, among the top 20 GO terms, the ‘‘catalytic activity” term was the most significantly enriched,followed by‘‘oxidation-reduc tion process” and ‘‘oxidoreductase activity” (Fig. S10C).

Fig. 5. Weighted gene co-expression network analysis (WGCNA) of common DEGs in leaves and stems.

3.6. KEGG pathway enrichment analysis

To characterize the complex biological behaviors of the transcriptome,all DEGs were subjected to a KEGG pathway enrichment analysis using the KOBAS v3.0 website. Respectively 3263 and 3216 DEGs were assigned to 125 and 128 KEGG pathways in the leaves and stems. The top 20 potential pathways of the leaves and stems were identified as significantly over-represented pathways(Fig.3).‘‘Phenylpropanoid biosynthesis”,‘‘plant hormone signal transduction”,‘‘fatty acid degradation”,and‘‘starch and sucrose metabolism” were significantly overrepresented pathways in leaves, while ‘‘starch and sucrose metabolism”, ‘‘phenylpropanoid biosynthesis”, ‘‘base excision repair”, and ‘‘DNA replication” were significantly overrepresented pathways in stems. ‘‘Phenylpropanoid biosynthesis” was the most enriched pathway in both tissues, while ‘‘plant hormone signal transduction” was significantly enriched only in leaves.

3.7. Identification of TFs

Of the 8074 leaf and 7752 stem DEGs,365 were found to belong to 42 TF families for leaves and 382 were found to belong to 39 TF families for stems(Fig.S11).Thirty-four TF families were common to the two tissues. Members of the MYB family (60 and 47 for leaves and stems, respectively) were the most abundant, followed by the AP2-EREBP (40 and 45), bHLH (33 and 48), WRKY (24 and 21), and NAC (21 and 28) families. Respectively 111 and 95 TFs were identified as specifically responsive in leaves and stems:C2C2-YABBY (3), E2F-DP (3), TifY (2), Alfin-like (1), CAMTA (1),LIM (1), RWP-RK (1), and TIG (1) were expressed specifically in leaves, and EIL (3), S1Fa-like (1), TAZ (1), Trihelix (9), TUB (3),and VARL(1) were expressed specifically in stems. The expression of 94 TFs common to leaves and stems was investigated. All of these TFs showed similar expression patterns in leaves and stems,and they could be divided into two groups: up- and downregulated (Fig. 4). Most NAC and WRKY TFs belonged to the upregulated group and most bHLH and bZIP TFs to the downregulated group.

3.8. WGCNA of common DEGs in leaves and stems

A total of 1694 DEGs were obtained from transcriptomes of alfalfa leaves and stems at five nodes. In order to further investigate the key genes of alfalfa lignin synthesis, the weighted gene co-expression network,including seven modules,was constructed using WGCNA(Fig.5).According to the color of each module in the figure, these modules were named turquoise, green, yellow, blue,brown,red,and gray,respectively.Of these seven modules,the turquoise module contains the most genes (429), followed by green(99), yellow (134), blue (397), brown (150), red (46), and gray (5)modules. The correlations between the modules and the NDF,ADF, and ADL contents ranged from -0.78 to 0.57. Among the seven modules, the yellow module had the highest correlation coefficient and showed significant (P＜0.001) associations with ADF and ADL contents, suggesting that the yellow module had the greatest relevance to lignin synthesis(Fig.S12;Table S3).Analysis of correlations between module membership (MM) and gene significance(GS) of the seven modules showed a significant correlation for the yellow module (Fig. S12B-D). Thus, the members of the yellow module may be closely associated with alfalfa nutritional quality.

Fig.6. Characterization of MsNST1.qRT-PCR analysis showed that the expression level of MsNST1 was higher in old than in young stems and leaves(A);two M.truncatula nst1 mutants(B)and two independent MsNST1 transgenics(C)were used in this study.L1 to L8 and S1 to S8 indicate leaves and stems from five different nodes(1,2,4,6,and 8,counting from the top of the stem),respectively.Asterisks above the bars in figures A and C indicate significant differences at the level of 0.05 compared with L1 and L2 and with the wild-type, respectively.

3.9. MsNST1 gene influences lignin synthesis

By qRT-PCR, the expression ofMsNST1(CL6270.Contig2_All)was higher in old than in young stems and leaves (Fig. 6A), andMsNST1is a homolog ofMtNST1(Fig.S13).These results suggested thatMsNST1might also influence the biosynthetic pathway of lignin. Accordingly,MsNST1was overexpressed inM. truncatula, and two independentMsNST1transgenics OE-7 and OE-14 were obtained.TwoM.truncatula nst1mutants(nst1-1andnst1-2)were used(Fig.6B,C).At first sight,the overall growth phenotypes of the two OE-MsNST1lines appeared similar to those of wild-type plants(Fig.7A,S14A);the stem length and diameter ofnst1mutants were smaller than those of wild-type and OE-MsNST1plants (Fig. S14B,C), but the branch number ofnst1mutants was the largest(Fig.S14D).These results indicate thatMsNST1can affect stem lignificationinM. truncatula.

The fourth internode of these fiveM. truncatulamaterials was used for histological analysis. In wild-type and OE-MsNST1plants,the G- and S-lignin were deposited in the primary xylem vessel(Xv), primary xylem (Px), and secondary xylem (Sx), but G-lignin was deposited only in the Xv, and S-lignin was not present in thenst1mutants (Fig. 7B-K). In comparison with wild-type plants,the OE-MsNST1plants clearly showed more G- and S-lignin in the Sx. The total lignin content of the stem (all internodes) in the OE-MsNST1plants was increased by around 15% compared with the wild-type, while that in thenst1mutants was decreased by around 20% (Fig. 8A). Lignin composition analysis revealed a reduction in G-lignin levels in the stems of thenst1-1mutant compared with the wild type,but a proportionately larger decrease in S-lignin such that the G/S ratio was dramatically increased in thenst1-1mutant. In contrast, a significant increase appeared in both G-and S-lignin of OE7,but there was no significant difference in G/S radio compared with the wild type (Fig. 8B-D). The relative expression levels of seven genes involved in lignin synthesis were detected, andCOMT,F5H, andHCTgenes showed significant differences amongM.truncatulawild-type,nst1mutants,and overexpression lines (Fig. 9). The nutritional value of these threeM.truncatulamaterials was measured using the near-infrared spectrometer. The NDF and ADF contents ofnst1-1were the smallest,and it also showed the highest crude protein content and RFV(Fig. S15). These results suggest thatMsNST1affects the synthesis of lignin and that down-regulating its expression may increase the nutritional value of alfalfa.

Fig. 7. MsNST1 affects lignin deposition in M. truncatula. The fourth internodes of 12-week-old WT, nst1 mutants, and OE-MsNST1 plants (A) were used to observe the characteristics of changes in lignin by histological analysis.G-lignin was stained with m-phloroglucinol ethanol and HCl(B-F),and S-lignin was stained with KMnO4,HCl,and ammonia solution(G-K).Co,collenchyma;Pf,primary phloem fiber;Pi,pith;Pp,primary phloem;Px,primary xylem;Sx,secondary xylem;Xv,primary xylem vessel.Scale bars, 1 cm.

4. Discussion

4.1. High lignification affects the nutritional quality of alfalfa

Cell walls contain large amounts of fiber, and these fiber fractions are often analytically quantified using the detergent fiber system:NDF,ADF,and ADL and used to evaluate the nutritive value of forages [43]. Generally, forages with a low nutritive value have greater NDF concentrations and more highly lignified cell walls,which reduce forage digestibility. A total of 10 samples were collected in the present study, and the NDF, ADF, and ADL contents continually increased in S1-S8,while only ADL contents gradually increased in L1-L8 (Fig. 1B).In a reported study[44], maize stems were divided into three parts:the upper(above internode 12 up to and including internode 15), middle (above internode 8 up to and including internode 11), and lower stem (above internode 5 up to and including internode 7), and the NDF, ADF, and ADL contents also increased continually from the top to the bottom of the stem.As shown in Fig.1B,the NDF,ADF,and ADL contents were higher in stems, except S1 and S2, than in leaves. Similar results also appeared in a previous study [27], and the nutritive value of two alfalfa cultivars were evaluated.NDF and ADL contents were higher in stems than in leaves, but an opposite trend was observed for crude protein.Therefore,reducing the lignification of stems should be investigated in alfalfa breeding research.

4.2. Sequence quality and qRT-PCR verification

In total, 199,583 unigenes were identified in the 30 sample libraries with a mean length of 1332 bp, longer than previously reported in alfalfa [45-47]. Of these all unigenes, 86% exhibited significant similarity (BLAST, E-value ≤10-5) to genes in at least one public database, a proportion higher than that identified in other plants, such asVicia sativa(66.10%) [48],Campeiostachys nutans(78.01%) [49], andElymus sibiricus(79.81%) [50]. The remaining unannotated unigenes may represent untranslated regions, non-coding RNAs, mis-assembly, or alfalfa-specific genes[48]. qRT-PCR verification showed that our sequencing data were reliable and could be used to identify genes that are differentially regulated in leaves and stems from nodes of alfalfa (Fig. S5).

Fig.8. Analysis of lignin content and composition in M.truncatula plants.The lignin content(A)and monomer content(B and D)in the stem material(all internodes)of five M.truncatula plants were determined,and a G/S radio was calculated(C).Asterisks above the bars indicate significant difference at P ＜0.05 compared with the wild type by Duncan’s multiple range test.

4.3. Role of plant hormones in alfalfa lignin synthesis

Plant hormones play an important role in plant growth and development. Gibberellic acid (GA), abscisic acid (ABA), and auxin(AUX) are the main endogenous hormones acting during plant development from flowering to harvest [51,52]. In the present study, the ‘‘plant hormone signal transduction” pathway was significantly enriched in alfalfa, and more so in leaves than in stems(Fig. 3). Similar results were reported in previous studies [9,53],showing that these plant hormones are indispensable signaling components during the development of leaves and stems in alfalfa.In previous studies [54,55], ABA and indole-3-acetic acid (IAA)were involved in lignification in horticultural plants. Under ABA treatment, lignin content increased and visible granulation was observed in the juice sacs of pummelo,and the expression of lignin biosynthetic gene includingCgPAL1,CgPAL2,Cg4CL1, andCgC3H,was induced, [13].PtrARF2.1is one of the poplar auxin response factor 2 members and is repressed by exogenous auxin treatment.Ectopic deposition of lignin in leaf veins and petioles was identified inPtrARF2.1-RNAi lines by histological staining[56].Among the 31 DEGs involved in plant hormone signal transduction in our study,20 were identified as auxin response genes, and their expression was higher in young than in old leaves and stems. Only one ABA receptor, PYL9-like protein (unigene11300_All), was observed,and its expression pattern was opposite to that of auxin response genes (Fig. S16). These results suggest that these auxin response genes and PYL9-like protein may be involved in signal transduction processes during the development of respectively young and old alfalfa leaves and stems.

4.4. Tfs involved in alfalfa lignin synthesis

The transcription factors involved in secondary cell wall formation have been well studied in plant species includingArabidopsis thaliana[57],Populus trichocarpa[58], andPanicum virgatum[59].However,the regulatory mechanisms of secondary cell wall formation are not well understood in alfalfa.In the past few decades,the development of transcriptome sequencing and bioinformatic pipelines has facilitated the use of large-scale data,and these tools were used to predict TFs involved in lignin biosynthesis. Seven TF families (MYB, bHLH, NAC, ERF, WRKY, C2H2, and bZIP), which coordinate their activity with lignin biosynthesis genes inA. thalianaand switchgrass,were predicted by comparative co-expression network analysis [59]. In the present study,to identify the regulatory mechanisms of lignin biosynthesis in alfalfa, 10 leaves and stems from five nodes were collected, and transcriptome analysis was performed. Respectively 42 and 39 TF families were found to be represented in leaves and stems (Fig. S11).Thirty-four TF families(94 DEGs)were common to these two tissues(Fig.4),and these TF families were also identified in a previous study [9].

Fig.9. Expression of lignin-synthesis genes in M.truncatula plants.(A-G)represent the respective expression levels of HCT,F5H,COMT,C3H,4CL,CCR,and CCoAOMT genes in three M. truncatula materials. Asterisks above bars indicate significant difference at P ＜0.05 compared with the wild type by Duncan’s multiple range test.

Two studies have found that NAC and MYB transcription factors act as first-and second-level master switches to regulate a series of downstream transcription factors and genes in the biosynthesis of the secondary cell wall[60,61].As the largest families of transcription factors in the plant kingdom,MYB proteins function in diverse biological processes during plant development [37], and are involved in regulating lignin biosynthesis[62].AtMYB58,AtMYB63,CgMYB58, andPtoMYB055are specific transcriptional activators of lignin biosynthesis inA. thaliana,C. grandis, andP. tomentosa[13,62,63],butAtMYB4,AtMYB52,andAtMYB75negatively regulate lignin biosynthesis [64-66]. However, NAC TFs function as master switches in the biosynthetic pathways for cellulose,xylan,and lignin by initiating a transcriptional signaling network that either affects MYB TFs or regulates the expression of structural genes[67].In the present study,the MYB TF family had the largest number of members in leaves (Fig. S11). Of the 94 TFs common to leaves and stems, 11 MYB and 8 NAC TFs were found, and most of these TFs were up-regulated in the old samples, such as‘‘Unigene31534_All”, ‘‘Unigene42900_All”, ‘‘CL11823.Contig2_All”,and ‘‘Unigene42901_All” (Fig. 4). These TFs may thus function in lignin synthesis regulation.

4.5. Regulatory network of lignin synthesis in alfalfa

Many lignin-synthesis genes have been reported,includingC3H,HCT,CAD,CCoAOMT, andCOMT, and these genes are recognized as being involved in phenylpropanoid metabolism [1,12,68].BpCCoAOMTis most highly expressed in young stems,and the total lignin content of antisenseBpCCoAOMTtransgenic tobacco decreased by 39% relative to control plants [69].C3HandHCTgenes are involved in the lignin synthesis pathway and the determination of lignin structure, content, and monomer ratio [70].The S unit and S/G ratio in theC3HRNAi transgenic poplar were increased, while its G unit was decreased, but the opposite trend appeared inHCT-down-regulated transgenic poplar [12]. In the present study, 3CCoAOMT, 3F5H, 3HCT, 6CAD, 6PAL, 4C4H, 10COMT, and 174CLgenes were found in leaf and stem DEGs(Table S4); most of these lignin synthesis genes were upregulated in old leaves and stems (below the 4th node) (Fig. 10),a finding consistent with the dynamic changes in ADL content among nodes (Fig. 1B). TheMsCCoAOMTandMsCOMTgenes have been shown to affect the structure of lignin and biosynthesis of G- and S-lignin [71], and a low-lignin alfalfa cultivar KK179 was cultivated by down-regulation ofCCoAOMT[28]. Our study has shown that a NAC TFMsNST1affects the synthesis of lignin inM.truncatula,including G-and S-lignin(Fig.8A,B).Compared with the wild type,the expression levels of three lignin synthesis genes were respectively reduced and increased innst1mutants and OEMsNST1plants (Fig. 9). The G/S radio was significantly increased in thenst1-1mutant relative to the wild type, while there was no significant difference in OE-MsNST1plants (Fig. 8C). As a possible explanation,the expression level ofF5Hwas reduced by around 90%in thenst1mutants relative to the wild type(Fig.9B),severely impairing the synthesis of S-lignin, and more carbon flow was provided in OE-MsNST1plants, promoting the synthesis of G- and S-lignin at the same time. These results suggest thatMsNST1is involved in the synthesis of lignin. The NDF, ADF, crude protein contents, and RFV in these three types of materials showed significant differences, but there was no significant difference in crude fat, crude ash, calcium ion, potassium ion, and phosphorus ion contents (Fig. S15). We accordingly speculate thatMsNST1is involved in carbon rather than mineral metabolism.Our study provides abundant gene resources for low-lignin alfalfa breeding research, which could be used to cultivate superior alfalfa for animal husbandry and biomass energy.

5. Conclusions

Alfalfa leaves and stems from five nodes were used as materials to investigate the molecular regulatory mechanism in lignin synthesis by RNA sequencing. Respectively 8074 and 7752 differentially expressed genes (DEGs) were identified in leaves and stems, and 1694 DEGs were common to both tissues. ‘‘Phenylpropanoid biosynthesis” was the most enriched pathway in both leaves and stems. WGCNA revealed 134 key regulatory genes for lignin synthesis of alfalfa. Respectively 365 and 382 transcription factors (TFs) were identified in leaves and stems, and 94 TFs were common to these two tissues. The NAC family transcription factorMsNST1gene showed high expression in old leaf and stem tissues.The deposition pattern of G- and S-lignin differed among the wild type,nst1mutants,and overexpression lines,and the G/S ratio was dramatically increased in thenst1-1mutant, compared with the wild-type, whereas there was no significant difference in overexpression lines. The expression of three lignin synthesis genes in these three materials differed, indicating thatMsNST1is involved in lignin synthesis. These findings provide a genetic basis and abundant gene resources for low-lignin alfalfa breeding research.

CRediT authorship contribution statement

Qiang Zhou:Investigation, Data curation, Writing - original draft.Pei Mao:Data curation,Resources.Dong Luo:Data curation,Formal analysis.Xutian Chai:Data curation, Formal analysis.Hao Deng:Data curation, Formal analysis.Qiangen Fang:Methodology.Longfa Fang:Writing - review & editing.Zhibiao Nan:Conceptualization, Resources.Jiangqi Wen:Methodology, Resources.Zhipeng Liu:Conceptualization,Writing-review&editing,Funding acquisition, Supervision.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments

This research was supported by the National Natural Science Foundation of China(32071862 and 31722055),the China Postdoctoral Science Foundation(2020M683609),and the Key Science and Technology Foundation of Gansu Province (19ZD2NA002).

Appendix A. Supplementary data

Supplementary data for this article can be found online at https://doi.org/10.1016/j.cj.2022.02.002.