De novo Assembly and Analysis of the Chilean Pencil Catfish Trichomycterus areolatus Transcriptome.

Trichomycterus areolatus is an endemic species of pencil catfish that inhabits the riffles and rapids of many freshwater ecosystems of Chile. Despite its unique adaptation to Chile's high gradient watersheds and therefore potential application in the investigation of ecosystem integrity and environmental contamination, relatively little is known regarding the molecular biology of this environmental sentinel. Here, we detail the assembly of the Trichomycterus areolatus transcriptome, a molecular resource for the study of this organism and its molecular response to the environment. RNA-Seq reads were obtained by next-generation sequencing with an Illumina® platform and processed using PRINSEQ. The transcriptome assembly was performed using TRINITY assembler. Transcriptome validation was performed by functional characterization with KOG, KEGG, and GO analyses. Additionally, differential expression analysis highlights sex-specific expression patterns, and a list of endocrine and oxidative stress related transcripts are included.

Introduction

Trichomycterus areolatus (described by Valenciennes, 1846) is a threatened pencil catfish found in Chile 1, 2. This freshwater stream fish has a broad distribution ranging from the Coquimbo to the Los Lago regions of Chile (30° S and 43° S respectively) 1, 3. The genus Trichomycterus is found throughout freshwater ecosystems in both Central and South America and includes over 120 species 4. Individual species are commonly restricted to specific river areas either by geographical restriction or habitat preference 5. As a result, T. areolatus displays an intricate integration with its environment; yet the species is fairly distributed throughout the region, providing an opportunity for its use as a model sentinel organism to study environmental stress in Chile's freshwater streams 6,7. As a benthic fish, T. areolatus displays interactions with riverine sediments and substrates 8. This species tends toward a generalist eating strategy, eating aquatic insects 9 as well as the surface organisms present on stream, plant, and rock surfaces 10; consequently, its food preferences are known to vary by season 4.

The Trichomycterus fish generally range between 50 - 150 mm 11. The mature fish in this species will reach average sizes of 56 and 51mm for males and females respectively 12. They display poor secondary sex characteristics, making the identification of sex by external means difficult 13. The spawning season takes place during October and November, based on observation of the female gonads rapidly increasing in size just before the season begins 13. The females are capable of releasing eggs several times throughout the spawning season. The males also exhibit numerous fertilization events throughout the spawning season; however, the male fish in this species do not participate in progeny care 14. The T. areolatus genome is described as diploid, with the typical individual displaying 2n = 54 chromosomes; however, intra-individual variation has been reported 15. This chromosome number appears common among the Trichomycterus genus and suggests conservation.

Genomes and/or transcriptomes have only been developed for fish that have either economic or scientific value including: Salmo salar (Atlantic salmon), Cyprinus carpio (common carp), Pimephales promelas (fathead minnow), and Danio rerio (zebrafish). There are limited transcriptomic resources available for fish native to South America, and none available of the Trichomycterids. Current molecular research with this organism is limited, and includes a karyotype 15, seasonal variation in biomarkers 8, local industrial discharge impacts 6, agricultural disturbances 16, microsatellite loci for conservation resources 17, and population genetics 18.

In this paper, we detail the assembly and analysis of the transcriptome for T. areolatus. Fish samples obtained from Chile were sequenced by next-generation sequencing and later assembled into a de novo transcriptome using Trinity. Transcriptome validation analyses are included and discussed. Transcriptomic data is deposited at NCBI SRA under accession SRP077018 and the completed assembly is freely available.

Results and Discussion

Transcriptome Characteristics

The constructed transcriptome assemblies are available at http://www.davislab.net/trichomycterus/. Assemblies available include: a complete transcriptome nucleotide assembly, a representative assembly (see Methods) consisting of 64,385 unique transcripts, and a translated protein assembly of the representative transcripts. De novo assembly was accomplished using TRINITY. Two samples (whole male and whole female adult fish) provided RNA for the assembled transcriptome, which was assembled from Illumina® paired-end reads. A resulting 41.8Gb of output (Table 1) was utilized to create the assembly. General statistics of the assembled representative transcriptome are included in Table 2.

