Assembled genomic and tissue-specific transcriptomic data resources for two genetically distinct lines of Cowpea (Vigna unguiculata (L.) Walp)

Plant materials and tissue collection

Cowpea lines IT86D-1010 and IT97K-499-35 were originally sourced from the International Institute of Tropical Agriculture (IITA) and their pedigrees are provided in Supplementary Figure 1. Lines have been maintained in CSIRO for more than 10 generations (not through single seed descent). Material (IT97K-499-35) used in the generation of the reference assembly (Vigna unguiculata v1.0) was sourced from Mike Timko at the School of Medicine, University of Virginia, who had previously received the material from IITA. UC Riverside took the IT97K-499-35 line through single seed descent and confirmed 100% homozygosity before bulking. Analysis is underway to compare the CSIRO lines and UC Riverside lines to quantitatively assess genetic similarity of the independently sourced seed stocks. The plants were grown as described by Salinas-Gamboa et al. (2016). Young unexpanded leaves were collected for DNA and total RNA extraction for both lines. The reproductive calendars developed for these varieties by Salinas-Gamboa et al. (2016) were used to harvest a set of five reproductive tissue types from both IT86D-1010 and IT97K-499-35. Anther tissues containing developing male gametophytes at pollen mother cell, tetrad and mature bicellular pollen stages were pooled to form a pooled male gametophyte (PMG) sample for both lines. In addition, ovules were extracted from both lines from floral buds to provide individual tissue samples containing differentiated megaspore mother cells (MMCs), female meiotic tetrads (FMT), and mature female gametophytes (MFG) at anthesis. Finally, early developing seeds (ES) were collected post-fertilization containing a mixture of zygotes and early globular embryos with proliferating endosperm.

Nucleic acid extraction and sequencing

DNA and RNA extractions were carried out using a Qiagen maxi DNA kit and Qiagen RNeasy plant mini kit, respectively, as per the manufacturer’s instructions. Illumina library preparation and sequencing of DNA and RNA was undertaken by the Australian Genome Research Facility (AGRF) with 2 × 100 bp standard insert paired-end sequencing using a Hiseq 2500 system. Shotgun sequencing libraries from single IT86D-1010 and IT97K-499-35 genomic DNA samples were prepared using the Illumina TruSeq Nano DNA Library Prep Kit as per the manufacturer’s instructions. Whereas the Illumina TruSeq Stranded mRNA Library Prep Kit was used to prepare poly(A) mRNA sequencing libraries from total RNA samples as per the manufacturer’s instructions. Three replicate libraries were prepared and sequenced for each of the RNA samples except for the IT97K-499-35 MMC and FMT samples, where two replicates were sequenced.

Sequence analysis

Raw genomic DNA sequencing reads from IT86D-1010 and IT97K-499-35 were separately assembled into contigs using Biokanga (version 4.3.6) in a multi-step process. First, raw reads were run through ‘biokanga filter’, where common adapter, primer, and vector contaminants were identified and trimmed. Redundant copies of identical paired-end read pairs were removed, and pairs with no sequence overlap to other raw sequence were also removed as they provided no value to the assembly. Filtered paired-end reads were then assembled into contigs, using ‘biokanga assemb’, with default parameterisation that allows 1 base substitution per 100bp of sequence overlap. Resulting contigs were run through a second ‘reassembly’ step with ‘biokanga assemb’, allowing up to 5 base substitutions per 100bp of sequence overlaps to provide reduction in redundant sequence representations. Finally, ‘biokanga scaffold’ added ordering to some contigs, by identifying raw paired-end pairs that match to ends of different contigs under assumptions of sequencing insert fragment size of 110–1500bp. Raw tissue-specific RNASeq reads were separately assembled into transcriptome contigs using Biokanga, with the same multi-step process as used for the genomic DNA reads (above), without the reassembly step to retain putative transcript isoforms. Summary quality metrics for all DNA and RNA sequence reads are provided in the data repository associated with this paper as Supplementary Data 1. These metrics include read length, average GC content and average quality score across the length of the read, the read midpoint and read end.

