Differential gene expression in Olympia oyster (Ostrea lurida) gonadal tissue

Abstract

Introduction

Materials and Methods

Sampling and library preparation

Transcriptome assembly

Transcriptome Annotation

RNA-Seq Analysis

Enrichment Analysis

Results

Transcriptome Annotation

RNASeq Analysis

Differential Expression Analysis

Discussion

References

 

Abstract

Genomic tools provide useful insights to an organism’s physiology, which aid in effective conservation strategies of vulnerable species, such as the Olympia oyster (Ostrea lurida). Establishment of baseline data greatly enhances our understanding of important biological processes and can provide valuable genomic resources for future studies. Gene expression was evaluated using sequenced O. lurida larval transcriptome (121,219 contigs) data. Transcriptome annotation compared known contigs from various genomic databases UniprotKB/Swissprot to the O. lurida transcriptome v3 using BLASTn and BLASTx algorithms, which resulted in a range of 164-35,687 matched contigs. Differential gene expression of male and female gonadal tissue was compared using DESeq analyses and showed 1789 genes expressed significantly higher in female and 495 in male. Genes were categorized by biological processes according to GO terms which led to the discovery of differentially expressed of genes involved in reproductive functions. Results from this study expand current available genomic resources for O. lurida and provide tools for more effective conservation strategies.

 

Introduction

Effective conservation is dependent on understanding an organism’s physiology. This can be done using molecular methods to analyze physiological processes on the genomic level. Genomic tools can provide insights to local adaptation, population structure, genes functions, responses to environmental conditions, immune responses, and many other processes important for conservation strategies (Allendorf et al. 2010). Gene expression analyses in particular are useful for creating more detailed conservation plans because gene expression tells organism’s initial response. Quantitative PCR methods were used to measure gene expression in red abalone under various pH conditions to predict the implications of ocean acidification (Zippay and Hoffmann 2010). Gene expression was also analyzed using whole transcriptome analysis to detect Crassostrea gigas immune response during vibrio infection (Lorgeril et al. 2011). Genomic resources prove to  be an important tool for conservation of vulnerable species.

Transcriptomic data can provide information about which genes are being expressed and further give insight about the interactions between the organism and its environment. Assessing gene expression in various tissues can indicate presence of disease or external stressors that may threaten a population’s viability (Adams 2008). Gene expression has been assessed to find important genes involved in fitness in Pacific salmon to preserve the economically important species under the pressures of climate change (Cooke et al. 2012). Transcriptome data has also been used to monitor physiological processes that have aided Crassostrea gigas adaptation to euryhaline conditions (Meng et al. 2013). Li et al. (2013) used transcriptome data to detect the presence and stages of white spot syndrome in Chinese shrimp. Transcriptome data can be used as a primary method to detect conservation threats to a vulnerable population or species.

 

Analysis of gene expression is useful in detection of an organism’s response environmental changes and population structure, however the addition of baseline data can greatly enhance these conservation targets. Monitoring transcriptomic data over time helps to deceifer evolutionary changes that are occurring (Wildt 2000). Baseline transcriptome analysis can also give insight to evolutionary forces, such as natural selection, acting on a population.  Intraspecific differential expression between sexes may be indicative of evolutionary selection leading to sexual dimorphism (Ellegren and Parsch 2007). Understanding which genes are differentially expressed between sexes gives insights to physiological trade offs that males and females make, which influences overall fitness and is important for conservation strategies for vulnerable species, such as the Olympia oyster.

Beginning in the mid 1800s, the native Pacific Northwest Olympia oyster (Ostrea lurida) population underwent a large decline, possibly due to overharvesting, and since has not fully recovered (White et al. 2009; Polson and Zacherl 2009). Restoration efforts have focused on methods to conserve this important commercial, cultural, and ecological species (Brumbaugh and Coen 2009). Transcriptomic resources are capable of assessing biological processes at the primary level and provide useful tools for understanding environmental interactions. Existing genomic resources for O. lurida include sequenced genome, partially annotated transcriptome and identified primers. Climate change presents a threat to the already declined population, making the establishment of a holistic genomic basis a pertinent issue.

 