Table 1

Transcriptome Tissue Sequencing Details. RNA samples were from 4.4 and 5.4ug of RNA for male and female fish, respectively. Sequencing was performed on an Illumina® Hi-Seq 2500. BP: base pair, GC: G-C nucleotide ratio.

Tissue	Total Reads	Total Output (bp)	GC Content (%)
Whole Female	328,721,780	32,872,178,000	48%
Whole Male	88,794,542	8,879,454,200	47%

Table 2

. The representative Trichomycterus areolatus transcriptome assembly was analyzed for general characteristics listed above. Putative protein coding transcripts were included and identified by TransDecoder. Redundant transcripts were removed by CD-HIT which collapses redundant and highly similar sequences into consensus sequences.

Total Transcripts	Mean Length (bp)	Median Length (bp)	N50	GC Content
64889	1484.85	857	2671	47.5%

Organism Phylogenetics

Assembled T. areolatus sequences were selected and aligned to publically available sequences from various fish to examine phylogenetic relatedness by identifying orthologs (Table 3). Of the fish compared, Ictalurus punctatus (Channel catfish) showed the highest relationship (87.6%) in a concatenated set of conserved sequences, followed closely by Cyprinus carpio (Common carp) at 86.5%. Intra-genus phylogenetic comparison of Trichomycterus was not possible due to a lack of published sequences.

Table 3

Phylogenetic Comparison to Other Fish. Transcripts produced in this study were concatenated and aligned to published sequences of fish species and a percent identity matrix was computed. This analysis utilized Trichomycterus areolatus transcripts: Ta_155828, Ta_53325, Ta_192266, Ta_56196, and Ta_56194 for cytochrome c oxidase subunit III, HSP70, NADH dehydrogenase subunit 5, estrogen receptor, and glutathione s-transferase kappa 1, respectively.

Species	Percent Identity
Ictalurus punctatus (Channel catfish)	87.6
Cyprinus carpio (Common carp)	86.5
Pimephales promelas (Fathead minnow)	84.5
Oncorhynchus mykiss (Rainbow trout)	82.5
Salmo salar (Salmon)	82.5
Carassius auratus (Gold fish)	76.4
Danio rerio (Zebra fish)	67.1

Assessing Full-length Transcript Coverage

In an effort to deduce the full-length nature of the transcripts within the representative assembled transcriptome, coverage histograms arising from alignments with non-redundant protein sequences of two model organisms—Salmo salar and Danio rerio—were produced. Model organisms were selected over more closely related organisms (Table 3) due to their more complete publically available transcriptomes. Figure 3 illustrates a histogram for the distribution of transcript length when the representative T. areolatus transcriptome was compared to related species. Length coverage exceeding 90% of the other organism's transcript length was found in 64.7% of T. areolatus protein sequences upon alignment with Salmo salar proteins and 69.1% of Danio rerio proteins (Figure 3). Accounting for genetic differences between these organisms, these data suggest that the produced representative transcriptome has a high degree of full-length transcripts.

Figure 3

Transcript Coverage of Two Model Organisms. Coverage of Salmo salar and Danio rerio predicted proteins by Trichomycterus areolatus predicted proteins. Predicted polypeptide sequences produced in this study were BLASTed against publically available non-redundant Salmo salar proteins (count = 112,089) and Danio rerio proteins (count = 81,931). The length of the local alignment region reported by the BLASTp algorithm was subsequently divided by the length of the query sequence. Compilation of these results indicated that a vast majority of Trichomycterus areolatus predicted protein sequences exhibited greater than 90% coverage of both Danio rerio (64.7%) and Salmo salar (68.9%) protein sequences, suggesting that the assembly produced a high degree of full-length transcripts.

Functional Analyses - Putative Transcript Functional Characterization