The assembled genomic DNA sequences of IT86D-1010 and IT97K-499-35 were annotated for predicted gene regions using Augustus v3.1.0 (Stanke & Waack, 2003). From the available Augustus training sets, tomato (Solanum lycopersicum, ITAG2.4) gene sequences were selected in Augustus on the basis of the greater percentage of cowpea RNA reads covered by the resulting gene predictions. Predictions from the Augustus approach also encompassed gene predictions from both DNA strands, partial gene predictions and predictions of untranslated regions (UTRs). The resulting protein sequence predictions, with a minimum length of 100 amino acids, were annotated through matches to the NCBI’s ‘nr’ protein sequence database (downloaded 8^th August 2017) using ‘blastp’ with an e-value threshold of 1e-50.

To complete sequence alignment analysis, the genomic DNA sequencing reads and the tissue-specific RNASeq reads from IT86D-1010 and IT97K-499-35 were pre-processed by ‘biokanga filter’ as described above, prior to alignment with the genomic sequence assemblies of IT86D-1010 and IT97K-499-35 and to the Vigna unguiculata v1.0 reference genome sequences. The software ‘biokanga align’ was used for these alignments and unique-best alignments for each paired-end sequence with an insert fragment size of 100–1000bp to a genomic sequence were reported, with at most 3 base substitutions per 100bp. Auto-end-trimming (read chimera detection) was permitted to 50bp where required. Detection and reporting of SNPs between DNA or RNA sequencing reads and assembled genomes was enabled where there was coverage of at least 5 reads.

Genomic sequence data for IT86D-1010 and IT97K-499-35

A total of 527 and 303 million pair-end DNA sequence reads from IT86D-1010 and IT97K-499-35, were generated, respectively. These were assembled into 39,123 contigs for IT86D-1010 and 57,690 contigs for IT97K-499-35 with average lengths of 15.6 and 9.8 kilobases (kb), respectively (Table 1). The contig assemblies generated were able to incorporate 68–73% of the raw DNA reads generated (Table 2). The majority (>87%) of the assembled genomic contigs from IT86D-1010 and IT97K-499-35 could be mapped to the V. unguiculata v1.0 reference genome (Table 3) with a minimum of 70% contig overlap. When the required contig overlap increased to 90%, contig mapping to the reference decreased to an average of 63% across assembled datasets. Possible causes for lack of mapping at higher stringencies are loss of contiguous alignment or loss of fidelity of assembly towards the end of the survey assembly contigs. In-silico gene prediction identified approximately 60,000 putative coding sequences in both IT86D-1010 and IT97K-499-35 and nearly 70% of these could be annotated to published protein sequences within the NCBI nr public database (Table 4).

Leaf and reproductive cell-type and seed transcriptomes and genomic comparisons

RNA sequencing of the six tissue transcriptomes for each variety generated read counts varying from 125 to 265 million pair-end sequences. These could be assembled into transcript sets varying in size between 35,000 to 74,000 transcript contigs averaging 1 kilobase in length (Table 5 and Table 6). In both cowpea varieties, leaf transcriptomes were the largest in terms of de novo assembled contig numbers and the anther transcriptomes were the smallest. In subsequent analyses RNA sequence read alignment to predicted gene models within the assembled genome resources were used to compare expression counts across tissues. The assembled genome resources for both cowpea varieties provided good coverage for the analysis of RNA sequence reads as approximately 70% of reads across all tissues could be aligned uniquely to all three genomic resources. Transcriptomes derived from IT86D-1010 displayed slightly greater alignment to the IT86D-1010 genomic resource, than the corresponding comparisons for IT97K-499-35 (Table 7). The majority of transcript contigs (80 to 88%) across all tissues in both cultivars could be mapped to the V. unguiculata v1.0 reference genome with a minimum of 70% contig coverage (Table 3). The remaining unmapped percentage could represent a range of scenarios including IT86D-1010 specific contigs, missing regions in the V. unguiculata v1.0 reference genome, tissue-specific extensions to the IT97K-499-35 resource or misassembled transcript contigs. Predicted gene models were considered expressed if they accrued at least 20 uniquely aligning RNASeq reads. In all tissues, approximately 30% of predicted gene models (Table 8) showed expression and 6% of predicted gene models displayed strong tissue-specific expression. We found that on average 90% of IT86D-1010 transcript contigs could be mapped within a IT86D-1010 genomic contig and that the median genomic contig size was 67 kb relative to median transcript contig size of 1.3 kb. This indicates that this resource contains substantial amounts of genomic sequence context around expressed genes in these tissues. This will be important for future explorations of cis-regulatory regions associated tissue-specific gene expression.