A chromosome-scale genome assembly of mungbean (Vigna radiata)

Sample collection and DNA extraction

We selected the mungbean cultivar KUML4 for genome sequencing and assembly. Fresh leaves were collected from a 4-week-old mungbean plant. The leaf samples were immediately frozen in liquid nitrogen and stored until use. High molecular weight (HMW) DNA was extracted from the frozen tissues using the Genomic-tip 100/G DNA preparation kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. After extraction, the concentration, integrity, and purity of the DNA were determined using the dsDNA HS assay on a Qubit fluorometer (Thermo Fisher Scientific, Waltham, USA) and pulsed-field gel electrophoresis (PFGE).

stLFR library construction and sequencing

The stLFR sequencing library was constructed using a total of 10 ng of HMW DNA following the MGIEasy stLFR Library Prep Kit’s instructions (MGI Tech, Shenzhen, China). In brief, the HMW DNA was first inserted with transposons. The transposon-inserted DNA was then incubated with clonal barcoded beads and subjected to random priming. Following that, the captured DNA molecules were disrupted to yield sub-fragments, each less than 1 kb in size. Adapters were ligated to the DNA fragments, and the library was developed by PCR. Finally, the library was sequenced on the DNBSEQ-G400 using the MGISEQ-2000RS Sequencing Flow Cell v3.0 (MGI Tech, Shenzhen, China).

Hi-C library preparation and sequencing

A high-throughput chromosome conformation capture (Hi-C) technique (Putnam et al., 2016) was performed by Biomarker Technologies Corporation (Beijing, China) to construct the chromosome-level assembly of mungbean. Briefly, the chromatin was crosslinked with fresh formaldehyde, and the fixed chromatin was then digested with HindIII restriction endonuclease. Subsequently, the digested DNA and the overhanging 5′ ends of the DNA fragments were repaired with biotinylated nucleotides. DNA fragments were then circularized by blunt-end ligation. After the cross-linking was reversed, the DNA was purified and sheared into smaller fragments of 300–700 bp in size, which were then pulled down with streptavidin beads. The purified fragments were used to generate a sequencing library that was then sequenced on the Illumina platform (PE150; Illumina, San Diego, CA, USA).

Genome assembly and Hi-C scaffolding