Multiple analyses exist that allow transcripts to be annotated and grouped by function. This includes Gene Ontology (GO), Eukaryotic Orthologous Groups (KOG), and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. These annotation and grouping algorithms were applied to examine the putative function of transcripts, and as a quality control technique to evaluate transcriptome completeness when compared to well-developed transcriptomes from Danio rerio, Salmo salar, and Cyprinus carpio. Figure 4 details the GO analysis, which is capable of classifying predicted gene products by cellular component, molecular function, and biological process 19. GO groupings between organisms suggest that the selected fish all share a pattern of functional distribution of gene products.

Figure 4

Gene Ontology (GO) Analysis of the GO functional analysis was performed on assigned proteins in order to evaluate transcript function and the overall completeness of the isolated transcriptome. GO terms were given for each of the T. areolatus predicted proteins as well as the proteomes of Salmo salar, Cyprinus carpio, and Danio rerio (retrieved from NCBI). The distribution of protein functions closely match one another, suggesting the assembled transcriptome is complete.

Similarly, Figure 5 details the KOG transcriptomic analysis; this database allows the user to compare against a collection of seven eukaryote genomes of a diverse set (e.g. human, fly, parasite, plant, worm, fungi) of model organisms 20. These annotated proteins are organized by function into clusters of eukaryotic orthologous groups. Orthologs typically have similar functions among different organisms and serve as an effective means of identifying putative functions of gene products 20. The KOG analysis of T. areolatus resulted in KOG groups classifying 81% of transcripts. Finally, KEGG analysis was performed to identify biologically relevant pathways of function for predicted protein products from the representative transcriptome 21. T. areolatus was compared to Danio rerio, Salmo salar and Cyprinus carpio as is shown in Figure 6. Distribution of group numbers also appears consistent as related fish were evaluated. Overall, each annotation and grouping algorithm showed close association with T. areolatus and the related interrogated species, providing evidence of the completeness and consistency of the transcriptome.

Figure 5

Eukaryotic Orthologous Groups (KOG) Characterization of Putative transcript functions were assessed and transcriptome completeness was evaluated using KOG analysis. The Trichomycterus areolatus transcriptome and mRNA nucleotide entries from NCBI of Cyprinus carpio, Salmo salar, and Danio rerio were assigned KOG terms. The three transcriptomes have similar distributions, supporting the completeness of the Trichomycterus areolatus transcriptome.

Figure 6

Kyoto Encyclopedia of Genes and Genomes (KEGG) Transcriptomic Analysis. KEGG analysis was performed to functionally describe transcript functions and evaluate transcriptome completeness. To serve as comparisons mRNA sequences for Danio rerio, Salmo salar, and Cyprinus carpio were retrieved from NCBI and characterized into KEGG pathways. The percent distribution shows a similar proportion among compared species indicating a complete transcriptome for Trichomycterus areolatus.

Differential Expression Analysis

Table 4 details unique transcripts which were differentially regulated between male and female (sex determined by necropsy) T. areolatus samples. Non-inherited transcripts that are unique to individual organisms (e.g. MHC molecules via rearrangement) were excluded.

Table 4

Differentially Expressed Transcripts. The resultant male and female Trichomycterus areolatus sequence files were interrogated to assess relative differential transcript expression. Unique transcripts, demonstrating changes greater than or equal to 10-fold, are identified by homology (if available) to known proteins, and only the most differentially expressed isoform is presented. Non-inherited transcripts that are unique to individual organisms (e.g. MHC molecules via rearrangement or similar immune transcripts) were excluded. Notably, many transcriptional differences are related to sex-specific expression.

Transcriptomic ID	Fold Change	Name
Male differentially upregulated transcripts
TRICH01_163992	15	parvalbumin beta-1
TRICH01_133698	12	hemoglobin subunit beta-2
TRICH01_83415	11	sperm acrosome membrane-associated protein 4
TRICH01_57563	10	endonuclease domain-containing 1 protein
Female differentially upregulated transcripts
TRICH01_101761	38604	microsomal triglyceride transfer protein large subunit
TRICH01_222696	96	complement C4
TRICH01_180471	96	alpha-2-macroglobulin
TRICH01_143497	56	3-hydroxyacyl-CoA dehydrogenase type-2
TRICH01_54871	42	vitellogenin 4
TRICH01_100604	38	pyruvate dehydrogenase phosphatase regulatory subunit
TRICH01_142598	35	CD59
TRICH01_51308	29	ribonucleoside-diphosphate reductase subunit M2
TRICH01_208172	22	coiled-coil domain-containing protein 36
TRICH01_51310	21	Jouberin
TRICH01_100599	19	PEX5
TRICH01_100601	17	histone-lysine N-methyltransferase ASH1L
TRICH01_100600	16	A-kinase anchor protein

