Chromosome-level genome assembly of Helwingia omeiensis: the first genome in the family Helwingiaceae

TL;DR

This paper presents the first high-quality genome assembly for the plant family Helwingiaceae, offering insights into its unique traits and evolutionary history.

Contribution

The first chromosome-level genome assembly for the family Helwingiaceae, enabling molecular studies of its unique floral traits.

Findings

The genome of Helwingia omeiensis is assembled at 2.75 Gb with 98.2% BUSCO completeness.

53,951 genes were identified and 99.7% were annotated in protein databases.

The genome provides a foundation for studying the evolution of Helwingiaceae-specific traits.

Abstract

Helwingia, a shrub of the monotypic cosmopolitan family Helwingiaceae, is distinguished by its inflorescence, in which flowers are borne on the midrib of the leaf—a trait not commonly observed in related plant families. Previous studies have investigated the development of this unusual structure using comparative anatomical methods. However, the scarcity of genomic data has hindered our understanding of the origins and evolutionary history of this uncommon trait at the molecular level. Here, we report the first high-quality genome of the family Helwingiaceae. Assembled using HiFi sequencing and Hi-C technologies, the genome of H. omeiensis is anchored to 19 chromosomes, with a total length of 2.75 Gb and a contig N50 length of 6.78 Mb. The BUSCO completeness score of the assembled genome was 98.2%. 53,951 genes were identified, of which 99.7% were annotated in at least one protein database. The high-quality reference genome of H. omeiensis provides an essential genetic resource and sheds light on the phylogeny and evolution of specific traits in the family Helwingiaceae.

Figures

Distribution profiles of 27-mer analysis of short reads.Fig. 2The Hi-C interactive heatmap of 19 pseudo-chromosomes of *H. omeiensis*.Table 2Statistics of chromoslme-level genome assembly of *H. omeiensis*.FeaturesStatisticsAssembly size (bp)2,748,316,618Largest sequence length(bp)153,486,838Counts of scaffold sequences1,584Scaffold N50 (bp)127,781,096Counts of scaffold N5010Scaffold N90 (bp)16,574,032scaffold L9023Contig N50 (bp)6,783,361Contig L50119Contig N90 (bp)927,000Counts of contig N90502GC content(%)35.19N Length193,200N content (%)0.007BUSCO completeness (%)98.2Fig. 3Circos plot of g

Venn diagram displaying the matches of genes of *H. omeiensis* in five public protein databases.

Full text

Background & Summary

Helwingiaceae is a monotypic family in the order Aquifoliales, comprising a single genus Helwingia. The innovative structure of this genus is that the flowers are borne on the midrib of the leaf, which is known as an “epiphyllous inflorescence”, setting them apart from other plants. In addition, the pith, leaves, and fruits of plants in this genus are traditionally used in herbal medicine to treat dysentery and as diuretic and anti-inflammatory remedies^1^. The genus includes four species, H. chinensis, H. himalaica, H. japonica, and H. omeiensis, which are all dioecious shrubs mainly found in eastern Asia^2,3^. Specifically, H. omeiensis is indigenous to Southwest China, and thrives in moist woodlands and on mountain slopes^2^.

Previous comparative anatomical studies suggested that changes in the position of flower primordium initiation and intercalary growth may contribute to the formation of this distinct structure^4–6^. With the development of high-throughput sequencing technologies, the genomes of three closely related species in the genus Ilex of the family Aquifoliaceae have been published^7,8^. However, despite the fact that RNA-seq data and the complete chloroplast genomes of three Helwingia species have been released^4,9,10^, a lack of genomic data remains a barrier to studying the evolutionary origin of the family.

In this study, we leveraged a combination of short reads, high-fidelity (HiFi) reads, and chromosome conformation capture (Hi-C) sequencing data to construct a chromosome-level genome assembly for H. omeiensis, providing the first genome resource for the family Helwingiaceae. The length of the genome assembly was 2.75 Gb, with a scaffold N50 of 127.8 Mb and a contig N50 of 6.78 Mb. We identified 1.98 Gb of repetitive elements, accounting for 72.21% of the assembled genome, as well as 53,951 protein-coding genes. The genome assembly and annotation of H. omeiensis will provide a critical foundation for exploring the genetic basis underpinning of this unique inflorescence structure and the phylogenetic relationships within the family Helwingiaceae.