This study aims to extend existing knowledge by annotating a newly sequenced transcriptome using multiple databases and to assess differential gene expression of adult male and female O. lurida. Annotation of the transcriptome may allow for connection of genes to their physiological function, which can be used in future studies to assess genomic changes associated with certain biological processes. Differential gene expression in gonadal tissue will give insight to patterns of sexual dimorphic traits. Additional data will extend the available genomic resources for O. lurida conservation.

Materials and Methods

Transcriptome assembly

        In order to obtain a comprehensive transcriptome a de novo assembly was performed using the four PE gonad libraries in combination with two larval libraries previously described (Timmins-Schiffman et al 2013), four additional larval libraries (reference) and four SE libraries corresponding to the samples above. Assembly was performed using Trinity. 

Sampling and library preparation

Gonad samples from two male and two female Olympia oysters in order to  construct four RNA-seq libraries. Briefly, RNA was isolated from using Tri-Reagent (Molecular Research Center). Messenger RNA was isolated using the Ambion Microio88p Poly(A) Purist kit. Library construct and sequencing (PE 50bp) was performed by the University of Washington High Throughput Genomics Unit at the University of Washington, using the Illumina HiSeq platform (San Diego, CA).  

Transcriptome Annotation

Transcriptome annotation was carried out by comparing contiguous sequences using BLAST algorithms. BLASTx and BLASTn (e-value<0.00001) were used to compare the O. lurida transcriptome to nucleotide and protein sequence databases (Table 1). Contigs were classified according to Swiss-Prot Gene Ontology (GO) associations, as well as respective parent categories (GO Slim).

Table 1 Genomic databases used to annotate the O. lurida transcriptome with the corresponding algorithm and source.

Database

Algorithm

Source

Conserved Domain Database (CDD)

BLASTx

http://www.ncbi.nlm.nih.gov/

Environmental Nucleotide Sequences

BLASTn

http://www.ncbi.nlm.nih.gov/

Nucleotide Sequence

BLASTn

http://www.ddbj.nig.ac.jp/

RefSeqGene

BLASTn

http://www.ncbi.nlm.nih.gov/

NCBI Genomic references

BLASTn

http://www.ncbi.nlm.nih.gov/

NCBI Transcript references

BLASTn

http://www.ncbi.nlm.nih.gov/

UniProtKB/Swissprot

BLASTx

http://uniprot.org/

REPRODUCTION  RNA-Seq Analysis

Gene expression in male and female gonadal tissue was analyzed using two methods, DESeq and RNA-Seq, to determine expression values, fold change, and differential expression significance. Transcriptomic libraries from two male and two female gonad libraries were mapped to the assembled transcriptome (v2) using CLC Genomics Workbench RNA-Seq analyses. Read mapping was based on Total Gene Reads with 0.9 minimum read length fraction and 10 max hits per read. RNA-Seq analyses were used to calculate expression values, total gene reads, and fold changes between male and female libraries. Total gene reads were analyzed to determine significant differential expression between male and female using the differential expression analysis of RNA-Seq data with replicates described by Anders and Huber (2012) in the DESeq package of the program R (R Core Team 2013).

Enrichment Analysis

Significant differentially expressed genes determined in DESeq analyses were characterized using DAVID analysis. Expressed genes with an adjusted p-value less than 0.05 were matched to the O. lurida transcriptome version 3 through contig identification. Matches were then described by biosynthetic pathway, functional categories, GO terms, and protein domains. Significantly expressed genes described by GO terms were further enriched using REViGO analysis to establish specific gene function from the Uniprot/Swissprot database.

Results

Transcriptome Annotation

