Vicia sativa L. (common vetch, NCBI:txid3908) (Figure 1) is an annual legume belonging to the Fabaceae family, and Vicia genus [1]. The Vicia genus contains about 180–210 species, including the economically important crop broad bean [2]. To date, no chromosome-level genome assembly has been reported within the Vicia genus. Interestingly, V. sativa has at least three different reported haploid chromosome numbers: n = 5, 6 or 7 [3], but n = 6 is the best characterized karyotype.

V. sativa is thought to have originated in the Fertile Crescent of the Middle East and is now widespread on every continent as both a crop and a weed [4]. V. sativa is a multipurpose legume; the plants are often grown for forage and the seeds can be used safely as a feed for ruminant animals. V. sativa seed contains up to 30% crude protein and is rich in essential amino acids and unsaturated fatty acids [5]. However, only a small amount of the seed can be safely fed to monogastric animals like chickens and pigs, because of the presence of the neurotoxic proteinaceous amino acids 𝛽-cyano-L-alanine and 𝛾-glutamyl-𝛽-cyano-alanine [6].

V. sativa is often used in crop rotation systems to increase nitrogen input to the soil. In a study of V. sativa/wheat rotation over a 4-year-period, cultivation of V. sativa during autumn increased soil water storage and subsequently increased biological yield and grain yield of wheat. Both yields were doubled in the third year compared with the second year of the rotation [7]. Furthermore, the symbiosis between soil rhizobia bacteria and V. sativa roots allows the plant to fix atmospheric nitrogen and later provide nitrogen for the following crop, hence reducing the use of expensive nitrogen fertilizer [8]. V. sativa exhibits excellent drought tolerance and is suitable for cultivation in arid areas. In one drought tolerance study, V. sativa could withstand a month of drought stress, with the leaf weight not decreasing significantly compared with the non-drought control [9]. V. sativa offers multiple usage and is a valuable crop in a sustainable agricultural system [10].

With the important value of V. sativa, vetch breeders have primarily selected for traits conferring high yield, pod shattering, flowering time, disease resistance against Ascochyta fabae, Uromyces viciae-fabae (rust) and Sclerotinia sclerotium [11]. Recently published transcriptome data has allowed agriculturally important traits to be uncovered at the gene expression level, such as pod-shattering resistance [12] and drought tolerance genes [13] in V. sativa. However, a lack of high-quality genome reference is currently impeding the genetic mapping of important genes and hindering further applications such as genome editing when compared with other crops.

In this study, we assembled a high-quality chromosome-level reference genome for V. sativa, which is the first chromosome-level reference genome in the Vicia genus. We performed genome annotation using RNA-seq data from five tissues to ensure most of the expressed genes were captured. We also included a phylogenetic analysis of V. sativa and legume relatives. We envisage that our V. sativa genome will be an important resource for evolutionary studies of this species. The well-annotated chromosome-level genome will also provide important information to facilitate genetic mapping, gene discovery and functional gene studies.

To prepare V. sativa for whole genome sequencing (WGS) using long-read and short-read data, seeds of cultivar Studenica (V. sativa subsp. sativa) were obtained from the South Australian Research and Development Institute (SARDI, South Australia, Australia). Seeds were sterilized and germinated in vitro on half-strength Murashige & Skoog (1/2 MS) basal medium with 1% sucrose for 3 days at 25 °C, in the dark. Bulk 3-mm-long primary root tips were then harvested and snap-frozen in liquid nitrogen for subsequent DNA extraction. DNA was extracted using the phenol:chloroform method [14], with an additional high-salt low-ethanol wash to improve DNA purity [15]. High-quality DNA was confirmed by electrophoresis on 1% agarose gel. The DNA was sent to the Australian Genome Research Facility (AGRF, Melbourne, Australia), and Novogene Co., Ltd (Hong Kong, China) for library preparation and sequencing on a PromethION (PromethION, RRID:SCR_017987) and Novo-Seq 6000 (Illumina NovaSeq 6000 Sequencing System, RRID:SCR_016387), respectively. We obtained 72 gigabase pairs (Gbp) of Nanopore long-read data, and 205 Gb paired-end short-read data (150 base pairs [bp] read length).

To produce V. sativa CHiCAGO sequencing data [16] and Hi-C sequencing data [17], 2 g of young leaf tissue was snap-frozen in liquid nitrogen and sent to Dovetail Genomics (USA) for library preparation and sequencing. CHiCAGO and Hi-C libraries were sequenced on an Illumina HiSeq X (Illumina HiSeq X Ten, RRID:SCR_016385) to produce 162 Gbp of CHiCAGO and 148 Gbp of Hi-C sequencing data, respectively.

