De Novo Assembly of the Liver Transcriptome of the European Starling, Sturnus vulgaris.

The European starling, Sturnus vulgaris, is a prolific and worldwide invasive species that also has served as an important model for avian ecological and invasion research. Although the genome sequence recently has become available, no transcriptome data have been published for this species. Here, we have sequenced and assembled the S. vulgaris liver transcriptome, which will provide a foundational resource for further annotation and validation of the draft genome. Moreover, it will be important for ecological and evolutionary studies investigating the genetic factors underlying rapid evolution and invasion success in this global invader.

Introduction

The European starling is one of the best-studied avian species, yet no genomic resources existed prior to the recent draft genome (December 2015) and no transcriptome data sets are currently available. Starlings are also an important species in invasion genetics research, due to their global introduction history 1, their invasion success and their known patterns of morphological and genetic variation across introduced ranges 2-4, making them ideal for studies of rapid evolution replicated across invasions. However, these studies require improved genetic resources.

Here we characterise the first starling transcriptome data set, produced using liver tissue from individuals sourced from the range edge in Western Australia and carrying genetic variants previously shown to be under selection 4. Two juvenile male starlings (S274, S290) were collected from Western Australia (S274: S 33.81542, E 120.95964; S290: S 33.79625, E 120.90078). Total RNA was extracted and mRNA libraries prepared following 5. Libraries were barcoded and run together on one lane of HiSeq 2500 (Illumina Inc, San Diego, USA), generating ~230 million 2 x 125 bp paired-end reads.