O. lurida transcriptome annotation using the Conserved Domain Database resulted in the highest number of contigs matched at 35,687 or 32.4% of the transcriptome (Table 2). Annotation with Uniprot/swissprot database resulted in the second highest contig matches at 4472 or 4.1% of the transcriptome and contigs identified with gene ontology (GO) terms. A total of 13,292 contigs were associated with GO biological processes; development processes and transport representing the most representation at 21% and 17%, respectively (Fig. 1). A total of 10,081 contigs were associated with GO cellular components; nucleus and plasma membrane representing the most representation at 23% and 21%, respectively (Fig. 2). A total of 3088 contigs were associated with GO molecular functions; nucleic acid binding activity and signal transduction activity representing the most representation at 31% and 28%, respectively (Fig. 3). All other databases resulted in less than 2.6% contigs matched.

Table 2 O. lurida transcriptome annotation results from BLAST analyses.

Database

Contigs Annotated

Percentage Contigs Annotated

Results

Conserved Domain Database (CDD)

35,687

32.42%

CDD Blast Table

Environmental Nucleotide Sequences

496

0.45%

Env_nt Blast Table

Nucleotide Sequence

2647

2.40%

Nt Blast Table

RefSeqGene

164

0.15%

Refseq_gene Blast Table

NCBI Genomic references

464

0.42%

Refseq_genomic Blast Table

NCBI Transcript references

2825

2.57%

Refseq_RNA Blast Table

UniProtKB/Swissprot

4472

4.06%

Swissprot Blast Table

Screen Shot 2014-04-14 at 12.35.20 AM.png

Fig. 1 O. lurida contigs classified by biological processes based on GO Slim categories. Other biological processes and other metabolic processes not shown.

Fig. 2 O. lurida contigs classified by cellular components based on GO Slim categories. Other cellular components, other cytoplasmic organelles, and other membranes not shown.

Fig. 3 O. lurida contigs classified by molecular function based on GO Slim categories.

RNASeq Analysis

RNASeq analysis mapped 56,966 unique gene reads showing 43,003 reads expressed higher in female and 12,879 reads expressed high in male. The other 617 reads have no difference in expression. Plotted fold changes show little variation for the majority of expressed genes (Figure 4). Expressed genes with the greatest fold change were molecular chaperone GroEL, collagen alpha-5 chain, and 3-methyladenine DNA glycosylase  in female (Table 3) and complement C5, actin-1, and DNA translocase FstK 1 in male (Table 4).

Screen Shot 2014-04-08 at 11.15.16 PM.png

Fig. 4 Male and female log expression of individual contigs from gonadal tissue.

Table 3 Ten contigs with highest relative expression in female annotated using UniprotKB/Swissprot database.

Screen Shot 2014-04-11 at 12.42.10 AM.png

Table 4 Ten contigs with highest relative expression in male annotated using UniprotKB/Swissprot database.

Screen Shot 2014-04-11 at 12.42.23 AM.png

Differential Expression Analysis

DESeq differential expression analysis resulted significant expression of 2284 genes with 1789 expressed higher in female and 495 higher in male (Figure 5). Gene annotation showed differential expression of GO biological processes relating to sperm-egg recognition, cell-cell recognition, binding of sperm to zona pellucida, and extracellular matrix organization (Table 5). Revigo analysis determined the majority of genes categorized by biological processes were involved in extracellular matrix organization (Figure 6). Differential expression also occurred in 6 GO cellular components and 14 GO molecular functions (Table 5). DAVID analyses showed a difference in N-Glycan biosynthesis pathway (Figure 6).

Screen Shot 2014-04-11 at 1.07.49 AM.png

Fig. 5 Plot of fold change between female/male expression. Positive fold change indicates reads expressed higher in males and negative fold change indicates reads expressed higher in females. Significant differential expression shown in red.

Table 5 Contigs with significant differential expression identified by GO terms.

Category

Term

p-value

Percentage

GO Biological Processes

Sperm-egg recognition

0.046

2.3

GO Biological Processes

Cell-cell recognition

0.046

2.3

GO Biological Processes

Binding of sperm to zona pellucida

0.046

2.3

GO Biological Processes

Extracellular matrix organization

0.058

3.5

GO Cellular Component

Postsynaptic membrane

0.0069

5.8

GO Cellular Component

Synapse part

0.011

7.0

GO Cellular Component

Extracellular region

0.021

14.0

GO Cellular Component

Cell junction

0.025

8.1

