Background & Summary

Hulless barley (Hordeum vulgare L. var. nudum) is a monophyletically originated variety of barley that has loose husk cover of the caryopses1 (Fig. 1). While many domesticated barley varieties are hulled and are mainly used for brewing malt and animal feed, hulless barley has been cultivated on a small scale and used as human food because of the ease in processing and edibility1. Although the cultivation of hulless barley is widely distributed, the frequency decreases from east to west2. The most frequently cultivated area is the Tibetan plateau, e.g. Nepal and Tibet, where hulless barley accounts for more than 95% of domesticated barley and is the staple food for people and an important livestock feed. Recently, hulless barley is also increasingly attracting attention as a potential crop for the development of value-added products and multiple food applications3.

Fig. 1
figure 1

Morphology of a Tibetan hulless barley. These pictures show (a) seedling; (b) heading stage; (c) mature stage; (d) filling stage spike; (e) grain of a hulless barley cultivated in Tibet.

Full size image

Obtaining genomic sequences is critical for efficient molecular breeding and understanding of the evolutionary history of crops. Recently, studies, including one from our own group, have made significant progress in sequencing the genomes of hulless barley. Using the short-read sequencing approach, the genomes of two hulless barley strains that were grown in Tibet were sequenced and assembled4,5. The results suggest that many stress-related genes, which were expanded in hulless barley, might have facilitated the adaptation to the high-altitude environment4 and may provide a useful genetic resource for improving barley. Furthermore, by sequencing a population of 437 accessions, a study also showed that the current Tibetan hulless barley cultivars were derived from eastern domesticated barley and were introduced to southern Tibet between 4,500 and 3,500 years ago6. However, due to its large genome size and rich in transposable element sequences (80.8–84%)4,7,8, the genome assemblies of hulless barley using short-read sequencing approach remain incomplete and fragmented. This constrains the molecular breeding in hulless barley and the use of hulless barley for food applications.

Here, using a long-read sequencing technique (Pacific Biosciences), we sequenced a Tibetan hulless barley cultivar (Lasa Goumang) that has been previously sequenced using short-read, in high coverage (>67X). Using both Pacbio long reads and available Illumina reads, we assembled a significantly improved genome, of which the N50 contig size reached ~1.56 Mb (Table 1). Based on this improved assembly, we re-annotated the protein-coding genes in hulless barley and anchored the scaffolds to a linkage map of the barley cv. Morex9. The improved genome assembly and annotations will not only serve as a key resource for exploring the economic and genetic values of hulless barley varieties, but will also advance researches in barley genomics and genetics.

Table 1 Titetan hulless barley genome size estimation and assembly statistics in previous study.
Full size table

Methods

DNA isolation, libraries construction and sequencing

The seedlings were germinated from seeds of a Tibetan hulless barley (cultivar Lasa Goumang, NCBI BioSample ID SAMN09914874). Tissues were flash-frozen in liquid nitrogen and stored in the freezer until DNA extraction. DNA was extracted using the cetyltrimethyl ammonium bromide (CTAB) method10. The quality of the extracted genomic DNA was checked using electrophoresis on 1% agarose gel and the concentration was quantified using a Qubit fluorometer (Invitrogen, Carlsbad, CA, USA).

Single-molecule real-time (SMRT) long-reads sequencing was performed at NextOmics Technology Corporation (Wuhan, China) with a PacBio sequel sequencer (Pacific Biosciences, Menlo Park, CA, USA). The SMRT Bell library was prepared using a DNA Template Prep Kit (1.0). In total, six 20-kb SMRT Bell libraries were constructed. Genomic DNA (~10 μg) was mechanically sheared to fragments of approximately 20 kb using a Covaris g-TUBE. The fragment size distribution was assessed using a bioanalyzer 2100 12 K DNA Chip assay. A blunt-end ligation reaction followed by exonuclease treatment was conducted to generate the SMRT Bell template. The size-selection of SMRT Bell templates was performed using a BluePippin size-selection system (Sage Science) to enrich large fragments (>10 kb). The quality and quantity of size-selected libraries were assessed on a bioanalyzer12Kb DNA Chip (Agilent) and a Qubit fluorometer (Life Technologies), respectively. The SMRT bell libraries were prepared using the binding kit 2.0 (PacBio p/n 100-862-200) according to the manufacturer’s instructions. The libraries were sequenced using a PacBio Sequel instrument on PacBio SMRT cells v2.0 (Pacific Biosciences, acquiring one movie of 360 min per SMRT cell). The MagBead loading (Pacific Biosciences) method was used to improve the enrichment of the larger fragments. In total, ~300 Gb subreads sequences (average length: 9,358 bp) were generated on 64 SMRT cells.

