MtBZR1 Plays an Important Role in Nodule Development in Medicago truncatula

Brassinosteroid (BR) is an essential hormone in plant growth and development. The BR signaling pathway was extensively studied, in which BRASSINAZOLE RESISTANT 1 (BZR1) functions as a key regulator. Here, we carried out a functional study of the homolog of BZR1 in Medicago truncatula R108, whose expression was induced in nodules upon Sinorhizobium meliloti 1021 inoculation. We identified a loss-of-function mutant mtbzr1-1 and generated 35S:MtBZR1 transgenic lines for further analysis at the genetic level. Both the mutant and the overexpression lines of MtBZR1 showed no obvious phenotypic changes under normal growth conditions. After S. meliloti 1021 inoculation, however, the shoot and root dry mass was reduced in mtbzr1-1 compared with the wild type, caused by partially impaired nodule development. The transcriptomic analysis identified 1319 differentially expressed genes in mtbzr1-1 compared with wild type, many of which are involved in nodule development and secondary metabolite biosynthesis. Our results demonstrate the role of MtBZR1 in nodule development in M. truncatula, shedding light on the potential role of BR in legume–rhizobium symbiosis.


Introduction
Medicago truncatula is a model fabacean species for nodulation and nitrogen fixation. During the initiation of symbiotic nodules, for example, the root nodule microsymbionts (collectively termed rhizobia) secrete nodulation (Nod) factors, which are lipochito-oligosaccharides functioning as signal molecules, in response to flavonoids from the host plant. The Nod factors are then perceived by the corresponding receptor kinases in the root tissue, which in turn trigger nodule organogenesis and root hair deformation [1]. After nodule formation, the host plant generates a class of nodule-specific cysteine-rich (NCR) peptides to induce the differentiation of endosymbionts into enlarged and polyploid bacteroids [2]. There are about 600 MtNCR genes in the M. truncatula genome, indicating their important roles in nitrogen fixation [3]. Hormones also play important roles during nodule development [4,5].
Flavonoids are known as secondary metabolites that play important roles in legume-rhizobium symbiosis [6]. Based on the molecular structure, flavonoids consist of flavones, flavonols, isoflavones, flavanones, flavanols, and anthocyanins [7]. The flavonoid biosynthesis pathway starts from downstream of the phenylalanine pathway, in which p-coumaroyl co-enzyme A (CoA) is produced and catalyzed into chalcone by chalcone synthase. In the flavonoid biosynthesis pathway, chalcone is catalyzed into dihydrokaempferol by chalcone isomerase and flavanone 3-hydroxylase, and then by flavonol synthases into flavonols. Flavonoids production is suggested to be highly localized at infection sites in legume [8].
Brassinosteroid (BR) is an essential plant hormone in plant growth and developmental processes, such as cell elongation, photomorphogenesis, seed germination, stomata differentiation and movement, flowering, and plant response to biotic and abiotic stresses [9][10][11]. AtBZR1 (Arabidopsis thaliana BRASSINAZOLE-RESISTANT 1) is a key transcriptional regulator in the BR signaling pathway [12][13][14]. AtBZR1 has a homolog, AtBES1 (At BRASSINOSTEROID INSENSITIVE 1-EMS SUPPRESSOR 1), which shares 88% similarity in the A. thaliana genome [15]. Because of the redundancy with AtBES1, the atbzr1 mutant showed no obvious phenotypic change, while the gain-of-function mutant atbzr1-D showed enhanced BR signaling in A. thaliana. Recently, mtbri1 was reported to show defective nodule development, indicating that BR signaling is involved in legume nodulation [16]. As MtBZR1 is an important transcription factor downstream of MtBRI1, it is intriguing to explore its role in nodule development in M. truncatula.
This study was undertaken with the aim to determine the function of MtBZR1 in M. truncatula R108 nodule development. In order to achieve this, we isolated mutants of MtBZR1 and characterized their phenotypes. Next, we performed transcriptomic analysis of the mtbzr1-1 null mutant, in comparison to the wild-type plants, to identify the deregulated genes. The results indicate that MtBZR1 is involved in the control of nodule development, but not its initiation, in M. truncatula. MtBZR1 was induced by the rhizobium Sinorhizobium meliloti 1021, and mtbzr1-1 displayed partially impaired nodule development compared to the control R108. In addition, transcriptomic analysis found that the genes involved in nodule development and secondary metabolite biosynthesis were under transcriptional reprogramming, consistent with the morphological changes in mtbzr1-1. These results indicate that MtBZR1 plays a role in controlling nodule development after the initiation in M. truncatula.