Among the most highly differentially regulated transcripts, the gene product of A2-macroglobulin (A2M) functions as a protease inhibitor and binds growth factors 22. A2M, increased in the female by 96-fold, has been shown to increase endogenous production of estradiol resulting in increased follicular cell proliferation and oocyte maturation 22.

Vitellogenin is an egg yolk lipoprotein precursor that serves as an essential material for oocyte development 23, 24, 25. Vitellogenin production is under estrogenic control, and was upregulated in the female T. areolatus 42 fold or more when compared to the male. The female liver (the primary vitellogenin production site) secretes the lipoprotein, which is then transported from the systemic circulation into oocytes within the developing ovary. Because vitellogenins are large lipoproteins, mobilization requires the microsomal triglyceride transfer protein for delivery to the oocyte. Microsomal triglyceride transfer protein is crucial in final yolk lipoprotein assembly 26, and is seen upregulated in the female by 38,604-fold.

Throughout the process of vitellogenesis, the female organism undergoes metabolic changes, shifting energy usage, in order to support the growth of newly forming oocytes. Moreover, lipid metabolism is essential during ovarian development, including lipid storage, oxidation, and as previously mentioned, mobilization. Under some conditions such as starvation, the female may utilize stored lipids for transfer to the growing oocyte, or subsequent oxidation to meet dietary needs when intake is low 27. The increased expression (56-fold) of 3-hydoxyacyl-CoA dehydrogenase suggests increased lipid breakdown and usage in the female most likely stimulated by increased energy requirements during vitellogenesis.

Similarly, the transcriptomic profile of the female demonstrates increased glycolytic activity, likely also to support the energy-intensive vitellogenesis. The pyruvate dehydrogenase complex is a highly regulated system that connects glycolysis with the TCA cycle. Activation of the complex is enhanced by pyruvate dehydrogenase phosphatase (PDP) 28. The observed 38-fold increased expression of pyruvate dehydrogenase phosphatase would promote pyruvate progression through the TCA cycle and yield additional cellular energy in the form of ATP.

The complement system in fish is a primary component of their innate immune system and is suggested to be more beneficial to organism defense than that of mammalian systems 29. Complement proteins are produced by the liver in an inactive form, where they are activated by proteolysis, ultimately leading to either opsonization, phagocytosis, or lysis of the pathogen. Vertebrate fish have multiple isoforms and isotypes of complement proteins, enabling the organism's immune system to recognize a wide range of pathogens, which is crucial due to slow lymphocyte proliferation, and limited antibody production and affinity in the fish immune system 30,31. Notably, catfish as well as other teleost species, practice external fertilization which may expose the developing embryo and developing fish to waterborne threats, including pathogens. Maternal complement proteins are transferred from the female to the egg and serve to protect the embryos until the immune system and lymphoid organs are competent enough to protect the developing fish 32. Thus, it is expected that maternal complement proteins, especially CD59 which acts as a stabilizer to prevent premature complement activation, would be upregulated. Accordingly, the CD59 gene is upregulated by 35-fold in the female. However, seasonal and temperature variations cause complement proteins to vary in expression, accounting for some elevated complement proteins in the male 33.