Methods

Plant materials

All of the fresh materials were collected from a female adult plant of Helwingia omeiensis cultivated in Mount Emei Botanical Garden, Sichuan Province, China (N29°35′40, E103°22′40), and the specimens were kept at the Museum of Sichuan University. The genomic DNA was extracted from young leaves, whereas RNA was extracted from mature leaves and terminal buds.

Library construction and sequencing

For short-read sequencing, the sample was randomly fragmented by an ultrasonic processor (Covaris S220; Woburn, MA, USA) to generate DNA fragments approximately 350 bp in length. The DNA fragments were subsequently constructed through end repair, the addition of a 3′ A tail and the ligation of adapters. Next, the library was sequenced with a DNBSEQ-G400 (BGI, Wuhan, China). The raw short reads were filtered by SOAPnuke v1.5.6^10^ to remove adapters and low-quality reads. A total of 87.36 Gb of clean data were obtained for H. omeiensis (Table 1).Table 1. Statistics of the sequencing data of the H.omeiensis genome.Library typesMoleculePlatformInsert sizeData size (Gb)Short-readDNADNBseq300–400 bp87.36HiFiDNAPacBio Sequel II20 Kb50.32Hi-CDNAIllumina HiSeq X Ten—221.52RNA-SeqRNAIllumina HiSeq X Ten—10.46

For HiFi (high-fidelity) sequencing, high-quality genomic DNA was sheared using Megaruptor^®^ 3 (Diagenode), and subreads with a length of 20 kb were further selected using Sage ELF to prepare the PacBio HiFi libraries in CCS mode on the Pacific Biosciences Sequel II System (Supplementary Figure S1). Finally, 50.32 Gb of long clean reads were generated (Table 1), with mean lengths of 13.0 kb and 14.5 kb, respectively.

Hi-C technology captures sequence interactions between all DNA segments within chromosomes to obtain information on interactions between segments of the genome for assisted genome assembly^11^. Fresh leaves of the same individual were used to construct Hi-C libraries, and the MboI restriction enzyme was used for DNA ligation. After tailing, pulldown, and adapter ligation, the DNA library was sequenced on an Illumina HiSeq X Ten System (BGI, Wuhan, China) with a strategy of 2 × 150 bp. After filtering low-quality reads, 221.52 Gb of clean Hi-C data were obtained (Table 1).

RNA sequencing