Search for AtBZR1 Homologs in M. truncatula
We searched the annotated M. truncatula genome for the homolog of A. thaliana AtBZR1. Medtr3g111680 encodes a 315-amino-acid (aa) protein that shares 64% identify with the A. thaliana AtBZR1 protein. According to the phylogenic analysis of 11 A. thaliana and M. truncatula homologs, Medtr3g111680 is the closest homolog of AtBZR1 in M. truncatula ( Figure 1A; Figure S1, Supplementary Materials). It is clustered with two other proteins, Medtr1g021990 and AtBZR2. In addition to the high similarity in the protein sequence, Medtr3g111680 also possesses a BZR domain with conservative sites ( Figure S2, Supplementary Materials). We, therefore, denoted Medtr3g111680 as MtBZR1.

Induction of MtBZR1 upon S. meliloti 1021 Inoculation
To characterize its role in plant development, an expression assay of MtBZR1 by quantitative real-time PCR (qRT-PCR) was carried out in various tissue types. As shown in Figure 1B, MtBZR1 expression was highest in the root tissue, and comparable in stem and leaf, and relatively low in flower and pod. As nodulation is a specific characteristic of legume, we then performed another expression assay in nodules of the M. truncatula R108. As shown in Figure 1C, MtBZR1 was induced more than 20-fold in nodules at seven and 14 days post inoculation (dpi) by S. meliloti 1021 strain, and then reduced to seven-fold at 21 dpi, compared with that at zero dpi in R108. These results suggest that MtBZR1 may be involved in nodule development in M. truncatula. . The expression at zero dpi was set as 1.0. In (B,C), MtUbiquitin gene was used as an internal control. Bars represent standard errors from three biological replicates.

Identification of the mtbzr1-1 Null Mutant
To illustrate the biological function of MtBZR1, we adopted a reverse genetic approach by searching the Noble Research Institute Mutant collection for the loss-of-function mutants. One mutant line (mtbzr1-1) that contains a single insertion in MtBZR1′s second exon was found and used for the following genetic analysis (Figure 2A). The insertion was confirmed by genotyping using the gene-specific primers and a Tnt1 primer ( Figure 2B,C), and homozygous plants were identified. As shown in Figure 2D, we detected no expression of MtBZR1 in the mtbzr1-1 mutant plant. There was no obvious difference in developmental phenotypes between mtbzr1-1 and wild-type plants ( Figure 2E).

Identification of the mtbzr1-1 Null Mutant
To illustrate the biological function of MtBZR1, we adopted a reverse genetic approach by searching the Noble Research Institute Mutant collection for the loss-of-function mutants. One mutant line (mtbzr1-1) that contains a single insertion in MtBZR1 s second exon was found and used for the following genetic analysis (Figure 2A). The insertion was confirmed by genotyping using the gene-specific primers and a Tnt1 primer ( Figure 2B,C), and homozygous plants were identified. As shown in Figure 2D, we detected no expression of MtBZR1 in the mtbzr1-1 mutant plant. There was no obvious difference in developmental phenotypes between mtbzr1-1 and wild-type plants ( Figure 2E).

Loss of Function in MtBZR1 Leads to a Decrease in Plant Dry Mass in S. meliloti 1021 Nodulation
The induction of MtBZR1 in nodulation implies its role in nodule development. Thus, we compared the growth between mtbzr1-1 mutant and wild type after S. meliloti 1021 inoculation. As shown in Figure 3A, the inoculated mtbzr1-1 mutant exhibited a smaller plant size compared with the wild-type R108. To confirm this observation, we measured the dry mass of the shoot and root tissues in mtbzr1-1 and wild-type R108 plants 21 days after inoculation. The dry mass of both shoot and root in mtbzr1-1 was significantly decreased, compared with that of the wild-type R108 control ( Figure 3B,C).