Relative increases in male-associated transcripts are also presented in Table 4. Parvalbumins are calcium binding proteins found in white, fast twitch skeletal muscle of most fish species. In muscular tissues it acts to sequester calcium, accelerating muscle relaxation. High expression of parvalbumins can promote quicker muscle relaxation, however overexpression leads to smaller mitochondrial densities in slow twitch muscles 34. More recently, the beta-1 parvalbumin isoform in carp seminal plasma has been characterized 35, and while the specific details of its function remain undefined, it is thought to play a role in sperm motility. Studies in carp have shown that sperm are not mobile without calcium, and initiation of motility is seen with an increase in intracellular calcium levels. The influx of calcium has been suggested to be the initiating factor of sperm motility, and therefore, the presence of parvalbumin as a calcium binding protein may play an essential role in fish sperm function. Additionally, studies in mammalian organisms have suggested parvalbumin regulation in calcium-mediated spermatogenesis and testosterone production 36. Accordingly, the male transcriptome demonstrates a 15-fold increase in beta-1 parvalbumin expression.

Though many studies do not identify significant sex differences of hemoglobin levels in fish 37, Falahatkar et al. 38 discovered decreased hemoglobin concentrations and hematological changes in juvenile Acipenser stellatus that were treated with estradiol for 7-9 months. They predict estradiol may have an inhibitory effect on erythropoiesis (the production of red blood cells). Hemoglobin transcription levels in the male were found to be 12-fold upregulated compared to the female.

During mammalian fertilization, the acrosome surrounding the sperm head releases acrosomal hydrolases that enables the sperm and egg to combine 39. The acrosomal-egg interaction, and release of hydrolytic enzymes, are thought to be stabilized by the outer acrosomal membrane matrix interactions, and the specificity of the binding interactions are essential in release of hydrolases 39. The sperm acrosome membrane associated protein 4 is retained following the egg-sperm binding, localizing to the inner acrosomal membrane in human sperm, and may play a role in fertilization. Of note, teleost organisms differ in fertilization strategies: their sperm do not possess outer acrosomal membranes and therefore, do not undergo the same acrosomal reaction seen in mammals. Interestingly, recent catfish transcriptome studies have identified a homolog of sperm acrosome membrane-associated protein 4 40, though the specific function in catfish has not been elucidated 41. A homolog to acrosome membrane associated protein 4 was observed in this male, expressed at 11-fold higher compared to the female.

Consideration as a Sentinel Organism

This transcriptome provides a resource to utilize T. areolatus as a sentinel organism or a “canary in the coal mine” for biological effects that may be experienced by local wildlife and nearby human populations. In North America, gene expression biomarkers have been readily applied to environmental monitoring for anthropogenic pollutants 42,43,44,45. For example, fathead minnows are commonly used environmental sentinel in studies on agricultural runoff 46,47,48, waste water treatment plant effluent 49,50 and industrial waste effluent 51. The growing reliance on transcriptomic tools in the field of environmental toxicology has been due to their increasing availability for non-model organisms as well as the mechanistic insight they provide for prediction of adverse outcomes at the whole organismal and possibly population level 52,53.

Chile is experiencing significant economic growth driven by agricultural and industrial development which puts considerable pressure on the freshwater resources across much of the country 54,55. However, scientifically documented adverse impacts on water quality are largely absent, particularly with respect to biologic impacts on sentinel species. Thus, a need exists to develop environmental sentinel organisms suitable for study across the country. T. areolatus is an ideal candidate, as it can be found across much of the range of many rivers in Chile from the mountain watersheds to the low elevation, coastal regions of the country. Furthermore, the fish can be found associated with the benthos in streams where the water is merely a few centimeters in depth. Finally, the previously mentioned benthic behavior of the fish causes it to be intimately associated with sediments, therefore it is susceptible to contaminants that are affiliated with the water and sediment 56. To provide a resource for the study of environmental systems associated with this organism, specific genes were identified for convenience (Table 5). These selected genes of interest will serve as markers of reproductive dysfunction and/or oxidative stress, in an upcoming manuscript currently in progress.

Table 5

Genes linked to endocrine disruption and/or oxidative stress were identified within the transcriptome assembly for convenience in developing Trichomycterus areolatus as an environmental sentinel organism. Danio rerio sequences were used as queries to BLAST the full assembly to identify putative homologs. Protein isoforms were differentiated based on query sequence annotation and bitscore.