GO Cellular Component

Mitochondrion

0.04

14.0

GO Cellular Component

Synapse

0.043

7.0

GO Molecular Function

Extracellular ligand-gated ion channel activity

0.001

5.8

GO Molecular Function

Ligand-gated channel activity

0.0032

5.8

GO Molecular Function

Ligand-gated ion channel activity

0.0032

5.8

GO Molecular Function

Ion channel activity

0.026

7.0

GO Molecular Function

Channel activity

0.028

7.0

GO Molecular Function

Substrate specific channel activity

0.028

7.0

GO Molecular Function

Passive transmembrane transporter activity

0.028

7.0

GO Molecular Function

Metal ion binding

0.032

33.7

GO Molecular Function

Gated channel activity

0.039

5.8

GO Molecular Function

Cation binding

0.04

33.7

GO Molecular Function

Ion binding

0.044

33.7

GO Molecular Function

Calcium ion binding

0.047

11.6

GO Molecular Function

Neurotransmitter receptor activity

0.081

3.5

GO Molecular Function

Neurotransmitter binding

0.088

3.5

Screen Shot 2014-05-07 at 11.44.09 PM.png

Fig. 6 Biological functions of differentially expressed genes in O. lurida gonadal tissue as determined by REVIGO analysis using GO terms. Size of circles indicates number of genes differentially expressed involved in the process.

Screen Shot 2014-04-14 at 2.25.16 AM.png

Fig. 7 N-Glycan biosynthesis pathway differentially expressed in adult O. lurida gonadal tissue.

Discussion

Genomic resources available for O. lurida have been expanded. Our study provides 110,147 sequenced contigs, increasing previously published transcriptome data from 41,136 sequenced contigs (Timmins-Schiffman et al. 2010). Annotation of the transcriptome connects the contig sequence to its biological function. Determination of biological processes categorized by Uniprot GO terms showed a slightly increased proportion of expressed genes related to developmental processes, although the ratios of other functional classifications were similar to existing data. The addition of contig data does not constitute a fully describe a complete transcriptome, but it does add a wealth of information to genomic resources for O. lurida conservation.

Males showed higher differential expression values, however a greater number of genes were differentially expressed in female O. luida opposed to male. This phenomenon has been observed in other organisms like Drosphila (Meiklejohn et al. 2003). Differential gene expression in gonadal tissue showed a higher expression of genes relating reproductive processes. Analysis of differential expression in gonadal tissues of murine models fathead minnow, and Northern cod also presented increased levels of genes relating to reproductive functions (Small et al. 2005; Filby and Tyler 2005; Nagasawa et al. 2014). Extracellular matrix organization was also observed to be differentially expressed, which has been observed in gonadal expression analysis of other organisms (Knox et al. 1994). The goal of the analysis was to identify which genes are being differentially expressed in gonad tissue, making the sex origin of the expressed gene irrelevant. Detection of differentially expressed genes in gonadal tissue provides a basis for normal expression in mature O. lurida that can be compared to monitor disease presence and environmental stressors.

Discuss the function of select genes in enriched biological processes

Data from this study is highly applicable to future studies focused on conservation of O. lurida. Expansion of the transcriptome sequences allows for the use of novel SNPs that can be helpful population structure analysis and monitoring of physiological changes at various life stages in differing environmental conditions (Pante et al. 2012; Belleghem et al. 2012). Transcriptome assembly has allowed for the detection of SNPs in Crassostrea gigas and Pacific herring to aid in monitoring important physiological signals for conservation purposes (Meng et al. 2013; Roberts et al. 2012). The findings in this study expand existing genomic resources available for O. lurida conservation.


References

Adams J (2008) Transcriptome: connecting the genome to gene function. Nature Education, 1, 195.

 

Anders S, Huber W (2012) Differential expression of RNA-Seq data at the gene level–the DESeq package.

 

Belleghem SMV, Roelofs D, Houdt JV, and F Hendrickx (2012) De novo transcriptome assembly and SNP discovery in the wing polymorphic salt marsh beetle Pogonus chalceus (Coleoptera, Carabidae). PLoS ONE, 7, e42605.

 