Loss of Function in MtBZR1 Leads to a Decrease in Plant Dry Mass in S. meliloti 1021 Nodulation
The induction of MtBZR1 in nodulation implies its role in nodule development. Thus, we compared the growth between mtbzr1-1 mutant and wild type after S. meliloti 1021 inoculation. As shown in Figure 3A, the inoculated mtbzr1-1 mutant exhibited a smaller plant size compared with the wild-type R108. To confirm this observation, we measured the dry mass of the shoot and root tissues in mtbzr1-1 and wild-type R108 plants 21 days after inoculation. The dry mass of both shoot and root in mtbzr1-1 was significantly decreased, compared with that of the wild-type R108 control ( Figure 3B,C).
To explore the reason why mtbzr1-1 growth was retarded, we compared the nodule development in the mutant with the wild-type plants. As shown in Figure 4A, the nodules on mtbzr1-1 roots were smaller than those from the wild type. The measurement of nodule dry mass confirmed the observation that the nodules from mtbzr1-1 were significantly smaller/lighter than those from the wild type, although the number of nodules was comparable between them ( Figure 4B,C). The genetic assay of mtbzr1-1 demonstrates that the loss of function in MtBZR1 leads to partially impaired nodule development in M. truncatula.
To further investigate the role of MtBZR1, the full length of the MtBZR1 coding sequence under control of the CaMV 35S promoter was introduced into the wild-type background. A total of six independent transgenic lines were generated and two lines with the highest MtBZR1 expression were used for further analysis ( Figure S3A compared the growth between mtbzr1-1 mutant and wild type after S. meliloti 1021 inoculation. As shown in Figure 3A, the inoculated mtbzr1-1 mutant exhibited a smaller plant size compared with the wild-type R108. To confirm this observation, we measured the dry mass of the shoot and root tissues in mtbzr1-1 and wild-type R108 plants 21 days after inoculation. The dry mass of both shoot and root in mtbzr1-1 was significantly decreased, compared with that of the wild-type R108 control ( Figure 3B,C).  To explore the reason why mtbzr1-1 growth was retarded, we compared the nodule development in the mutant with the wild-type plants. As shown in Figure 4A, the nodules on mtbzr1-1 roots were smaller than those from the wild type. The measurement of nodule dry mass confirmed the observation that the nodules from mtbzr1-1 were significantly smaller/lighter than those from the wild type, although the number of nodules was comparable between them ( Figure  4B,C). The genetic assay of mtbzr1-1 demonstrates that the loss of function in MtBZR1 leads to partially impaired nodule development in M. truncatula.
To further investigate the role of MtBZR1, the full length of the MtBZR1 coding sequence under control of the CaMV 35S promoter was introduced into the wild-type background. A total of six independent transgenic lines were generated and two lines with the highest MtBZR1 expression were used for further analysis ( Figure S3A When we carried out rhizobial inoculation of the S. meliloti 1021 strain on the 35S:MtBZR1 lines, as shown in Figure 3A, the 35S:MtBZR1#1 plant showed a similar stature as the wild-type R108 plant, and also had a similar shoot and root dry mass as wild-type R108 ( Figure 3B,C). The phenotype and dry mass of nodules were also comparable between the 35S:MtBZR1#1 plant and the wild-type control ( Figure 4A-C).  When we carried out rhizobial inoculation of the S. meliloti 1021 strain on the 35S:MtBZR1 lines, as shown in Figure 3A, the 35S:MtBZR1#1 plant showed a similar stature as the wild-type R108 plant, and also had a similar shoot and root dry mass as wild-type R108 ( Figure 3B,C). The phenotype and dry mass of nodules were also comparable between the 35S:MtBZR1#1 plant and the wild-type control ( Figure 4A-C).

Transcriptomic Profiles of Nodules in mtbzr1-1
To understand the involvement of MtBZR1 in nodule development, we performed a transcriptomic assay of mtbzr1-1. The nodules were harvested from three biological replicates of mtbzr1-1 and wild-type plants 21 days after rhizobial inoculation by S. meliloti 1021 strain. For each replicate, RNA sequencing generated more than 20 million raw reads, which were filtered and aligned against the annotated reference genome for gene quantification. The correlation efficiency between biological replicates was as high as 99.9%, indicating a good repeatability ( Figure S4A, Supplementary Materials). The differentially expressed genes (DEGs) were identified between the mtbzr1-1 and wild-type sample by DEGseq. Genes with more than two-fold expression changes and a false discovery rate (FDR) smaller than 0.001 were identified as DEGs. A total of 1319 DEGs were found in mtbzr1-1, among which 618 were elevated and 701 were suppressed ( Figure 5; Figure S4B and

Transcriptomic Profiles of Nodules in mtbzr1-1
To understand the involvement of MtBZR1 in nodule development, we performed a transcriptomic assay of mtbzr1-1. The nodules were harvested from three biological replicates of mtbzr1-1 and wild-type plants 21 days after rhizobial inoculation by S. meliloti 1021 strain. For each replicate, RNA sequencing generated more than 20 million raw reads, which were filtered and aligned against the annotated reference genome for gene quantification. The correlation efficiency between biological replicates was as high as 99.9%, indicating a good repeatability ( Figure S4A, Supplementary Materials). The differentially expressed genes (DEGs) were identified between the mtbzr1-1 and wild-type sample by DEGseq. Genes with more than two-fold expression changes and a false discovery rate (FDR) smaller than 0.001 were identified as DEGs. A total of 1,319 DEGs were found in mtbzr1-1, among which 618 were elevated and 701 were suppressed ( Figure 5; Figure S4B and Table S2, Supplementary Materials).

Discussion
In this study, we characterized MtBZR1, a transcriptional regulator of BR signaling in M. truncatula, via genetic and molecular approaches. MtBZR1 expression is elevated in young nodules (seven and 14 dpi) in comparison to roots. Furthermore, the dry mass of shoots and roots was reduced in mtbzr1-1 compared with wild type. Considering the smaller nodules in mtbzr1-1 and suppression of genes involved in nodulation, our results demonstrate that MtBZR1 plays an important role in legume-rhizobia symbiosis.
NCRs play important roles in nodulation in fabaceans [3]. We found that 41 of 43 differential expressed NCR genes were under suppression in mtbzr1-1 compared with wild type. The NCR genes are important for the infection, and most of these genes are expressed exclusively in nodules and through the processes of nodule development [17,18]. In our study, the suppression of these NCR genes is consistent with nodule dry mass reduction of the S. meliloti 1021 strain in mtbzr1-1. Within the promoter regions of NCRs, conserved motifs such as auxin response factor binding sites were reported [3]. Auxin response factor (ARF) binding sites (TGTCTC) are found to be located within 1000-bp promoter regions after searching for conserved cis-elements and motifs in NCR gene promoters [19]. The conserved ARF binding motifs, which are also found in the promoters of NCR downregulated in mtbzr1-1 in this study, suggest the important role of auxin during the symbiosis process between M. truncatula and S. meliloti 1021. In the pathway of auxin signal transduction, the auxin receptor TIR1 (TRANSPORT INHIBITOR RESPONSE 1) and the series of downstream ARF-related transcription factors are target genes of AtBZR1 in A. thaliana [20,21]. Hence, the nodulation-deficient phenotype of mtbzr1-1 is probably due to MtBZR1 s regulation of the elements downstream of other signaling pathways.
Flavonoids function as signaling molecules in legume-rhizobia recognition. Specific flavonoids are produced to activate Nod production of rhizobia's compatible symbionts in conditions with low nitrogen. In accordance with this, the chalcone synthase knockdown in M. truncatula disrupts the initiation of nodulation, although the deformation of root hair remains [22]. A diverse group of metabolites derived from the flavonoid biosynthesis pathway accumulate in host roots and nodules. Despite this, only a small number of them clearly demonstrated their role in rhizobial infection. One of these flavonoids is methoxychalcone, which significantly enhances inducing activity on the nodulation genes [8]. It is produced from coumaroyl-CoA by chalcone synthases, whose coding genes are downregulated in the mutant mtbzr1-1, and then by chalcone-O-methyltransferases. The flavonoids are related to nodulation gene expression in nodules. For example, the infection zones of nodules contain a high level of flavonoids and the nodulation gene expression levels are high [6]. The dysregulation of flavonoid production would impair nodule initiation development in the mtbzr1-1 mutant. Our data suggest that the suppression of NCR genes and flavonoid metabolism probably lead to the impairment of the development of nodules in the cross-talk pathway in mtbzr1-1.
The mutant of the BR 5α-reductase gene in the BR synthesis pathway in pea displayed defective root development and exhibited reduced nodule organogenesis [5]. Furthermore, mtbzr1-1 developed normal shoots and roots probably due to gene redundancy with its homologs in the genome, similarly to atbzr1 in A. thaliana. Interestingly, mtbzr1-1 showed defects in nodule development, indicating its different role between shoot development, root development, and nodulation. In the KEGG analysis, the MtFLSs were downregulated in mtbzr1-1. Their homolog in A. thaliana, At5g05600, functions as a flavonol synthase and a 2-oxoglutarate/Fe (II)-dependent oxygenase that hydroxylates jasmonic acid (JA) [23]. Consist with this, we found that genes involved in linoleic acid metabolism (ko00591) were enriched with DEGs. Linoleic acids are precursors for JA, which inhibits the calcium spiking induced by Nod factors [23,24]. As an important transcription factor in BR signaling, MtBZR1 might function in hormone cross-talk in nodule development.
In previous studies, the strain S. meliloti 1021 was proven inefficient in nitrogen fixation on A17 and R108 plants [25,26], and its inefficiency might affect the quantitative differences in phenotypes investigated in experiments reported here. However, the mtbzr1-1 mutant displayed significantly less nodule mass although the numbers of nodules were comparable with the wild-type plants. Despite the possible defective differentiation of bacteroids, the differences in nodule phenotypes between mtbzr1-1 and wild type were caused by the loss of function of MtBZR1. These results indicate that MtBZR1 is involved in the nodule development after its initiation. Inoculation with an efficient strain, e.g., S. meliloti 102F34, would likely result in more pronounced variation in nodule phenotype, which will facilitate screening for new mutants in the future.

Plant Materials and Growth Conditions
The plant materials used in this work were the Medicago truncatula ecotype R108. The mtbzr1-1 mutant was requested from the Noble Research Institute Tnt1 Mutant collection [27]. The mtbzr1-1 mutants were identified using gene-specific and Tnt1-specific primers (listed in Table S1, Supplementary Materials). The seeds of wild type, mtbzr1-1, and 35S:MtBZR1 were scarified mechanically and vernalized at 4 • C on moist filter paper for one week before transfer to soil. The plants were grown in a growth chamber at 22 • C under a 16-h light/8-h dark cycle with 150 µmol/m 2 /s light intensity and 70-80% relative humidity.

Root Nodule Induction
For rhizobial inoculation, wild-type, mtbzr1-1, and 35S:MtBZR1 seedlings were transferred into plastic trays filled with pearlite/sand (in a 3:1 ratio). The plants were grown in the growth chamber and watered every three days with nutrient solution without nitrogen as previously reported [28]. The rhizobium S. meliloti 1021 strain was cultured in TY (tryptone, yeast extract, and sodium chloride) medium (supplemented with 10 mg·mL −1 tetracycline, 200 mg·mL −1 streptomycin, and 6 mmol·L −1 calcium chloride) at 28 • C in a shaker to reach an OD 600 value of 1.0 [28]. Five-day-old seedling of wild type, mtbzr1-1, and 35S:MtBZR1 were inoculated with 5 mL of rhizobial suspension diluted to OD 600 = 0.1. Nodule numbers and dry masses were determined 21 days after inoculation.

Identification of MtBZR1 Gene and Phylogenetic Analysis
To identify the homologs of A. thaliana BZR1 in M. truncatula, we used the amino-acid sequences of the two AtBZRs (AtBZR1 and AtBZR2) and four AtBZR1 homologs (AtBEH1 to AtBEH4) from A. thaliana for a BLASTP search in the M. truncatula (Available online: http://bioinfo3.noble.org/ doblast/). A total of seven homologous sequences, Medtr3g111680, Medtr5g019550, Medtr1g021990, Medtr2g075410, Medtr2g075400, Medtr5g026210, and Medtr7g077410, were obtained. To study the phylogenetic relationships, the six amino-acid sequences from A. thaliana and seven from M. truncatula were aligned in the ClustalW program (Available online: https://www.genome.jp/tools-bin/clustalw). The MYB transcription factors related to AtMYB5, AtTT2, MtMYB5 (Medtr3g083540), and MtMYB14 (Medtr4g125520) acted as the outgroup sequences. Then, the phylogenic tree was constructed by the neighbor-joining, maximum-parsimony, or minimum-likelihood method with 1000 bootstrap replicates in MEGA6 [29].

Plasmid Construction and Plant Transformation
The coding sequence (CDS) of MtBZR1 was amplified from a complementary DNA (cDNA) sample prepared from wild-type root RNA with MtBZR1F and MtBZR1R primers (listed in Table  S1, Supplementary Materials). The PCR product was purified and cloned into the pENTR/D-TOPO vector (Invitrogen, Carlsbad, CA, USA) and MtBZR1 s CDS was confirmed by Sanger sequencing. The resulting pENTR-MtBZR1 plasmid was used to construct 35S:MtBZR1-GFP for plant transformation by recombination with a destination vector pEarleyGate103 (Invitrogen, Carlsbad, CA, USA) [30]. For M. truncatula transformation, the 35S:MtBZR1-GFP plasmid was introduced into Agrobacterium strain EHA105 and used to transform the wild-type plant by the leaf disc method as previously described [31]. The transgenic plants were genotyped with the 35S primer and reverse primer of the gene. The positive ones were accessed for the expression level of MtBZR1 with its qRT primers for a step further.
The 35S:MtBZR1-GFP construct or the 35S:GFP control plasmid were introduced into Agrobacterium tumefaciens strain EHA105 and used for an infiltration assay of Nicotiana benthamiana leaves as previously reported [32]. Leaves from four-week-old N. benthamiana plants were infiltrated with A. tumefaciens at an OD 600 value of 0.4 containing 35S:MtBZR1-GFP or 35S:GFP, respectively, and incubated in the dark at 25 • C for 16 h. Leaf discs infiltrated with A. tumefaciens were observed using a Zeiss LSM880 confocal microscope at 488-nm wavelength for the GFP fluorescence (Zeiss, Oberkochen, Germany).

RNA Isolation and Quantitative Real-Time PCR (qRT-PCR)
To analyze the expression level of MtBZR1 in different organs, samples of root, stem, leaf, flower, and pod were harvested from two-month-old wild-type plants. For time course analysis of MtBZR1 expression after rhizobial inoculation, nodule samples from wild-type plants were harvested at seven, 14, and 21 days post-inoculation (dpi). To analyze the expression level of MtBZR1 in mtbzr1-1 and 35S:MtBZR1 plants, leaves were harvested from one-month-old wild-type, mtbzr1-1, and 35S:MtBZR1 plants. For expression analysis of NCR genes in mtbzr1-1 nodules, nodule samples from wild type and mtbzr1-1 were harvested at 21 dpi. Each biological replicate included 100 mg of nodules collected from 25 plants. To minimize the differences caused by nodule age, we pooled 100 mg of nodules from those plants for RNA extraction. All samples were frozen immediately in liquid nitrogen. RNA was extracted using RNeasy Mini Kit (Qiagen, Beijing, China), and treated with DNase I (Qiagen, Beijing, China). The quantity and quality of extracted RNA were measured by a Nanodrop 2000 Spectrophotometer (NanoDrop Technologies, Waltham, MA, USA). One microgram of total RNA was transcribed into cDNA using ThermoScript RT-PCR system (ThermoFisher, Waltham, MA, USA). After dilution (in a 1:50 ratio), the cDNA samples were used for qRT-PCR on a CFX Connect Real-Time PCR Detection system (Bio-Rad, Hercules, CA, USA) using SYBR Green PCR Master Mix (Roche, Basel, Switzerland) in three biological replicates. An MtUbiquitin gene, Medtr3g110110 [33], was used as an internal control, and relative gene expression was calculated according to the ∆∆C T method [34]. The primer sequences used for qRT-PCR were designed using Primer Express 3.0 software (https://primer-express.updatestar.com/) and are listed in Table S1 (Supplementary Materials).

Transcriptomic Analysis
For the transcriptomic analysis of mtbzr1-1, three biological replicates of nodule tissue were harvested from wild-type or mtbzr1-1 plants 21 days after rhizobial inoculation. RNA samples were sequenced on a BGISEQ-500 platform as eukaryote samples following the manufacturer's instructions at BGI Genomics Institute (BGI-Shenzhen, Shenzhen, China). About 22 million raw reads were generated for each sample. Adapter sequences, low-quality reads, and reads containing more than 5% unknown nucleotides were first filtered out by trimmonmatic (v0.37) [35]. Clean reads were then aligned to the M. truncatula reference transcriptome (version 4.0) using Bowtie [36]. Gene expression (Table S5, Supplementary Materials) was quantified using RSEM (RNA-Seq by Expectation Maximization) [37], and differentially expressed genes with fold change >2 and false discovery rate (FDR) <0.001 were identified using R package DEGseq [38]. The enrichment of Gene Ontology (GO) terms and KEGG pathways was assessed using a hypergeometric test with p-values adjusted by the FDR method.

Statistical Analysis
Data processing and statistical analysis were performed using Microsoft ® Office Excel 2010 (Available online: http://www.microsoft.com/). A Student's t-test was used to evaluate if the differences were significant in the analysis of gene expression level, nodule numbers, and dry masses between samples.

Conclusions
In this study, we characterized MtBZR1, the homolog of A. thaliana BZR1 in M. truncatula, and found evidence that MtBZR1 plays an essential role in nodule development. The transcriptomic analysis of the loss-of-function mtbzr1-1 mutant suggests that MtBZR1 regulates NCR and secondary metabolite biosynthesis possibly via hormone signaling cross-talk. This work demonstrates the biological function of MtBZR1 and a potential role of BR in nodule development, which implies sophisticated hormone cross-talk during nodulation in M. truncatula.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/20/12/ 2941/s1. Figure S1.  Figure  S4. Repeatability and differential analysis of RNA-sequencing data. (A) The heatmap of Pearson's coefficient between sequencing samples. (B) The volcano plot of differential analysis between mtbzr1-1 and wild-type plants. The induced and suppressed genes are indicated as red and green points, respectively. Figure S5. The expressions of MtNCRs, MtNCR301, MtNCR237, MtNCR327, MtNCR145, MtNCR332, and MtNCR224 in mtbzr1-1 and wild-type plants. MtUbiquitin was used as an internal control. Bars represent standard errors from three biological replicates. Figure S6. The phenylpropanoid biosynthesis pathway with induced and suppressed genes indicated in red and green frames, respectively. Figure S7. The flavonoid biosynthesis pathway with induced and suppressed genes indicated in red and green frames, respectively. Figure S8. The isoflavonoid biosynthesis pathway with induced and suppressed genes indicated in red and green frames, respectively. Figure S9. The expressions of MtFLSs and MtIOMTs in mtbzr1-1 and wild-type plants. MtUbiquitin was used as an internal control. Bars represent standard errors from three biological replicates. Table S1. Primers used in this study. Table S2. The 1319 differentially expressed genes in mtbzr1-1. Table S3. The expression of genes involved in the flavonoid biosynthesis pathway. Table S4. The expression of genes involved in the isoflavonoid biosynthesis pathway. Table S5. The expression of all genes in WT and mtbzr1-1.

RSEM
RNA -seq by expectation maximization  TIR1  TRANSPORT INHIBITOR RESPONSE 1  TT2  TRANSPARENT TESTA 2  TY tryptone, yeast extract, and sodium chloride WT wild type