To prepare V. sativa RNA sequencing (RNA-seq) data, RNA was purified from the first two fully expanded leaves, shoot apexes with young leaves up to 1 cm long from 4-week-old plants, roots from 5-day-old seedlings and 4-week-old leaf-derived callus tissues using the Spectrum™ Plant Total RNA Kit (Sigma Aldrich). Additional DNase I treatment was used to remove DNA contamination (On-Column DNase I Digestion, Sigma Aldrich), and ribosome removal treatment to enrich for the non-ribosomal RNA fraction (Ribo-Zero rRNA Removal Kit for Plant Leaf or Plant Seed/Root, Illumina) [18]. Directional RNA libraries were prepared for each tissue using the NEBNext Ultra™ Directional RNA Library Prep Kit for Illumina (New England Biolabs) following the manufacturer’s protocol. Libraries were sent to Novogene Co., Ltd (Hong Kong, China) for sequencing on Novo-Seq 6000 (Illumina) to obtain 150-bp paired-end read data. In total, we obtained 74.6 Gbp of RNA-seq data. A summary of the long and short-read sequencing data is provided in Table 1.

We first performed a genome size estimation for V. sativa. To do this, short-Illumina (paired-end 150 nt) reads were trimmed using TrimGalore v0.4.2 (Trim Galore, RRID:SCR_011847) with default parameters and 25-mers were counted using Jellyfish v2.2.6 (Jellyfish, RRID:SCR_005491) [19]. The 25-mer count distribution data was used as an input to GenomeScope (GenomeScope, RRID:SCR_017014) [20] for genome size estimation with the read length set to 150 and max k-mer coverage set to 1 million. GenomeScope estimated a genome-wide heterozygosity level of 0.09% (Figure 2) and a genome size of 1.61 Gbp; approximately 160 megabase pairs (Mbp) smaller than the genome size estimated by flow cytometry (1.77 Gbp) [21].

Next, we conducted contig assembly from the Nanopore long-reads using Canu v1.7 (Canu, RRID:SCR_015880) [22] under default parameters with the expected genome size set at 1.77 Gbp. Canu was used to perform read trimming and sequencing error correction for the input data before performing contig assembly. The assembled contigs were polished using clean WGS short-reads with Pilon v1.22 (Pilon, RRID:SCR_014731) [23] under default parameters. We repeated the polishing step and observed a further improvement in contig quality (Table 2). Contig quality was assessed using BUSCO v5.2.2 (BUSCO, RRID:SCR_015008) [24] for the completeness of the genome, and after two rounds of polishing, complete BUSCOs increased from 69.9% to 97.8% (Table 2). Overall, we obtained 9,990 assembled contigs, which were 1.93 Gbp, with an N50 value of 685 kilobase pairs (Kbp).

To generate a chromosome-level assembly for V. sativa, Hi-C proximity [25] ligation data and CHiCAGO [26] were used to order and orient the contigs along chromosomes. The scaffolding process was carried out by Dovetail Genomics (Santa Cruz, CA, USA) using Dovetail™ Hi-C library reads to connect and order the input set of contigs. After scaffolding with HiRise (v2.1.7) [51], the assembled genome sequence initially comprised 1.8 Gbp, with a scaffold and contig N50 of 51.4 and 0.6 Mbp, respectively. A high fraction of the assembled sequences (93%) were contained in four pseudo-chromosomes (Figure 3A); however V. sativa has six pairs of chromosomes [1]. We observed that two of the four pseudo-chromosomes had weak interactions, suggesting misjoining of two contigs (Figure 3A).

In parallel to the HiRise analysis, we performed a second chromosome-level assembly using 3D-DNA (3D de novo assembly, RRID:SCR_017227) [27]. 3D-DNA scaffolding was performed by first mapping Hi-C reads to the contig assembly using Juicer v1.6 (Juicer, RRID:SCR_017226) [28], which then generated 304,484,352 uniquely mapped pair-end reads and of which 51.1% (155,477,211) of the uniquely mapped reads were identified as valid Hi-C contacts. The 3D-DNA v180114 pipeline was integrated to anchor contigs to pseudo-chromosomes based on valid Hi-C contacts. The output file was used to generate a Hi-C heatmap for manual inspection using Juicebox Assembly Tools v1.11.08 (Juicebox, RRID:SCR_021172). This revealed six high-quality pseudo-chromosomes (Figure 3A).