Brumbaugh RD and LD Coen (2009) Contemporary approaches for small-scale oyster reef restoration to address substrate versus recruitment limitation: a review and comments relevant for the Olympia oyster, Ostrea lurida, Carpenter 1864. J Shellfish Research, 28, 147-161.

 

Cooke SJ, Hinch SG, Donaldson MR, Clark TD, Eliason EJ, Crossin GT, Raby GD, et al. (2012) Conservation physiology in practice: how physiological knowledge has improved our ability to sustainably manage Pacific salmon during up-river migration. Phil Trans R Soc Lond B Biological Science, 36, 1757-69.

 

Ellegren H and J Parsch (2007) The evolution of sex-biased gens and sex-biased gene expression. Nature Reviews Genetics, 8, 689-698.

 

Filby Al and CR Tyler (2005) Molecular characterization of estrogen receptors 1, 2a, and 2b and their tissue and ontogenic expression profiles in fathead minnow (Pimephales promelas). Bio Reproduction, 73, 648-662.

 

Knox JD, Cress AE, Clark V, Manriquez L, Affinito KS, Dalkin BL, and RB Nagle (1994) Differential expression of extracellular matrix molecules and the alpha 6-integrins in the normal and neoplastic prostate, 145, 167-174.

 

Li S, Zhang X, Sun Z, Li F, and J Xiang (2013) Transcriptome analysis on Chinese shrimp, Feneropenaeus chinensis, during WSSV acute infection. PLos ONE, 8, e58627.

 

Lorgeril J, Zenagui R, Rosa RD, Piquemal D, and E Bachere (2011) Whole transcriptome profiling of successful immune response to Vibrio infections in the oyster Crassostrea gigas by digital gene expression analysis. PLoS ONE, 6, e23142.

 

Meng J, Zhu Q, Zhang L, Li C, Li L, et al. (2013) Genome and transcriptome analyses provide Insight into the euryhaline adaptation mechanism of Crassostrea gigas. PLoS ONE, 8, e58563.

 

Nagasawa K, Presslauer C, Kirtiklis L, Babiak I, and JMO Fernandes (2014) Sexually dimorphic transcription of estrogen receptors in cod gonads throughout a reproductive cycle. J Mol Endocrinol, 52, 357-371.

 

Pante E, Rohfritsch A, Becquet V, Belkhir K, Bierne N, and P Garcia (2012) SNP detection from De Novo transcriptome sequencing in the bivalve Macoma balthica: marker development for evolutionary studies. PLoS ONE, 8, 9.

 

Polson MP, Zacherl DC (2009) Geographic distribution and intertidal population status for the olympia oyster, Ostrea lurida carpenter 1864, from Alaska to Baja. J Shellfish Research 28, 68–77.

 

R Core Team (2013) R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria. ISBN 3-900051-07-0, URL http://www.R-project.org/.

 

Roberts SB, Hauser L, Seeb LW, Seeb JE (2012) Development of genomic resources for Pacific herring through targeted transcriptome pyrosequencing. PLoS ONE, 7, e30908.

 

White J, Ruesink JL, Trimble AC (2009) The nearly forgotten oyster: Ostrea lurida carpenter 1864 (olympia oyster) history and management in Washington State. J Shellfish Research, 28, 43–49.

 

Wildt DE (2000) Genome resource banking for wildlife research, management, and conservation. ILAR, 4, 228-234.

 

Zhang G, Fang X, Guo X, Li L, Luo R, et al. (2012) The oyster genome reveals stress adaptation and complexity of shell formation. Nature, 490, 49-54.

 

Zippay ML and GE Hofmann (2010) Effect of pH on gene expression and thermal tolerance of early life history stages of red abalone (Haliotis rufescens). J Shellfish Research 29, 429-439.

Extra Materials

Ostrea lurida Transcriptome Version 3

Analysis

File Output

RNASeq Analysis

RNASeq Male/Female Expression Output

DESeq Analysis

Swissprot DESeq Output and Swissprot Annotation

Enrichment: Revigo