Gene Name	Gene Symbol	Transcriptomic ID	Transcript Length (bp)	Query ID	Bit Score
Androgen Receptor	AR	TRICH01_58265	4455	NP_001076592.1	788
Aromatase	CYP19a1	TRICH01_14384	284	AAB65788.1	759
Aryl Hydrocarbon Receptor	AHR	TRICH01_166983	3366	NP_001019987.1	593
Aryl Hydrocarbon Receptor 2	AHR2	TRICH01_225654	2827	NP_571339.1	905
Cytochrome P450 1A1	CYP1a1	TRICH01_117246	2028	NP_571954.1	820
Estrogen Receptor Alpha	ESRa	TRICH01_56196	4380	AAK16740.1	729
Estrogen Receptor Beta 1	ESRb1	TRICH01_211240	4700	CAC93848.1	673
Estrogen Receptor Beta 2	ESRb2	TRICH01_95151	3578	CAC93849.1	796
Follicle Stimulating Receptor	FSHR	TRICH01_111037	3596	AAP33512.1	996
Forkhead Box L2	FOXL2	TRICH01_143206	1866	AAI16586.1	370
Heat Shock Protein 70	HSP70	TRICH01_53325	2625	AAF70445.1	1216
Heat Shock Protein 90 Alpha 1	HSP90a1	TRICH01_121146	2863	NP_571403.1	1292
Heat Shock Protein 90 Alpha 2	HSP90a2	TRICH01_121144	2926	AAI63166.1	1278
Metallothionein	MT	TRICH01_130427	554	AAS00513.1	53
Superoxide Dismutase	SOD	TRICH01_28552	2500	NP_571369.1	261
Thyroid Receptor Alpha	THRa	TRICH01_196629	2500	AAA99811.1	760
Thyroid Receptor Beta	THRb	TRICH01_20904	2406	AF109732_1	732
Vitellogenin 1	VTG1	TRICH01_101739	2851	AF406784_1	1224

The availability of a complete transcriptome for T. areolatus provides a valuable resource relative to the development of this species as a sentinel organism. The unique niche that T. areolatus fills in Chile will allow it to be useful when considering water and sediment contamination.

Methods

Tissue Collection and RNA Preparation: Whole fish were collected from the Choapa River basin in the Coquimbo region (region VII) of Chile in July 2015 under fishing authorization #2017 by the Chilean Subsecretary of Fisheries and Aquaculture. Specifically, a whole male and female with clear sexual differentiation were sampled from a downstream and upstream site respectively. The river sampling site coordinates are shown in Figure 2. This watershed is used intensively for agricultural production.

Figure 2

Choapa River Basin Tissue Sampling Sites-The male and female fish samples were collected from the downstream site “A” (altitude 243m, Lat; Lon: -31.749639; -71.160722) and upstream site “B” (altitude 792m, Lat; Lon: -31.89675; -70.783056) respectively. This river basin is proximate to heavy agricultural practice and downstream of heavy metals mining (e.g. copper). The stream itself is predominately supplied by glacial melt from the Andean mountains. Historically, this system drained into the Pacific Ocean but is becoming an isolated system due to the deleterious effects of climate change.

Fish samples were prepared and immediately submerged in RNAlater® (Ambion) according to specific manufacturer recommendations to preserve RNA integrity. The samples were mechanically homogenized in the presence of Qiagen Lysis Buffer RLT; immediately following, a Qiagen RNeasy Mini Plus isolation kit was used to isolate and purify organism total RNA. The resultant RNA was quantified by Thermo Scientific™ NanoDrop 2000c and verified for integrity with a bleach denaturing agarose electrophoresis gel 57. Prior to sequencing, the purified RNA was stored at -80° Celsius with minimal handling and freeze-thawing cycles.

High-throughput sequencing: Prior to sequencing, a TruSeq® RNA Sample Preparation Kit was used to prepare the library for sequencing. An Illumina® HiSeq2500 next generation sequencer was used to generate paired-end 101bp reads. Sequencing was performed at the University of Nebraska Medical Center sequencing core. The whole female and male tissue reads were used exclusively in transcriptome construction.