We compared the HiRise and 3D-DNA assembled pseudo-chromosomes by performing a whole genome alignment with Mummer v4.0.0 (MUMmer, RRID:SCR_018171) [29]. The alignment showed a strong synteny between the HiRise and 3D-DNA pseudo-chromosomes (Figure 3B). However, the two longest HiRise pseudo-chromosomes aligned to four 3D-DNA pseudo-chromosomes suggesting two HiRise pseudo-chromosomes were misjoined (Figure 3B). The putative misjoined HiRise pseudo-chromosomes also coincided with low Hi-C interactions (Figure 3A).

To further support that these two HiRise pseudo-chromosomes were misjoined, we compared the synteny of the HiRise and 3D-DNA pseudo-chromosomes to the high-quality V. faba genetic linkage map [30] as no genetic linkage map is available for V. sativa. When we compared the order of 1536 sequenced V. faba DNA markers to their homologous regions in our HiRise and 3D-DNA pseudo-chromosomes, we observed a clear synteny between V. faba and V. sativa. However, two out of four of the HiRise pseudo-chromosomes appeared to be misjoined, for example, the markers on HiRise pseudo-chromosomes one, mapped to two V. faba linkage groups (Figure 4). After combining the karyotype, Hi-C interaction and synteny data to V. faba, we concluded the 3D-DNA assembly was most likely correct and subsequently used this assembly in further analysis. Finally, we used purge_dup pipeline v1.2.5 (purge dups, RRID:SCR_021173) [31] to remove low coverage scaffolds, partial overlaps and haplotigs. The final version of the genome assembly contains six pseudo-chromosomes (Table 3), in which a total of 1.65 Gbp contigs are anchored to these pseudo-chromosomes (Figure 5), and remain 10 unassigned contigs (overall 334,511 bp length). The overall genome size is 1,653,553,227 bp, with a GC content of 35.6% (Table 4).

To annotate the V. sativa genome assembly, we masked repeat regions of the genome, then mapped the RNA-seq data to the masked genome and performed gene prediction. First, the repeat families found in the V. sativa genome were identified de novo and classified using the software package EDTA v1.9.6 [35] with the sensitive model setting. EDTA integrated multiple programs, including LTR_FINDER (LTR_Finder, RRID:SCR_015247) [36] and RepeatModeler (RepeatModeler, RRID:SCR_015027), which generated a non-redundant transposable element (TE) library used to annotate the TE regions on the genome. The TE library generated from EDTA was also used as an input to RepeatMasker v4.1.2 (RepeatMasker, RRID:SCR_012954) to identify and perform “hard-masking” and “soft-masking” for the repetitive region on the genome. A total of 83.9% of the genome was masked, and 64.4% of the genome was detected as LTR elements (Table 6).

After genome masking, a combination of ab initio prediction and transcript evidence from the RNA-seq was used for gene prediction. Briefly, each RNA-seq data set was trimmed for low quality bases using TrimGalore v0.4.2, and mapped to the hard-masked-genome by using STAR v2.7.9 (STAR, RRID:SCR_004463) [37] to generate BAM files. Then the soft-masked genome and the BAM files generated from STAR were used for gene prediction using BRAKER v2.1.6 (BRAKER, RRID:SCR_018964) [38]. A total of 53,218 predicted protein-coding-genes were identified (Table 7). To assess the completeness of these protein-coding-genes, BUSCO v5.1.3 with fabales_odb10 gene sets were used which then identified 5127 (95.6%) complete, 395 (7.4%) duplicated, 70 fragmented (1.3%) and 169 missing (3.1%) orthologs.

Putative functions of the predicted protein-coding-genes were characterized by comparing the predicted proteins against the SwissProt and National Center for Biotechnology Information (NCBI) non-redundant database using Diamond v2.0.11 (DIAMOND, RRID:SCR_016071) [39] with e-value cut-off of 1 × 10−5. Protein motifs and domains were annotated by comparing the predicted proteins against the InterPro database using Interproscan v5.52-86.0 (InterProScan, RRID:SCR_005829) [40]. The predicted proteins were also assigned with Gene Ontology (GO) terms corresponding to the InterPro entries using Interproscan v5.52-86.0. In addition, we compared the predicted proteins against the EggNOG database v5.0 (eggnog, RRID:SCR_002456) [41] using eggNOG-mapper v2.1.4-2 (eggNOG-mapper, RRID:SCR_021165) [42] and assigned them with Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and KEGG orthologous groups (KO). As a result, we were able to annotate 47,580 (89.4%) predicted protein-coding genes with at least one function term (Table 8).