Raw reads were pooled and processed using Trimmomatic v0.33 6 using the following parameters, ILLUMINACLIP:TruSeq3-PE.fa:2:30:10:4 HEADCROP:13 AVGQUAL:30 MINLEN:36, reducing the dataset to ~36 million paired reads and ~9 million 'orphaned' reads (Table 1). We used both the filtered paired-end and 'orphan' reads in the subsequent de novo assembly with Trinity v2.1.1 7 (using default settings and retaining only those contigs > 300bp) producing 59,557 transcripts, encompassing 48,279 unigenes. After assembly, filtered paired-end reads were mapped to the assembled transcripts with Bowtie2 v2.2.4 8 (82% of reads map back to the de novo assembly; 69 % are properly paired), gene expression was estimated using RSEM v1 2.14 9 (mean transcript expression of 10.7 FPKM) and open reading frames (ORFs) were predicted for all assembled transcripts using Transdecoder (https://github.com/TransDecoder/TransDecoder).

Table 1

Transcriptome assembly and annotation statistics compared to other passerine transcriptomes.

	S. vulgarisa	C. chlorisb	Z. albicollisc
Raw sequencing reads	230403632	nr	-
Reads used in assembly	45309889	~500000000	-
Unfiltered de novo assembly
Number of unigenes	48279	nr	-
Number of transcripts	59557	66072	313060
n50 transcript length (bp)	1765	803	3979
sum transcript length (Mb)	64993660	39395826	334636954
median transcript length (bp)	626	367	345
mean transcript length (bp)	1091	596	1069
GC %	48.28	46.34	45.43
Filtered de novo assembly
Number of unigenes	18678	-	-
n50 longest unigene/transcript	2232	-	-
Sum longest unigene/transcript	26825376	-	-
Median longest unigene/transcript length (bp)	979	-	-
Mean longest unigene/transcript length (bp)	1436	-	-
Number of transcripts	23945	-	-
n50 transcript length (bp)	2328	-	-
Sum transcript length (Mb)	37637538	-	-
Median transcript length (bp)	1178	-	-
Mean transcript length (bp)	1572	-	-
Annotation statistics
Unfiltered transcriptome
Transcripts with Blastx match	33041 (55%)	nr	-
Transcripts with Blastp match	24715 (41%)	23,151 (35%)	-
Transcripts with GO terms	27576 (46%)	nr	-
Filtered transcriptome
Transcripts with Blastx match	19701 (82%)	-	-
Transcripts with Blastp match	17898 (75%)	-	-
Transcripts with GO terms	17462 (73%)	-	-

aThis study; bdata from Meitern et al. (2014); cdata calculated from NCBI GBBC00000000.1; nr, not reported; percentage in parentheses

Transvestigator v0.1 10 was used to filter and prepare the raw transcriptome assembly for NCBI Transcriptome Shotgun Assembly (TSA). Briefly, we filtered poorly supported and lowly expressed transcripts if they had a transcript per million (TPM) value < 0.5 (7,873 transcripts) and an isoform expression level less than 5% of the parent unigenes expression (7,708 transcripts). Additionally, only those transcripts passing the above filter and containing a predicted ORF were kept (18,103 transcripts removed). Assembled sequences were queried against the UniVec database (Accessed January 2017) and any sequences with a strong match were removed (1/1,000,000 chance of a random match for queries of 350 Kb, terminal match score ≥ 24, internal match score ≥ 30). Lastly, preparation for TSA submission included ensuring there was only one ORF per transcript, the ORF was on the positive strand and contained within transcript coordinates, and that start and stop codons were properly created. The final filtered transcriptome contains 23,945 protein-coding transcripts from 18,678 'unigenes' (accessible via the NCBI TSA accession: GFDQ00000000). Overall, the filtering resulted in the removal of 1,545 transcripts that may have contained some protein-coding potential. As non-coding RNAs have key biological roles (rRNA, tRNA etc), may be significant in terms of adaptive processes and are useful for genome annotation we also provide the unfiltered transcripts (includes those not passing ORF filters, but meeting support and expression thresholds and those with ORFs not meeting expression thresholds) in the associated data repository 13. Functional annotation was conducted utilising the Trinotate pipeline (http://trinotate.sourceforge.net) - these are also provided in the associated repository.

To assess the relative quality and completeness of the assembly, we compared core assembly statistics and evaluated completeness using the BUSCO (benchmarking universal single-copy orthologs) vertebrate gene set 11 to two recent passerine transcriptomes available in the NCBI TSA that also used the Trinity assembler: Carduelis chloris, GBCG00000000.1 12 and Zonotrichia albicollis, GBBC00000000.1 (Table 1, 2). Our assembly compares favorably in terms of core statistics (Table 1). Similarly, BUSCO assessment also reveals comparable completeness (S. vulgaris, 50% complete; C. chloris, 30% complete; Z. albicollis, 62% complete; Table 2) even with a much smaller number of transcripts assembled.

Table 2

BUSCO evaluations of completeness against the vertebrate gene set compared to other passerine transcriptomes.

	S. vulgaris	C. chloris	Z. albicollis
Complete	1523 (50%)	904 (30%)	1860 (62%)
Single	1409 (47%)	900 (30%)	1563 (52%)
Multi	114 (4%)	4 (~0%)	297 (10%)
Fragment	346 (11%)	361 (12%)	249 (8%)
Missing	1154 (38%)	1758 (58%)	914 (30%)

The recent S. vulgaris draft genome assembly (NCBI Bioproject: PRJNA304638) includes a predicted transcriptome produced by the NCBI Eukaryotic Genome Annotation Pipeline. We compared the final de novo assembled liver transcriptome set (valid ORF and evidence of expression) to the draft genome and predicted transcriptome to aid characterization of gene features and validate gene models through both a standard and reciprocal best-hit BLAST approach and mapping of our mRNA reads to the genome. As our transcriptome assembly is based on a single tissue we only expect a partial representation of the predicted genes to be confirmed. The standard BLASTn approached revealed 23,557 (98.4 %) of our final transcripts had a significant match (e-value <10 -5) to the genome, with 20,768 (86.7 %) having a significant match to the predicted transcriptome. The reciprocal best-hit approach reduced the number of significant matches to the predicted transcriptome to 10,265 (42.9%), which may in part be due to predicted protein coding genes mapping to multiple isoforms present in our transcriptome assembly. While 11.7% (2,790) of our final transcripts did not produce a significant match to the predicted transcriptome they did produce a significant match to the draft genome assembly - underlying the importance of including RNA-Seq data in the annotation process. We used bbmap v35, (https://sourceforge.net/projects/bbmap/) with default settings to map the mRNA reads used to construct the de novo assembly to both the draft genome and predicted transcriptome, with 84.2% and 69.1% of reads mapping respectively. The higher percentage of reads mapping to the genome indicates there may be missed features in the predicted transcriptome. Additionally, standard BLASTn of the unfiltered de novo assembly to the genome revealed that 58,192 (97.8%) transcripts had a significant match, whereas only 38,021 (63.8%) had a significant match to the predicted transcriptome. This gives further support to the previous conclusion of missing features in the predicted transcriptome and suggests that our filtering (only keeping transcripts with a valid ORF and evidence of expression) may have removed some true transcripts (potential ncRNAs) that are encoded in both the genome and predicted transcriptome. Consequently, we have included these transcripts in the associated public repository 13.

This resource will be useful for the validation of gene models predicted in the recent draft genome assembly and for future research into the physiology and ecology of this and other closely related species. In particular, these data will enable a greater understanding of the genetic factors underlying rapid evolution and invasion success in this global invader.

Availability of supporting data

The datasets supporting the results presented here are available at http://dro.deakin.edu.au/view/DU:30091025 and doi: 10.4225/16/5893999965ac2 13. All raw sequencing data used in this study is available in the NCBI SRA and associated with the BioProject accession, PRJNA335913. The final transcriptome assembly has been deposited at DDBJ/EMBL/GenBank under the accession GFDQ00000000. The version described in this paper is the first version, GFDQ01000000.