The stLFR reads were used to construct the preliminary draft genome assembly following the single-tube long fragment read data analysis software stLFRdenovo version 1.0.5, with default settings (https://github.com/BGI-biotools/stLFRdenovo/releases/tag/v1.0.5, accessed on 12 January 2024). This preliminary assembly was further scaffolded into a chromosome-level assembly using the high-throughput chromosome conformation capture (Hi-C) technique. The preliminary draft genome assembly and Hi-C reads were used as the input data for HiRise version 2.1.9, a software pipeline designed specifically for using proximity ligation data to scaffold genome assemblies (Putnam et al., 2016). The Hi-C sequences were mapped to the draft input assembly using BWA version 0.7.17 (Li & Durbin, 2009) (https://github.com/lh3/bwa). HiC-Pro version 2.10.0 (Servant et al., 2015) was then applied for data filtering and quality assessment of the Hi-C reads. The separations of Hi-C read pairs mapped within draft scaffolds were analyzed by HiRise to produce a likelihood model for genomic distance between read pairs, and the model was used to identify and break putative misjoins, score prospective joins, and make joins above a threshold. Manual adjustment and curation were carried out using a Hi-C visualization tool to correct potential misplacements or orientation errors. The mungbean genome assembly was deposited in National Center for Biotechnology Information (NCBI) under the accession number JAVCZA000000000. The raw reads stLFR data were deposited in the SRA database at the NCBI with the accession number SRR29729411.

Genome size estimation

To estimate the genome size, k-mer analysis was performed on the high-quality paired-end (PE) read distribution using Jellyfish software version 2.2.10, and the k-mer distribution was plotted with GenomeScope version 2.0 (k = 21; http://genomescope.org/genomescope2.0/; latest accessed: September 18, 2024) (Vurture et al., 2017).

Genome assembly evaluation

After assembly, we evaluated the quality of the final genome assembly by aligning stLFR reads using BWA version 0.7.17 (Li & Durbin, 2009). We also used Benchmarking Universal Single-Copy Orthologs (BUSCO) version 5.4.4 (Manni et al., 2021) to assess the completeness of the assembly. The BUSCO pipeline was used to test for the presence and completeness of orthologs using the plant-specific database from Embryophyta OrthoDB release 10 (Kriventseva et al., 2015). Additionally, the LTR Assembly Index (LAI) was computed and assessed using LTR_retriever version 3.0.1 (Ou & Jiang, 2018).

Repetitive sequence identification

The repetitive sequences in the assembled genome were identified with de novo and homology-based methods. First, we constructed a de novo repeat library of the mungbean using the software package RepeatModeler version 2.0.1 with default parameters (https://github.com/Dfam-consortium/RepeatModeler, accessed on 30 April 2024). This package includes distinct discovery algorithms, RECON version 1.08 and RepeatScout version 1.0.5 (Bao & Eddy, 2002; Flynn et al., 2020; Price, Jones & Pevzner, 2005). The programs were employed to identify the boundaries of repetitive elements and to build consensus models of interspersed repeats. After we obtained the repeat library, we aligned the repeat sequences to NCBI GenBank’s non-redundant protein database using BLASTX with an e-value threshold of 10⁻⁶ to verify that the library did not contain sequences belonging to large families of protein-coding sequences. Then, using the RepBase plant repeat database (Jurka et al., 2005) (https://www.girinst.org/), homology-based repetitive sequences of the mungbean genome were identified with RepeatMasker version 4.0.9_p2 with default parameters (Tempel, 2012). Additionally, we employed Tandem Repeats Finder (TRF) version 4.09 (Benson, 1999) with default parameters to predict tandem repeats element. HelitronScanner version 1.0 (Xiong et al., 2014) and miteFinder version 3.0 (Hu, Zheng & Shang, 2018) were used to identify helitrons and miniature inverted-repeat transposable elements (MITEs), respectively.

Gene annotation

Regarding gene annotation, we employed the BRAKER2 version 2.1.6 (Brůna et al., 2021) which performs gene prediction by using GeneMark-ES/ET/EP/ETP version 4.65 and AUGUSTUS version 3.4.0. BRAKER2 version 2.1.6 also enables the integration of transcript and protein homology information into the predictions. Six sets of RNA-Seq data from various tissues, SRS652317 (roots, stems, leaves), SRS1308070 (seed), SRR27190622 (petal blooming flower), SRS3960651 (V6-stage leaves), SRS9581523 (Thai SUT1 4-days-old root) and SRS702243 (5-day seedling hypocotyl) were retrieved from SRA database of NCBI. The RNA-Seq data were mapped onto the V. radiata assembly of this study with HiSAT2 version 2.2.1 (Kim et al., 2019) and used for BRAKER2 analysis. Protein data from Viridiplantae OrthoDB11 (Kuznetsov et al., 2023) and V. radiata cultivar VC1973A (GCF_000741045.1) were included in the analysis. The predicted genes were filtered and the longest isoforms of overlapping transcripts were selected using TSEBRA version 1.0.3 (Gabriel et al., 2021) with the default evidence-weight setting. BUSCO version 5.4.4 (Manni et al., 2021) was used with eudicots_odb10 to evaluate the quality of the annotation. Functional annotation of the predicted genes was performed using OmicsBox version 1.3.11 (https://www.biobam.com/download-omicsbox/, last accessed on 10 May 2024). We aligned the protein sequences against the UniProtKB/Swiss-Prot (Bairoch & Apweiler, 2000) and GenBank non-redundant (NR) databases via BLASTP (Gish & States, 1993) with an e-value cutoff of 10⁻⁵. Gene ontology (GO) terms were retrieved and assigned to V. radiata query sequences, while enzyme codes (EC) corresponding to V. radiata gene ontology were extracted and mapped to Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway annotations. We also predicted tRNAs using the tRNAscan-SE tool version 2.0 (Chan Patricia et al., 2021) (http://lowelab.ucsc.edu/tRNAscan-SE/). As rRNAs are highly conserved, we used rRNA sequences from closely related species as references and used BLAST to predict rRNA sequences. Other ncRNAs, including miRNAs and snRNAs, were predicted by searching against the Rfam version 14.1 database (Griffiths-Jones et al., 2005) (http://infernal.janelia.org/) using the Infernal software version 1.1.5 with the default parameters. The genome annotation of mungbean was deposited in the Figshare database (https://figshare.com/articles/dataset/Vigna_radiata_annotation_files/27094231).

Comparative genomics and phylogenetic analyses

OrthoFinder version 2.4.0 (Emms & Kelly, 2019) was used to identify orthologous gene groups among the three mungbean cultivars (KUML4, Vrad_JL7, and VC1973A v2) and 14 other plant species (Table S1). These included nine legumes [black gram (Vigna mungo (L.) Hepper), créole bean (Vigna reflexo-pilosa Hayata), cowpea (Vigna unguiculata (L.) Walp.), adzuki bean (Vigna angularis (Ohwi) Ohwi and Ohashi), common bean (Phaseolus vulgaris L.), barrel medic (Medicago truncatula Gaertn.), soybean (Glycine max (L.) Merr.), chickpea (Cicer arietinum) and peanut (Arachis duranensis)], three cucurbit species (cucumber (Cucumis sativus L.), melon (Cucumis melo L.) and watermelon (Citrullus lanatus L.)), one rosid species (Arabidopsis), and one monocot species (rice (Oryza sativa L.)). The three mungbean cultivars were further clustered and subjected to Venn diagram analysis to explore the common and cultivar-specific genes across the three mungbean genomes. Single-copy orthologous protein sequences were aligned with MUSCLE version 3.8.1551 (Edgar, 2004), and alignment gaps were removed using trimAl version 1.4 rev15 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009). Alignment blocks were concatenated using catsequences program (https://github.com/ChrisCreevey/catsequences, last accessed on 22 May 2024), and the best substitution model for each block was estimated using ModelTest-NG program version 0.1.7 (Darriba et al., 2019). A maximum-likelihood phylogenetic tree was constructed using RAxML-NG program version 1.0.2 (Kozlov et al., 2019) with the best-fit substitution model, and 1,000 bootstrap replicates. O. sativa was designated as the outgroup of the phylogenetic tree. Divergence times of species within the phylogenetic tree were estimated with the MCMCTree program version 4.0 in the Phylogenetic Analysis by Maximum Likelihood software (PAML4 package) (Yang, 2007) using the relaxed-clock model with the known divergence time between cucumber and melon (8.4–11.8 million years ago; Sebastian et al., 2010). Finally, we employed CAFE software version 4.2 (Han et al., 2013) to predict the significant expansion and contraction of gene families across the phylogeny (p-value < 0.05) with the gene birth-death (λ) parameters estimated using the maximum-likelihood method. The expanded and contracted genes were then subjected to GO and KEGG functional annotation. Genome collinearity analysis was performed using the MCScan tool (Python version, https://github.com/tanghaibao/jcvi/) from jcvi v1.3.885 (Tang et al., 2024; Wang et al., 2012). MCscan was employed to define syntenic blocks between the V. radiata KUML4, V. radiata JL7, and V. mungo genomes. Collinearity between the three genomes was visualized using collinearity plots.

Genome assembly and evaluation

To construct the chromosome-scale assembly of the V. radiata genome, we integrated datasets from stLFR reads and Hi-C sequences (Table S3). This initial genome assembly was generated from a total of 501,995,196 reads totaling 121.48 Gb of stLFR sequencing data (225.12 × coverage based on the estimated genome size obtained from k-mer analysis), resulting in a draft genome of 464.02 Mb with 6,685 scaffolds. The longest scaffold was 17.68 Mb, and the N50 length was 4.64 Mb. Additionally, we used 50.32 Gb of Hi-C data to generate chromosome-level scaffolds. The preliminary assembly was scaffolded using the long-range Hi-C technique, resulting in a final assembly containing 468.08 Mb in 5,834 scaffolds with an N50 length of 40.75 Mb (Table S2 and Fig. S1). We estimated the size of the V. radiata KUML4 genome based on the k-mer frequency histogram derived from short-read stLFR sequences (Fig. S2). The assembly size (468.08 Mb; Fig. 1) corresponded well with the analysis of k-mer distribution of PE reads (445.98 Mb). Our assembly size was about 2.35% smaller than the estimated size of the Vrad_JL7 genome (479.35 Mb; Table 1) (Liu et al., 2022).

The final assembly was anchored into 11 pseudochromosomes larger than 25 Mb, corresponding to the haploid chromosome number in V. radiata (2n = 2x = 22). These pseudomolecules covered 432.35 Mb, or 96.94% of the estimated genome size and were named in accordance with Liu et al. (2022). To evaluate the quality and accuracy of the V. radiata KUML4 assembly, we first mapped the stLFR reads back to the genome assembly and obtained a mapping rate of 98.09%. Second, we evaluated the assembly using 1,614 Benchmarking Universal Single Copy Orthologs (BUSCO) genes from Embryophyta to assess the completeness of the conserved orthologs in the V. radiata genome. Our assembly demonstrated high levels of completeness, with proportions of complete (C), complete and single-copy (S), complete and duplicated (D), fragmented (F), and missing (M) genes as follows: C: 98.3% [S: 95.0%, D: 3.1%], F: 0.8%, M: 1.0% (Table S2). Finally, the long terminal repeat (LTR) assembly index (LAI) was employed to assess the completeness of the assembled LTR as an indicator of the quality of the whole genome assembly. The LAI score of the assembly was 11.57 (Table 1), meeting the reference genome level according to developer classification standards (Ou, Chen & Jiang, 2018).

Gene annotation

A combination of several types of evidence, including ab initio prediction, RNA-seq supported evidence, and protein homology evidence, was integrated into the genome annotation pipeline to obtain reliable gene predictions. The genome annotation contained 28,593 predicted gene models, of which 27,667 (96.76%) were protein-coding genes. The BUSCO completeness for the predicted genes was 94.1%. The average gene length was 2,793 bp with 4.65 exons per gene, and the average exon and intron length were 236 bp and 464.93 bp, respectively (Table 2). We observed that average GC content of the V. radiata genome assembly was 33.27% (Table S2), which was close to the average GC content in introns (32.16%), whereas the average GC content in exons was higher at 44.50% (Table 2). A total of 18,334 predicted genes had been assigned gene ontology (GO) terms. The top three most prevalent GO terms associated with biological process were protein phosphorylation, regulation of DNA-templated transcription, and transmembrane transport, while the most prevalent terms associated with the cellular component were membrane followed by nucleus, and cytoplasm. The largest categories of genes assigned to molecular function were ATP binding followed by metal ion binding, and DNA binding (Fig. S3). Of 28,593 predicted gene models, 96.76%, 70.52%, and 64.12%, could be identified functionally using the NCBI nonredundant (NR) protein, SwissProt and GO databases, respectively (Table S4). For the non-coding RNAs, we also identified 28,017 microRNAs, 765 transfer RNAs, 483 ribosomal RNAs and 12,497 small nuclear RNAs in the V. radiata genome (Table S5).

Identification of repetitive elements

The V. radiata KUML4 genome assembly contained a total of 253.85 Mb of repetitive elements, accounting 54.76% of the assembly (Table 3). These repeat elements comprised 25.72 Mb (10.12%) of DNA transposons, 8.04 Mb (3.17%) of simple sequence repeats, 119.09 Mb (46.92%) of retrotransposons, 35.49 Mb (13.98%) of tandem repeats, 18.53 Mb (7.30%) of miniature inverted-repeat transposable elements (MITEs), and 1.60 Mb (0.63%) of helitron elements. Retrotransposons constituted the majority of known repeats, comprising nearly half of the total repeat content in the V. radiata KUML4 genome. Within the retrotransposon class, long terminal repeat (LTR) retrotransposons were the most frequently detected, accounting for 46.23% of the repetitive sequences. The most abundant LTR superfamilies, Gypsy and Copia, accounted for 15.65% and 9.65% of the mungbean genome, respectively (Table 3), which differed from the 16.65% and 13.77% in the Vrad_JL7 genome (Liu et al., 2022).

Comparative genomic analysis

To explore the evolutionary relationships among the three mungbean cultivars (KUML4, Vrad_JL7, and VC1973A v2) and other plant species, the gene sets from the five Vigna species (V. radiata, V. mungo, V. reflexo-pilosa, V. unguiculata, and V. angularis), five other legumes (P. vulgaris, M. truncatula, G. max, C. arietinum, and A. duranensis), four other eudicots (A. thaliana, C. sativus, C. melo, C. lanatus), as well as one monocot (O. sativa) were used for gene family clustering. Of the total 670,738 proteins in 17 genomes from 15 species, 641,506 proteins (95.6%) were clustered into 36,504 orthologous groups. Among the 27,625 protein-coding genes in KUML4 that had orthologs in the other species analyzed, only 256 (0.93%) were specific to KUML4, while 575 (2.08%) were shared among legume species (Fig. 2A). Comparison of gene families among the three mungbean genomes revealed that 238, 1,150, and 308 orthologous groups were unique to KUML4, Vrad_JL7 and VC1973A v2, respectively, while 18,152 orthologous groups were shared among KUML4 and the two previously published mungbean genomes (Vrad_JL7 and VC1973A v2; Fig. 2B). Sequence information from single-copy orthologous genes was used to construct a maximum-likelihood tree with O. sativa as an outgroup, and the divergence time was estimated based on the topology and branch length. The phylogenetic analysis revealed that the three V. radiata cultivars (KUML4, Vrad_JL7 and VC1973A v2) and V. mungo were closest and diverged approximately 4.17 million years ago (MYA). The divergence of KUML4 from Vrad_JL7 occurred around 1.38 MYA. The ancestor of V. radiata and V. mungo formed a sister clade to the ancestor of V. angularis and V. reflexo-pilosa. The estimated divergence time between the two clades was about 7 MYA (Fig. 2C). Additionally, our analysis of gene family expansion and contraction across the three mungbean cultivars and other plant species identified 2,764 significant gene families. A total of 39 significantly expanded gene families, and 932 significantly contracted families were identified in KUML4 (Fig. 2C).

Functional categories indicated that a large number of expanded gene families were associated with secondary metabolite biosynthesis, including berberine bridge enzyme-like 13, isoflavone 2′-hydroxylase, UDP-glycosyltransferase 73D1. In addition, they were enriched in the defense responses, and signal transduction related to plant hormones like ethylene and cytokinin, including aspartic proteinase CDR1-like, GDSL esterase/lipase EXL3, ethylene-responsive transcription factor ERF098-like, zeatin O-glucosyltransferase-like, and cytokinin riboside 5′-monophosphate phosphoribohydrolase (Table S6). The contracted gene families were mainly involved in the cysteine-rich receptor-like protein kinase 10, proline-rich receptor-like protein kinase PERK4, shaggy-related protein kinase epsilon isoform X1, and MDIS1-interacting receptor-like kinase (Table S6).

We also explored genome synteny among V. radiata KUML4, V. radiata JL7, and V. mungo. Our genome assembly of KUML4 reveals a high level of synteny with the Vrad_JL7 genome, as illustrated by the large syntenic blocks in Fig. 3. Although some structural variations were detected, including inverted regions, these were observed on chromosomes 1, 2, 3, 6, 9, 10, and 11 between the KUML4 and Vrad_JL7 genomes.