In addition, we also identified and annotated non-coding RNA in the V. sativa genome. tRNA was identified using tRNAscan-SE v2.07 (tRNAscan-SE, RRID:SCR_010835) [43], rRNA was identified using Rnammer v1.2 (RNAmmer, RRID:SCR_017075) [44] and other types of non-coding RNA were identified by using Infernal v1.1.4 (Infernal, RRID:SCR_011809) [45] based on the Rfam database (Rfam, RRID:SCR_007891) [46]. Overall, 3966 of noncoding genes were annotated, including 158 miRNA, 649 rRNA and 1777 snRNA (Table 9).

We identified the orthogroups, phylogenetic positions and divergence times between V. sativa and 11 other plant species. The source of the protein-coding sequences used in our analysis are listed in Table 10. First, protein sequences of V. sativa, Pisum sativum, Medicago truncatula, Trifolium pratense, Phaseolus vulgaris, P. lunatus, Vigna unguiculata, Chamaecrista fasciculata, Faidherbia albida, Cercis canadensis, Carya illinoinensis, and Arabidopsis thaliana [47–54] were clustered into orthogroups using Orthofinder v2.5.4 (OrthoFinder, RRID:SCR_017118) [55] with default parameters. A total of 10,009 single-copy and 43,209 multi-copy genes were identified in the V. sativa annotation (Figure 6B, Table 11), forming 19,096 orthogroups (Figure 6A, Table 10). Comparing orthogroups amongst V. sativa, P. sativum, M. truncatula, P. vulgaris, F. albida, we identified 2309 orthogroups that are specific to V. sativa (Figure 6A). Orthofinder was further used to perform phylogenetic reconstruction with the multiple sequence alignment approach (using the -msa command). The generated species tree has a support value of one on all nodes (Figure 7), indicating the high reliability of the revealed phylogenetic relationships.

To estimate divergence times between V. sativa and other important legume species (Table 10), coding sequences of 64 randomly selected single copy orthologs (see Supplementary File 1 [57]) were aligned using MACSE v1.2 [58]. Low-quality regions of each alignment were trimmed using Trimal v1.4.1 (trimAl, RRID:SCR_017334) [59], resulting in high-quality alignments amounting to 139,956 bp. Individual alignments were then imported into Beast v2.6.3 (BEAST2, RRID:SCR_017307) [60] for phylogenetic dating. Substitution models were selected using BEAST Model Test [61] for each alignment and were allowed to coalesce using unlinked relaxed log-normal molecular clocks [62]. A calibrated Yule prior [63] was used to inform tree building and speciation with four node calibrations (Table 12). First, a log normal distribution of 89.3 MYA (5% quantile 97.9 million years ago [MYA], median 106 MYA, 95% quantile 121 MYA) [64] was used to inform the root prior (Brassicaceae, Fabaceae split). Three fossil calibrations were then set using CladeAge [65]: (i) Fabaceae (Figure 8 red dot; 65.3 MYA) [66], (ii) Caesalpinioideae (Figure 8 blue dot; 58 MYA) [67, 68], and (iii) Papilionoideae (Figure 8 green dot; 55 MYA) [69]. Furthermore, a net diversification rate was set to 0.1–0.134 to construct a distribution around the literature value of 0.117 [70], turnover rate was set to 0.823–0.883 to construct a distribution around the literature value of 0.853 [70], and sampling rate of 0.000034–0.013 [71] was set to determine CladeAge prior distributions. Final chain length of the Markov Chain Monte Carlo (MCMC; 600 million) was determined through continuous examination of the log file using Tracer until proper mixing was observed. This allowed us to determine a robust estimate for the most common recent ancestor (MRCA) of V. sativa and P. sativum at 10.6 (95% Highest Posterior Density: 9.9–11.4) MYA (Figure 8). Gene family expansion and contraction analysis using CAFE v4.2.1 (Computational Analysis of gene Family Evolution, RRID:SCR_018924) [72] with a single 𝜆 revealed 5195 gene families that have undergone gene expansion (3024) or contraction (2171) since the MRCA of V. sativa and A. thaliana (Figure 6C).

To identify whole genome duplication events (WGD), WGDI v0.5.1 [73] was used to identify gene collinearity between V. sativa, M. truncatula and P. vulgaris. The Ks (synonymous substitutions per synonymous site) value was calculated based on the identified collinearity gene to construct a frequency distribution map. The Ks distribution indicated that V. sativa, M. truncatula and P. vulgaris share the same ancestral WGD event. The estimated time of this WGD event (∼58 MYA) [74] and the corresponding Ks value (∼0.93, Figure 6D) reveal that the average mutation rate of V. sativa genome is 8.02 ×10−9 per site per year.