The soybean dwarf mutant M₂ line was selected from EMS-induced mutants of cultivar Zhongpin 661 (Li et al., 2017; Li et al., 2018b). Plants were germinated on plate and transplanted into soil in greenhouse with the same growth conditions. Progenies of dwarf mutant and WT control were planted in Changping base of Chinese Academy of Agricultural Science in the summer of 2014. After harvesting seeds, the stems of plants were collected for cellulose content determination. Six plant individuals were phenotyped for plant height, stem, leaf size, branch and node number at different developmental stages. Fresh leaves from the third node were sampled for genome resequencing, RNA-seq analysis, and methyl-seq at V5 stage.

Genomic DNA was isolated from collected leaves of the mutant and Zhongpin 661 with a genomic DNA purification kit (Thermo Fisher Scientific Inc., United States) according to the manufacturer’s protocol. DNA samples were quantified using a Quawell Q5000 spectrophotometer (Quawell Technology, Inc., United States). The DNA from the dwarf mutant and Zhongpin661 were used to construct paired-end sequencing libraries, which were sequenced on an Illumina HiSeq-2500 platform independently. After removing adapters and low-quality sequence reads, the clean reads were further rechecked for quality with FastQC (https://www.bioinformatics.babraham.ac.uk/projects/fastqc/). High-quality reads were aligned and mapped to the Glycine max Wm82.a2.v1 reference genome (Schmutz et al., 2010) downloaded from Phytozome (http://phytozome.jgi.doe.gov/pz/portal.html#!info?alias=Org_Gmax)with BWA (version 0.7.1) with default parameters (Li and Durbin, 2009). GATK (Genome Analysis Toolkit, version 4) was used to call SNPs and small InDels (< 50 bp) across the mutant pool and wild type pool with a minimum read coverage great than two (Mckenna et al., 2010).

Around 3 μg RNA extracted by TriZol (ThermoFisher, Cat. 15596026) was used for sequencing library preparation with NEBNext^® Ultra™ RNA Library Prep Kit for Illumina^® (NEB, USA) according to the manufacturer’s recommendations. The libraries were sequenced on an Illumina HiSeq 2500 platform and 100 bp paired-end raw reads were generated.

Raw reads were firstly processed to filter out adapter, ploy-N, and low-quality reads. The clean reads with high quality were mapped to the soybean genome (Schmutz et al., 2010) with TopHat 2 (Kim et al., 2013). The read counts of each gene were adjusted with EdgeR through one scaling normalized factor (Mortazavi et al., 2008). Differential expression analysis between mutant and WT was performed with the DEGseq R package (1.12.0) (Wang et al., 2010). The P-values were adjusted using the Benjamini & Hochberg method. Corrected P-value less than 0.005 and fold changes greater than two were set as the threshold for significantly differential expression.

MethylC-seq which differentiates cytosine from 5-methylcytosine was used to examine the methylation in a DNA sample (Urich et al., 2015). DNA was extracted by CTAB and then quantified using Qubit^® DNA Assay Kit in Qubit^® 2.0 Fluorometer (Life Technologies, CA, USA). For library construction, fragmented DNA of 200-300 bp in length was obtained by sonication of a total amount of 5.2 μg genomic DNA on Covaris S220, followed by end repair and adenylation. Then these DNA fragments were treated twice with bisulfite with EZ DNA Methylation-Gold Kit™ (Zymo Research). Then the single-strand fragments were enriched by PCR amplification with KAPA HiFi HotStart Uracil + ReadyMix(2X). The libraries were sequenced on an Illumina HiSeq-2500 platform to generate 100-bp paired-end raw reads.

Raw reads were processed to remove adaptor sequence, oligo-N containing sequence (>10%), and low quality (PHRED score<=5, and percentage of the low-quality bases >=50%) to create clean reads. These clean reads were aligned to the indexed soybean reference genome (As above) using Bowtie2 (Langmead and Salzberg, 2012) and Bismark (version 0.12.5; (Krueger and Andrews, 2011). The alignments were both performed in bisulfite-converted format with C-to-T and G-to-A.

To estimate the methylation level, we identified the methylation site using s i , j + ~ B i n ( s i , j + + s i , j − , r i , j ) , in which the sum( s i , j + ) of methylated counts was modeled as a binomial (Bin) random variable with methylation rate(r _i,j ). A sliding-window approach with window size of 3,000 bp and step size of 600 bp was used to calculate the sum of methylated and unmethylated read counts (Smallwood et al., 2014). Methylation level (ML) for each C site shows the fraction of methylated Cs, and is defined as  ML(C)=reads(mC)/[reads(mC)+reads(C)] . ML was then corrected by  ML _(corrected)=(ML−r)/(1−r) , where r represents the bisulfite non-conversion rate according a previous description (Lister et al., 2013).

Gene Ontology (GO) enrichment analysis of differentially expressed genes was conducted with the GOseq R package, in which gene length bias was corrected (Young et al., 2010). Corrected P-value less than 0.05 were regarded as significantly enriched.

The cellulose associated compounds in the stem were measured for ZP661 and dwarf mutant according to the national and agriculture standards. Briefly, the crude fiber was measured with methods described in GB/T6434-2006/ISO6865:2000 (GB/T6434, 2006); the acidic detergent fiber in GB/T20805-2006 (GB/T20805, 2006) and neutral detergent fiber were measured with methods described in GB/T20806-2006 (GB/T20806, 2006).

To predict CpG island region of DEGs, genomic sequence of a given gene, spanning from promoter sequence (-2000 bp) to 3’ UTR, was extracted from the soybean genome (Wm82.a2.v1). Then CGI region in the extracted sequence was predicted with the online tool Newcpgreport (http://www.ebi.ac.uk/Tools/seqstats/emboss_newcpgreport/). Methylation site in each context within CGI region was counted for further analysis.

A soybean dwarf mutant (dm) was screened from an M2 population derived from EMS mutagenized cultivar Zhongping661 (ZP) seed lines. This dm showed a stunted growth with a much shorter plant height ( Figure 1C ), thinner stem ( Figure 1C ), and smaller leaf size ( Figure 1D ) than the control ZP. Interestingly two branches and three leaves were juxtaposed at the fourth node of dm, and the lengths of the two branches were almost the same as the main stem ( Figures 1A–C ), as compared with control with only one trifoliolate and without branch on the fourth node. The node number of M3 plants of this dwarf mutant was significantly greater than that of the control ( Figure 1E ). No other difference was found between the mutant dm and control ZP.

The whole genomic sequence of soybean ZP is not available so the generation of the DNA sequence of this cultivar is needed. To understand the DNA variation caused by EMS, DNA in the dm and control plant ZP were sequenced with an Illumina HiSeq 2500 platform to generate ~10 x average depth of short sequencing reads ( Table 1 ). We focused on SNP and InDel variants resulting from EMS mutagenesis. After mapping reads to soybean reference genome Williams82, ~210,000 SNPs in total were identified with GATK (Mckenna et al., 2010) from both the dm and control compared with Williams82, suggesting lots of differences from Williams82. However, only 32,231 SNPs (~34 kb/SNP) were presented between the dm and the control, which should come from mutagenesis. Approximate 94.6% of these SNPs remained heterozygous. Transition (A↔G or T↔C) overtook transversion (A↔T or G↔C) with Ti/Tv =1.43 in the dwarf mutant ( Supplementary Table 1 ). Relative to annotated gene models in Wm82.a2.v1, most of these SNPs (76.78%) were distributed among intergenic regions; however, only a small number of SNPs (0.29%) caused gene structure variation via synonymous or nonsynonymous substitution, transcriptional start site gained or lost, stop codon gained or lost, and splice site acceptor or donor change (first 2 bp or last 2 bp of an exon) etc. In total, these SNPs resulted in 600 nonsynonymous substitutions and 18 alternative splicing within genes ( Figure 2 ). Most of these genes were unequally scattered on 20 chromosomes, with the most mutated genes (96) on Chr03 ( Figure 2D ; Supplementary Table 2 ). It is interesting that around 11.7% of SNPs are located within 5 kb of 3’ prime downstream of genes.

A total of 2716 small InDels (< 50 bp), including 1356 insertions and 1350 deletions, were detected between the dm and the control ( Figure 2C ). The indels were distributed on 19 chromosomes except for Chr12 ( Figure 2E ). Interestingly, the utmost indels are located on Chr03, which coincides with SNPs. Most of these indels are located on intergenic regions (~59.2%). For genic indels, around 27.98% preferentially occurred in a gene’s downstream region (< 5 kb) on the 3’ prime, and ~9.5% within introns. Only ~1.4% are located within exon, which led to codon change and alternative splicing of 33 genes ( Supplementary Table 3 ).

Of these mutated genes with SNPs and indels, 15 contained two types of critical mutations of nonsynonymous substitution, alternative splicing, and codon change ( Figure 3A ). Taken together, genomic variation resulted in 636 mutated genes. Of those, 44.2% (281 out of 636 genes) genes had known GOs, and their molecular function of GOs were enriched in binding activity toward protein, ADP, nucleotide/nucleoside, purine ribonucleotide/ribonucleoside, and adenyl ribonucleotide/ribonucleoside ( Figures 3B, C ).

To obtain transcriptional regulation on the dwarf phenotype, we performed RNA-Seq analysis of transcriptomes in leaves from the third node of the dm and the control at the six unfolded trifoliolate leaves (V6) stage. Over 60 Mb clean Illumina paired-end short reads were generated from the dwarf mutant and control ZP independently. After mapping these reads to the soybean reference genome (Williams82.a2.v1), the gene’s expression level was evaluated by calculating Fragment Per Kilo bases per Million reads (FPKM) ( Table 1 ). In total, 424 down-regulated and 263 up-regulated differentially expressed genes (DEGs) were identified in the dm relative to the control with at least two folds of expressional change and corrected P-value< 0.005 ( Supplementary Figure S1 ). GO analysis revealed that these up-regulated DEGs were significantly (Fisher test P< 0.01 and FDR< 0.01) enriched in 34 GOs, which were mostly related to cell glucan, cellulose and polysaccharide metabolism, biosynthesis, and energy associated activities( Figures 4A, B ; Supplementary Table 4 ). In contrast, down-regulated DEGs were mostly enriched in 16 GO categories (Fisher test P< 0.01 and FDR< 0.01), mostly related to fatty acid synthesis Figures 4C, D ; Supplementary Table 5 ).

To evaluate epigenetic regulation, we measured the genome-wide DNA methylation status in leaves, the same material used for RNA-seq and DNA sequencing, of the dwarf mutant dm and control. In total, ~137.5 Mb and ~127.7 Mb clean reads were generated for each line, meaning ~16 x coverage of soybean genome on average ( Table 1 ; Supplementary Figure S2 ). Comparison analysis found that the DNA methylation level in the dwarf mutant was slightly lower than that in control although ~11% of all cytosines in genome had been methylated in both samples ( Figure 5 ). Generally, ~57.6% (C:58.96%, M:56.24%), ~32.8% (C:34.12%, M:31.47%) and ~2.6% (C:2.48%, M:2.73%) of cytosines were methylated in CG, CHG and CHH context, respectively ( Figure 5 ). The dwarf mutant had slightly fewer methylated cytosines at CG and CHG sites but more at CHH than control, in terms of both methylation levels or relative percentage of methylation in each context ( Figure 5 ). A total of 704 differentially methylated regions (DMR) spanning from 45 bp to 21.89 kb were identified between dm and control, among which 119 and 343 genes were predicted on chromosomes and unanchored scaffolds, respectively. The unanchored scaffolds are mostly repeat sequence-like transposons, so those genes will most likely not affect the dwarf phenotype. These 119 chromosomal genes’ transcription could be affected by methylation alternation, theoretically inhibited, which we further investigated.

To better understand the crossing molecular regulation in the dwarf mutant dm, we further checked the genes affected by any two combinations of variants, transcriptome, and methylome. Eleven genes, ~1.6% of DEGs, were found to have both nonsynonymous mutations and differential expression, which suggested that these genes’ expression was most likely directly regulated by mutations ( Figure 6A ). Thus, these genes were worthy of further investigation. The remaining 675 DEGs (98%), except one DEG, didn’t have nonsynonymous mutation and didn’t belong to DMR; therefore, they may be indirectly regulated by others, e.g., the cascaded signal of mutated genes. Three DMRs were found to also contain genes with nonsynonymous mutation(s), which were also important candidates for this dwarf phenotype ( Figure 6A ). All these mentioned genes were further checked in terms of expression level and functional prediction. Among these, a cellulose synthase gene (Glyma.15G157100) was included, which is known to catalyze beta-1,4-glucan cellulose microfibril for the secondary cell wall formation in the plant cellulose biosynthesis and in glycan metabolism (Mcfarlane et al., 2014). This gene expression in the dwarf mutant was increased by 3.5 folds ( Table 2 ). However, a previous study showed that a mutation in this gene reduced cellulose levels and caused a dwarf phenotype (Holland et al., 2000). So a further chemical analysis of cellulose related compounds was conducted and revealed that acidic detergent fiber, neutral detergent fiber, crude fiber, cellulose and lignin in the stem of the dwarf mutant were 19.6%, 13.3%, 2.5%, 20% and 17.3%, and were higher than that of control, respectively, but hemicellulose content was 2.8% lower than that in control ( Figure 6B ). We hypothesized that the mutation may affect the cellulose balance in dm, which led to dwarfism.

Methylation of CpG island (CGI) affects the targeting capability of transcriptional factors, which can change gene expression (Candaele et al., 2014; Zhang et al., 2014b). To better elucidate the possible correlation between CGI methylation and gene expression level in our study, CGI was predicted for each DEG (-2000 bp to 3’UTR) with Newcpgreport. CGI was found in 71 up-regulated and 105 down-regulated DEGs, of which ~33.7% (32) and ~66.3% (63) showed opposite methylation patterns ( Figure 7A ), namely up-regulated genes with down methylation level and vice versa. These results suggested that these DE genes might be regulated by methylation status alteration on CGI. CHH was the most abundant methylation within CGI regions followed by CG and CHG ( Figure 7B ).

Within these CGI-containing DEGs, Glyma.13G312900, encoding a BEL-1 protein, was modified at all three levels of transcript, methylation, and variants. BEL-1 was previously reported to interact with STM, ATH1, PNY, and PNF protein in the shoot apical meristem organogenesis and differentiation (Rutjens et al., 2010). The CGI sites in exon 1 of this gene were down methylated in all three contexts (CG, CHG, and CHH, Figure 8 ) and transcriptionally up-regulated in dwarf mutant. At the variant level, three nonsynonymous mutations (Ser107Gly, Ala111Thr, and Thr112Ala) were identified in this gene in the dwarf mutant relative to control, but they were located outside of the conserved domains of POX domain and Homebox KN domain ( Figure 8 ). However, the function of BEL-1 still could be impaired. These three nonsynonymous mutations were heterozygous in control and homozygous in the dwarf mutant, which may be recessive mutations. Recessive mutations had a higher chance of loss of function than homozygous mutations. Together, these results suggested that Glyma.13G312900 may cause the dwarf mutant phenotype and is of high priority in our future investigation.