RNA isolation and Iso-Seq sequencing

For RNA samples, plants were grown in a climate chamber in the laboratory (Lhasa). Roots, stems and leaves were sampled seven weeks after germination. To have sufficient materials for each RNA sample, we pooled plant tissue from 10 plants. Samples were placed on dry ice during sample collection and stored in −80 °C freezer until RNA isolation. In total, five pooled samples (one root, two stem and two leaf samples) were collected. Samples were ground with liquid nitrogen and total RNA was extracted using TRIzol reagent (Invitrogen) according to manufacturer’s protocol. RQ1 DNase (Promega) was used to remove DNA.

cDNA libraries were prepared using the ClontechSMARTer® cDNA synthesis kit according to the manufacturer’s recommendations. One μg total RNA was used for each of the five samples. Barcoded oligo dT was used to barcode samples. The cDNA products were purified with AMPure PB beads and quality control (QC) was performed on BioAnalyzer 2100 (Agilent). The purified cDNA libraries were pooled in an equal molar ratio. The pooled cDNA (~3.8 μg) was size fractionated using the Sage ELF system. Subsequent re-amplification was performed to yield four libraries (size of 1–2, 2–3, 3–6 and 5–10Kb) to minimize artifacts during large-scale amplification. The pooled PCR products were purified using AMPure PB beads. One to five μg of purified amplicons were subjected to Iso-Seq SMRT Bell library preparation (https://pacbio.secure.force.com/SamplePrep). A total of 17 SMRT cells were sequenced on the PacBio RS II platform using P6-C4 chemistry with 3–4 h movies. In total, 19.68 Gb sequence data (~1.5 million reads) were obtained (Table 2). The average subreads length was 4.0 kb and the average subreads quality was 0.9.

Table 2 Iso-Seq library information and sequencing results.
Full size table

Genome assembly

In our assembly workflow, raw bam files from PacBio Sequel were first converted into subreads in fasta format using the PacBio software BAM2fastx. Then we used the falcon package (https://github.com/PacificBiosciences/falcon) to construct the primary assembly. Error correction was performed using an overlap-based strategy and the error-corrected reads were used to construct the contigs (parameters: length_cutoff = 5000; length_cutoff_pr = 10000). To correct errors in the primary assembly, we used the arrow pipeline from the SMRT link 4 toolkit to polish the genome (https://www.pacb.com/products-and-services/analytical-software/smrt-analysis/). The PacBio reads were aligned to the primary assembly using pbalign and variantCaller was used to call variants.

SSPACE11 was used to construct scaffolds from contigs. First, we aligned previously sequenced Illumina mate-pair libraries (20 kb and 40Kb fragment long)4 to the assembled contigs using bowtie v1.1.212 and constructed scaffolds with SSPACE-STANDARD-3.0. Second, we used PacBio long-reads to further improve the scaffolding using SSPACE-LongRead13. After scaffolding, the assembly contains 1,856 scaffolds, with a N50 contig size of 1.56 Mb and an N50 scaffold size of 4.0 Mb. The assembled genome size is 4.0 Gb (Table 3).

Table 3 Comparison of the new genome with previously published assemblies of the Tibetan hulless barley genome.
Full size table

While the genome size of barley cv. Morex is ~5.1 Gb9, our previous work has suggested that the genome size of Tibetan hulless barley is ~4.5 Gb using k-mer analysis4 (Table 1). However, it is well-known that genome size estimation from both k-mer approach and flow cytometry can have substantial standard deviations (e.g., 10%)14,15. To draw concrete conclusion on genome size differences between hulless barley and cv. Morex, additional experiments are required. However, this is beyond the scope of this study.

We further generated pseudochromosomes using the assembled scaffolds and linkage map of barley cv. Morex9. We used blastn to map the marker sequences of cv. Morex genome to the scaffolds. Only uniquely mapped markers with coverage greater than 0.8 and identity greater than 0.95 were considered. To anchor the scaffolds to pseudochromosomes, ALLMAPS16 was used (Fig. 2a). The synteny comparison between the newly assembled Tibetan hulless barley and barley cv. Morex (Fig. 3a) was performed using CoGe platform17 (https://genomevolution.org/coge/).

Fig. 2
figure 2

The workflows of genome assembly and annotation used in this study. (a) Genome assembly pipeline; (b) Protein-coding gene annotation pipeline. Software and tools were indicated at lines, data and database information were shown in rectangles.

Full size image
Fig. 3
figure 3

Genome comparison between Tibetan hulless barley and barley cv. Morex. (a) The plot shows the LAST alignments of predicted protein-coding genes in barley cv. Morex assembly and Tibetan hulless barley assembly. (b–d) Dot plots show the MUMMER alignments of Tibetan hulless barley scaffolds and assembled barley cv. Morex bacterial artificial chromosome sequences. Different coverages were shown in B (65%), C (85%) and D (99.9%). The differences can be due to either true sequence divergences or assembly errors.

Full size image

Repetitive sequences annotation

Repetitive DNA sequences are highly abundant in many organisms and their variations in abundancy resulted in remarkable genome size variations in plant18. In many Gramineae crop plants, repetitive elements represent more than 80% of their genome8,19,20. Repetitive elements can be classified as simple repeats and transposable elements (TE). Using tandem repeats finder21, we annotated ~155 Mb (3.89%) sequences as simple repeats. To annotate TE, we used both homology-based and de novo TE annotation tools: RepeatMasker22, RepeatProteinMask22, RepeatModeler (http://www.repeatmasker.org/RepeatModeler/) and LTR_FINDER23. For RepeatMasker, Repbase 21.0124 was used. In total, ~87.5% of the assembled genome were identified as TE (Table 4).

Table 4 Repetitive element annotation statistics.
Full size table

Protein-coding gene prediction

For the annotation of protein-coding genes, we used a previously established gene annotation pipeline (Fig. 2b) with minor modifications4. For de novo gene prediction, we first extracted the full-length transcripts from the Iso-Seq data using SMRTLINK. In total, 39,442 full-length transcripts were obtained and were subsequently aligned to the assembled genome using GMAP25. Among the 38,013 aligned transcripts, we removed all transcripts that had coverage less than 0.9 or sequence identity less than 0.85. This resulted in 14,099 high-quality full-length transcripts, which were used for open reading frame (ORF) prediction by TransDecoder (https://github.com/TransDecoder). In total, 13,936 (98.8%) transcripts contained at least one open reading frame (ORF) that is larger than 50 amino acids. These ORF containing transcripts were then assigned to 9,360 genes, which were considered as authentic genes. The authentic genes were then used for training the gene prediction models using AUGUSTUS v3.2.326. Based on the trained models, AUGUSTUS predicted 128,400 putative genes.

For homology based gene prediction, we used the protein sequences of seven monocot species (Triticum urartu (progenitor of wheat A genome)20, Triticum tauscii (progenitor of wheat D genome)27, Brachy podiumdistachyon28, Hordeum vulgare8, Oryza sativa29, Sorghum bicolor30 and Zea mays31) from public databases. All protein sequences were aligned to the hulless barley genome using tblastn32. The gene structure was predicted using GeneWise33 with the input protein sequence as reference.

To provide further evidence for evaluating the predicted gene models, we assembled the transcriptome using available RNA-seq Illumina short-reads from different libraries34. The transcriptome was assembled using both reference-guided approach (mapping: hisat235, assembly: stringtie36,37) and de novo approach (Trinity pipeline38). The reference-guided approach resulted in 47,490 transcripts and the de novo approach resulted in 722,803 transcripts.

The full-length transcripts from Iso-Seq, assembled transcripts from short-reads were used as evidence to evaluate the predicted gene models using EVidenceModeler39. For the data integration, evidence from different sources was assigned to different weight parameters: 20 for Iso-Seq assembly, 8 for short-reads assembly, 5 for homology-based prediction, 2 for AUGUSTUS gene prediction. In total, 129,269 transcripts were obtained, and the structural optimization was performed using PASA. We removed transcripts that either do not show any homology to sequences in nr or uniport database (blast results with identity > = 50% and coverage > = 50%) or have no protein sequences containing any Pfam40 domain (hmmer results with e-value < = 1e-5). For each gene, only the transcript with the longest protein sequences was kept. The tandem duplicated genes were identified using MCScanX. Genes (4,530) that contain large TE sequences (90% coverage) were discarded. The pipeline generated 61,303 genes.

We further classified these genes into high-confidence (HC) genes (40,457), which are likely true protein-coding genes, and less reliable low-confidence (LC) genes (20,846), which potentially are fragmented genes, pseudogenes and/or non-coding genes. This was done in a two-step procedure as described previously8.

We annotated putative functions of the 61,303 protein sequences using public databases, including nr, KEGG41, SwissProt42, Trembl42, GO43, PFAM40, and InterPro44. Blastp was used to compare the predicted protein sequences with the protein databases(e-value < = 1e-5). Blast2GO45 was used to annotate the GO terms using nr database (downloaded in December 2017) with default parameters. The protein domains were annotated using PfamScan46 and InterProScan47 based on InterPro protein databases, including CATHGene3D48, HAMAP49, PANTHER50, PIRSF51, PRINTS52, ProDom53, PROSITE54, SMART55, SUPERFAMILY56 and TIGRFAMs57.

Non-coding gene prediction

tRNAs were annotated using tRNAscan-SE v1.3.158 and rRNAs were annotated using blastn with the rRNA sequences from Arabidopsis thaliana and Oryza sativa (5S rRNA: AJ307354, 5.8S rRNA: AJ232900, 18S rRNA: X16077, 28S rRNA: AH001750). In addition, we also used INFERNAL to predict the miRNA and snRNA.

Data Records

The genomic Pacbio sequencing data (SRS3725794) and Iso-Seq sequencing data (SRS4809149) are available in NCBI Sequence Read Archive under SRP15912959. The available Illumina genome sequencing data that was deposited under SRP0550424,60 was used in our genome assembly and validation processes. The final genome assembly and annotation was deposited at NCBI GenBank under SDOW0000000061 and NCBI Assembly under GCA_004114815.162. The previously generated RNA sequencing data (deposited in NCBI’s Sequence Read Archive under SRP07487063) was used in our genome annotation steps. All the files in this project, such as the assembled scaffolds, repeat annotation, gene predictions and gene function annotations were uploaded to figshare64.

Technical Validation

We evaluated the quality of the new assembly using three independent approaches. First, we mapped 682.57 Gb previously generated genomic Illumina paired-end reads4,60 to the assembled genome. Overall, ~99.9% of paired-end reads were mapped to the genome concordantly. Second, using the mapped short reads, we estimated the quality value (QV) of the assembly using a previously described method65 in which erroneous bases in the genome assembly were identified based on the variant calling software Genome Analysis Toolkit (GATK)66. The estimated base pair error rate is 6.3E-06, suggesting high accuracy of the assembly at base-pair level. Third, we mapped available RNA-seq Illumina reads to the new assembly34 using bowtie2 (parameters:--sensitive --score-min L,0,-0.1 -k 200 --no-discordant --gbar 99999999 --dpad 0 -p 24 --no-mixed -X 1000 --mp 1,1 -I 1 --np 1). The mapping rate with the new assembly (more than 73.6%) is more than 10% higher than the previous assembly. Fourth, we also mapped the PacBio long reads to the new assembly using blasr (parameters: -m 4 --minMatch 8 --minPctIdentity 70 --bestn 1 --nCandidates 20 --maxScore -500). Overall, 91.2% PacBio long-reads can be mapped back to the assembly. Fifth, we downloaded the assembled bacterial artificial chromosome (BAC) sequences from barley cv. Morex67 and mapped them back to the new hulless barley assembly using mummer68 (Fig. 3b–d). Among 299 BACs that were larger than 200 kb and can be mapped back to barley cv. Morex, 74.6% showed high collinearity with the new hulless barley assembly (coverage greater than 60%). Sixth, we also used BUSCO69 to assess the genome completeness. Among 1,440 conserved eukaryotic core genes, 1,378 (95.7%) were complete, 15 were fragmented and 47 were missing in the hulless barley genome assembly. Together, the results suggested that the newly assembled hulless barley genome is of high quality and will serve as a key resource for future research in barley genetics and genomics.