Differential gene expression analysis

Shrimp
strains and tissue collection

The red, yellow and transparent strains were obtained from Xiamen Fu
shrimp Agricultural Development Co., Ltd., Xiamen, China. All three strains
were cultured in plastic aquarium tanks (30×30×45
cm). Approximately 200 adults were kept in each plastic aquarium tank, and were feed with ornamental shrimp feed (Shirakura, Japan) twice a day at
08:00 and 17:00 h. About 50% water change was carried out once per three days. Few
main water quality parameters remained stable (water temperature 25℃,
dissolved oxygen 7.0-8.0mg/L, pH 7.0-7.4). With the help of anatomy microscope,
the integument of the second abdominal segment was stripped on ice. Other
tissues were removed using 0.7% normal saline solution. The 10 individuals’ integuments
of the same strain were mixed together. The mixture was washed quickly with
0.7% normal saline, blotted on filter paper, put into liquid nitrogen, and then
stored at -80℃ until RNA extraction was performed. There were 9 mixtures from each
strain.

RNA
preparation, cDNA library and and Illumine RNA-seq

Total RNA was prepared from a mixture of 10 individuals using TRIzol
reagent (Invitrogen, USA) according to the manufacture’s protocol. RNA
degradation and contamination were determined using 1% agarose gel electrophoresis.
RNA concentration and purity were checked using a Nano-Photometer spectrophotometer
(IMPLEN, CA, USA), and RNA integrity was assessed using an RNA Nano 6000 Assay
Kit that was provided with a Bio-analyzer 2100 system (Agilent Technologies,
CA, USA). Equal amounts of total RNA from 3 mixtures of the same strain was
mixed to a pool. There were 3 RNA pools of each strain. Each RNA pool was used
to generate a sequencing librarie, and was regarded as a biological sample in
the present paper. 3 μg total RNA per biological sample was used as the input material for
RNA sample preparation. After RNA sample preparation, sequencing libraries were
generated using the NEBNext Ultra RNA Library Prep Kit for Illumina (New
England Biolabs, E7530) following the manufacturer’s recommendations. In
general, mRNA was enriched by oligo(dT) beads. Then the enriched mRNA was
fragmented using fragmentation buffer and reverse transcripted into first
strand cDNA with random primers. Second strand cDNA was synthesized by DNA
polymerase I, RNase H, dNTP and buffer. Then the cDNA fragments were purified with
QiaQuick PCR extraction kit, end repaired, poly(A) added, and ligated to illumine
sequencing adapers. Illumine RNA-seq was carried out by Gene Denovo
Bitochenology Co. (Guangzhou, China) with an Illumina Hiseq 4000
platform.

De
novo assembly of reference sequences

All raw sequencing reads were trimmed by removing reads containing
adapters, reads containing more than 10% of unknown nucleotides, and low
quality reads containing more than 50% of low quality bases (Q-value<10). The
Trinity software was used to assemble the clean reads of 9 libraries using
default parameters (Grabherr et al 2011). The assembled sequences were then
filtered using the CD-hit program to reduce redundancy (Li and Godzik, 2006). Sequences
shorter than 200 bp were discarded. The resulting sequences, which are called
unigenes, were considered as the final non-redundant transcripts and were used
as the reference sequences for subsequent analysis.

Transcriptome
annotation and ontology

The unigenes were used as query sequences to search against the NCBI
non-redundant (NR) protein database (http://www.ncbi.nlm.nih.gov), the Swiss-Prot
database (http://www.expasy.ch/sprot), Kyoto Encyclopedia of Genes and Genomes
(KEGG) database (http://www.genome.jp/kegg) and Cluster of Orthologous Group.
Searches were conducted using the BlastX program with an E-value cutoff of 1e-5.
Gene ontology (GO) was performed by importing the NCBI NR Blast results into
Blast2GO software. GO terms were then assigned to each sequence automatically.
The annotation output was categorized by cellular component, molecular function
and biological process.

Differential
expression gene analysis

Clean reads from each RNA sample were aligned to the reference
transcriptome using short reads alignment tool Bowtie2 by default parameters
(Li et al., 2009). The gene abundances were calculated and normalized to RPKM
(reads per kb per million reads). The expression level of each transcript in
each sample was then normalized using edgeR. Transcripts with fold change
values larger than 2 and p values lower than 0.05, were included in subsequent
analyses as the differentially expressed genes (DEGs) (Robinson et al., 2010).

GO
and KEGG enrichment analysis of DEGs

All
transcripts were mapped to GO terms in the Gene Ontology database (http://www.geneontology.org/ ) and
KEGG database (http://www.genome.jp/kegg/). The GO analysis of transcripts and
highly expressed genes was plotted using the web-based program WEGO program (http://wego.genomics.org.cn/
cgi-bin/wego/index.pl) with default parameters. GO terms with
corrected P-values less than 0.05 were considered significantly enriched. KOBAS
software was used to assess the statistically significant enrichment of DEGs in
the KEGG pathway.