Data Processing and Assembly: The resulting sequence reads were first processed with FastQC 58 to evaluate sequence quality. Next, PRINSEQ was used to both trim reads by quality score, and apply a dusting technique 59. Transcriptome assembly was performed de novo with TRINITY assembler 60, 61, 62. The resultant initial transcriptome assembly was searched for predicted coding sequences at least 50 amino acids long using TransDecoder 61, BLASTed against NCBI Refseq metazoan protein sequences from March 2016 (retained if the bit score was 50 or higher), and filtered for microbial sequence contamination using the same BLAST 63 method. Noncoding RNA was removed based on homology to Rfam sequences of Danio rerio, and highly similar sequences were collapsed using CD-HIT to form the “representative transcriptome”. This Transcriptome Shotgun Assembly project has been deposited at DDBJ/ENA/GenBank under the accession GEVC00000000. The version described in this paper is the first version, GEVC01000000.

Percent Identity Matrix: Sequences acquired from the current assembly were aligned with sequences of previously published sequences from different species to evaluate phylogenic relationships within fish. Alignments were trimmed with Gblocks 64, 65, concatenated manually into a single ordered sequence, and aligned with ClustalW 66 to produce a percent identity matrix.

Transcript Coverage against Model Organisms: A FASTA file was compiled containing all protein sequences derived from the T. areolatus transcriptome; the sequences were BLASTed against all non-redundant Salmo salar and Danio rerio protein sequences retrieved from NCBI. The BLASTp algorithm was used to establish local alignments with E-values smaller than 1e-5. The sequences demonstrating the highest-scored alignments were kept. The length of each high-scoring alignment was subsequently compared to the overall length of the reference sequence to obtain the coverage, after which a count of the unique sequences present in a particular range of coverage was obtained through the Linux command line. The results were then compiled into Figure 3.

Gene Ontology: BLASTp was used to assign the top hit for the TransDecoded T. areolatus proteins and the proteomes of Salmo salar, Cyprinus carpio, and Danio rerio (retrieved from NCBI January 29th, 2016) against the NCBI non-redundant database (GI list: Arabidopsis thaliana, Caenorhabditis elegans, Danio rerio, Dictyostelium discoideum, Drosophila melanogaster, Escherichia coli, Gallus gallus, Homo sapiens, Mus musculus, Rattus norvegicus, Saccharomyces cerevisiae, and Schizosaccharomyces pombe retrieved June 10th, 2016). The resulting file was scanned by BLAST2GO (version 2.8.0) 67 with the b2g_feb16 GO database; once terms were assigned, level 2 GO ID terms were utilized for a comparison with all groups.

KOG: Analysis included translating publically available transcripts from chosen fish species using TransDecoder version 3.0 61. These results were aligned using RPS-BLAST (E-value of ≤ 1e-5) to the NCBI KOG database (version 3.14).

KEGG: The KEGG analysis was conducted by uploading organism transcriptomic sequences to the KEGG Automatic Annotation Server (KAAS) where they are processed by BLAST and GHOST comparisons against a database of KEGG genes. Ortholog assignments were returned and graphed to illustrate a transcriptome pathway comparison.

Differential Expression Analysis: Differential transcript expression analysis of T. areolatus transcript files were aligned to the assembled reference sequences. Reads were mapped with Bowtie 68, and RSEM 69 was used to quantify differential expression values. Transcripts Per Kilobase Million (TPM) values were calculated to produce fold changes. Transcripts were then annotated by homology to NCBI Refseq proteins.

Genes of Interest: To identify genes related to endocrine and/or oxidative stress, query protein sequences from published sequences available on NCBI were obtained. Sequences from Danio rerio were used with BLASTp to find top hits in the putative coding sequence transcriptome (T. areolatus). Top hits were reviewed for the highest bit score (in combination with E-value) and included in Table 5. Protein isoforms were differentiated based off of chosen model organism annotation and could vary between query species (e.g. estrogen receptor alpha, beta).