Mature leaves and young terminal buds of the same individual were collected for RNA extraction. The RNA-seq library was constructed using the Illumina standard protocol (San Diego, CA, United States) and sequenced on the Illumina HiSeq X Ten platform (BGI, Wuhan, China). The raw data were filtered by Cutadapt v1.16^12^ to remove adapters and low-quality reads. After quality control by FastQC v0.11.8 (https://github.com/s-andrews/FastQC), 10.46 Gb of paired-end short clean reads were generated from the RNA-seq library (Table 1).

Genome survey and de novo assembly

Jellyfish v2.1.4^13^ was used to quickly count K-mer frequencies ranging from 17 to 31, and then GenomeScope^14^ predicted genomic features using a K-mer-based statistical approach (Supplementary Table S1). The H. omeiensis genome was estimated to be 2.54 Gb in size, with a heterozygosity rate of 1.19% and repetitive sequences accounting for 54.85% of the total length of the genome (Fig. 1). Using 50.32 Gb of clean HiFi reads with hifiasm v0.19.6-r595^15^, we generated a genome assembly of 2.92 Gb in size with a contig N50 of 6.21 Mb. Following that, Chromap v0.2.5-r473^16^ was utilized to align Hi-C clean reads to the contig assembly, and according to the strength of interactions between pairs of reciprocal sequences, YaHS v1.2a.1^17^ was used to anchor contigs onto 1,584 scaffolds. Next, using Juicebox v1.11.08^18^, we visualized the Hi-C contact maps of the scaffold assembly and made final refinements to the genome assembly. With reference to chromosome counts indexed in the Chromosome Counts Database (CCDB)^19^ (https://ccdb.tau.ac.il/) and the whole-genome Hi-C interaction heatmap, we identified the 19 longest scaffolds as pseudo-chromosomes (Fig. 2). TGS-GapCloser v1.2.1^20^ filled 75 of the 1,011 gaps in the scaffold assembly based on HiFi reads. The final assembly had a total length of 2.75 Gb, with a contig N50 of 6.78 Mb. The length of 19 pseudochromosomes was 2.38 Gb, with a maximum chromosome length of 153.79 Mb (Table 2). Since there is no reference genome for this species, we numbered the chromosomes in order from largest to smallest (Fig. 3 and Table 3).Fig. 1. Distribution profiles of 27-mer analysis of short reads.Fig. 2. The Hi-C interactive heatmap of 19 pseudo-chromosomes of H. omeiensis.Table 2. Statistics of chromoslme-level genome assembly of H. omeiensis.FeaturesStatisticsAssembly size (bp)2,748,316,618Largest sequence length(bp)153,486,838Counts of scaffold sequences1,584Scaffold N50 (bp)127,781,096Counts of scaffold N5010Scaffold N90 (bp)16,574,032scaffold L9023Contig N50 (bp)6,783,361Contig L50119Contig N90 (bp)927,000Counts of contig N90502GC content(%)35.19N Length193,200N content (%)0.007BUSCO completeness (%)98.2Fig. 3Circos plot of genomic characteristics and annotation of the H. omeiensis genome assembly (1 Mb window size). The data from the outer to the inner circles are as follows: (a) length of 19 pseudochromosomes, (b) GC density, (c) gene density, (d) Gypsy retrotransposon density, (e) Copia retrotransposon density, and (f) DNA transposon density.Table 3. Statistics of 19 chromosomes of H. omeiensis genome assembly.No.Length (bp)Number of gapsGC percentage (%)chr11534868383735.66chr21491337244534.64chr31414812712135.36chr41411236493335.04chr51407775984035.28chr61393091074835.52chr71381972172535.42chr81315425325035.34chr91305453433935.49chr101277810963735.20chr111268495563635.61chr121239010922835.34chr131218415273035.29chr141203482313234.64chr151191691702635.52chr161059073202934.74chr171029474812435.28chr18914279181835.36chr19791234872035.16

Gene annotation

To perform a comprehensive prediction of protein-coding genes, the GETA v2.5.6 pipeline (https://github.com/chenlianfu/geta) was used for automatic genome-wide annotation. First, RepeatModeler v2.0.3^21^ and DeepTE^22^ were used for self-training and to construct a repeat library. On this basis, RepeatMasker v4.1.2-p1^23^ was employed to predict and combine repetitive elements for homology-based methods. The analysis revealed that 72.21% of the genome was composed of repetitive sequences, including 46.39% long-terminal repeat (LTR) retrotransposons and 19.43% DNA transposons (Table 4).Table 4. Classification of repetitive sequences of H. omeiensis genome.TypeNumber of elementsSequence length (bp)Percentage (%)RetroelementsLTR elementsGypsy701,5311,006,106,96736.60Copia185,233172,265,8916.26Caulimovirus12,35913,943,8580.50Cassandra410573,7240.02DIRS12878,3530.00Total929,6261,275,113,39046.39LINE26,02443,464,8431.58SINE13,2252,251,4760.08Total968,8751,320,734,58048.05DNA transposons1,318,664534,237,27719.43Rolling-circles6,34612,191,6570.44Simple repeat366,17317,496,4500.63Low complexity62,0604,145,2990.15Satellite13,48911,312,6720.41snRNA1223,9950.00tRNA4340,9220.00rRNA1,2444,406,5400.16Unknown441,975132,712,7824.82Total3,178,8811,984,642,02172.21

After masking repetitive sequences in the genome, three strategies (homology-based, RNA-seq-guided, and ab initio methods) were used for the annotation process. For the RNA-seq-guided method, the RNA sequencing data were provided to HISAT2 v2.1.0^24^ and SAMtools v1.11^25^ to map the data to the repeat-masked genome. Then, TransDecoder v5.5.0 (https://github.com/TransDecoder/TransDecoder) was used to predict the open reading frame (ORF), and filter out the gene models with identities greater than 80% at the amino acid level between pairs to obtain nonredundant results. Protein sequences from Vitis vinifera, Arabidopsis thaliana, Solanum lycopersicum, Daucus carota, and Ilex latifolia were aligned to the query genome as homologous proteins using GeneWise v2.4.1^26^ to estimate protein-coding genes (Supplementary Table S2). Ab initio prediction was carried out with AUGUSTUS v3.4.0^27^, which guided by previous prediction results. Based on the GETA pipeline, all the outputs were validated using HMMER v3.3.2^28^ and NCBI-BLAST + v2.13.0 + before being integrated into a complete and nonredundant set of gene annotations.

Following the alignments by DIAMOND v2.0.15^29^, gene functions were indicated using the Nonredundant Protein Sequence Database (NR)^30^, InterPro^31^, UniProt^32^, and EggNOG^33^ with an e-value of 1e-5. In addition, GO annotation was performed by KOBAS^34^ (http://kobas.cbi.pku.edu.cnwas) aligned with the Arabidopsis thaliana database.

Data Records

All the raw sequencing reads of H. omeiensis were uploaded to the NCBI database under accession number SRP435213^35^. The genome assembly had been submitted to Genome Warehouse in China National Center for Bioinformation under accession number GWHEQHK00000000^36^ and European Nucleotide Archive (ENA) with accession number GCA_964187755.2^37^. The annotation files of the genome are available in the figshare database: 10.6084/m9.figshare.22817414.v3^38^.

Technical Validation

Evaluation of the genome assembly and annotation

To assess the integrity of the assembly, short reads were mapped to the genomes using minimap2^39^, giving a mapping rate of 96.59% and a genome coverage of 99.85%. The alignment rate of RNA sequencing reads was 96.95% and 94.10% for two H. omeiensis samples by HISAT2 v2.1.0 (Supplementary Table S3)^24^. The completeness and accuracy of the final genome assembly were checked by Benchmarking Universal Single-Copy Orthologs (BUSCO) v5.4.2^40^ with eudicots_odb10. The results showed that 98.2% of orthologs of eudicots could be identified in the assembly (Supplementary Figure S2). Moreover, the values evaluated by Merqury v1.3^41^ based on short reads also showed high consensus quality (accuracy > 99.99%, QV > 58) and low base-level error rates (1.37 × 10^−6^). In addition, the LTR Assembly Index (LAI) score of the whole-genome assembly was calculated to be 24.52, exceeding that of rice (MSUV7) and Arabidopsis (TAIR10), reaching the ‘gold quality’^42^. These results demonstrated that the assembly is reliable and has high base-level accuracy, high completeness, and high contiguity.

Via multiple annotation approaches, we identified 53,951 protein-coding genes in the H. omeiensis genome (Table 5). BUSCO analysis showed the completeness of predicted genes was 94.5% (Supplementary Figure S2). The functional analysis revealed that 99.7% of the protein-encoding genes could be annotated in at least one of five public databases (Fig. 4).Table 5. Predicted protein-coding genes of H. omeiensis.FeaturesNumberSize mean (bp)% of the genomegene53,9515,145.910.1exon208,603341.62.6intron154,6521,334.47.5mean exons per gene3.9——mean introns per gene2.9——five_prime_utr25,149171.80.2three_prime_utr16,528282.90.2Fig. 4Venn diagram displaying the matches of genes of H. omeiensis in five public protein databases.

Supplementary information

Supplementary Information