\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[ngerman,greek,polish,english]{babel}



\usepackage{float}








\begin{document}

\title{Influence of RNA-Seq library construction, sampling methods, and tissue
harvesting time on gene expression estimation}



\author[1]{Ylenia Chiari}%
\author[2]{Leif Howard}%
\author[1]{Nickolas Moreno}%
\author[3]{Scott Relyea}%
\author[3]{James Dunnigan}%
\author[4]{Matthew Boyer}%
\author[5]{Marty Kardos}%
\author[1]{Scott Glaberman}%
\author[2]{Gordon Luikart}%
\affil[1]{George Mason University}%
\affil[2]{University of Montana Missoula Division of Biological Sciences}%
\affil[3]{Montana Fish Wildlife and Parks}%
\affil[4]{University of Montana}%
\affil[5]{National Oceanic and Atmospheric Center}%


\vspace{-1em}



  
  \date{May 22, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
RNA sequencing (RNA-Seq) is a popular method for measuring gene
expression in non-model organisms, including wild populations. While
RNA-Seq can measure gene expression variation among wild-caught
individuals and can yield important biological insights into organism
function, sampling methods may also influence gene expression estimates.
We examined the influence of multiple technical variables on estimated
gene expression in a non-model fish, the westslope cutthroat trout
(Oncorhynchus clarkii lewisi), using two RNA-Seq library types: 3'
RNA-Seq and whole mRNA-Seq. We evaluated effects of dip netting versus
electrofishing, and of harvesting tissue immediately versus 5 minutes
after euthanasia on estimated gene expression in blood, gill, and
muscle. We detected 30\% more genes with whole mRNA-Seq than with 3'
RNA-Seq and found that 58\% of genes were significantly differently
expressed between 3' RNA-Seq and whole mRNA-Seq. Our findings indicate
that 3' RNA-Seq and whole mRNA-Seq are robust to the technical variables
related to the field sampling approaches tested here with a lack of
differential gene expression among sampling methods and tissue
collection time after euthanasia. However, we found that gene expression
varied based on which RNA-Seq library type was used on the same set of
samples. Our study suggests researchers could safely rely on different
fish sampling strategies in the field and save money and analyze more
individuals using 3' RNA-Seq, but should use whole mRNA-Seq when working
with a species without good genomic resources, and when maximizing the
number of genes identified and detecting alternative splicing are
important.%
\end{abstract}%



\sloppy


\textbf{1 Introduction}

RNA sequencing (RNA-Seq) is increasingly common in ecological and
evolutionary studies focusing on variation in gene expression (Alvarez
et al., 2014, Conesa et al., 2016; Ekblom \& Galindo, 2011), It has been
used in research on physiology, conservation, and to assess organismal
response to environmental variables (Todd et al., 2016; Corlett, 2017,
Rey et al., 2020). RNA-Seq is highly accurate for quantifying expression
levels, requires less RNA sample compared to microarrays, does not
necessarily require a reference genome (e.g., Cahais et al., 2012), can
uncover sequence variation in transcribed regions, and shows high
reproducibility (Wang et al., 2009). However, gene expression data can
be strongly influenced by biological and non-biological factors such as
experimental and stochastic variation (Auer \& Doerge, 2010; Qian et
al., 2014; Todd et al., 2016). Given the recent surge in RNA-based
studies, it is therefore critical to identify and quantify
non-biological sources of variation in gene expression estimates.

Tissue sampling methods can be an important experimental cause of
variation in estimated gene expression (Mutch et al., 2008; Passow et
al., 2019). Delay in sample preservation after collection, for example
by increasing storage time in buffer at room temperature for more than
10 days, may result in higher RNA degradation and introduce bias in
estimated gene expression (e.g., Gayral et al., 2011; Romero et al.,
2014). This is a consequence of mRNAs being produced in relatively short
or rapid bursts in response to internal or external stimuli and having
short half-lives (Ross, 1995; Staton et al., 2000). Similarly, the use
of different anesthetics, methods of tissue preservation, different RNA
extraction methods, and timeframe between sample collection and RNA
isolation can all impact RNA quality and gene expression (e.g., Debey et
al., 2004; Huitink et al., 2010; Jeffries et al., 2014; Mutter et al.,
2004; Olsvik et al., 2007; Passow et al., 2019).

Variation in gene expression due to stochastic variation in cellular and
molecular processes can result in random differences among individuals
of the same population for the same genes without necessarily being a
consequence of micro-environmental variation or other biological factors
(e.g., maternal effects and potentially heritable variation). For
studies with few biological replicates, this variation may be
misinterpreted as biologically relevant (Hansen et al., 2011; Kaern et
al., 2005). Detection of stochastic variation in gene expression may be
achieved through careful sampling design (e.g., individuals vary at only
one treatment) and by increasing the number of sampled individuals (Kim
et al., 2015, Liu et al., 2014) to gain statistical power (Ching et al.,
2014). However, RNA-Seq experiments are often limited in the number of
sampled individuals due to cost, with consequent loss of statistical
power and potentially misleading results (Bi \& Liu, 2016; Li et al.,
2013).

Higher sequencing cost has led to the development of RNA library
construction protocols that allow processing and sequencing a larger
number of samples in a more cost-effective manner (Meyer et al., 2011;
Morrissy et al., 2009; Wu et al., 2010). 3' RNA-Seq methods only primes
the 3' poly-A tail, thus reducing the sequencing effort and cost (Lohman
et al., 2016; Ma et al., 2019). Independent of sample size, library
construction and RNA sequencing techniques however may also produce
variability in detection of transcripts, detection of differentially
expressed genes among treatments and observed differences in gene
expression between whole mRNA-Seq and 3' RNA-Seq (e.g., Crow et al.,
2020; Jarvis et al.; 2020; Ma et al., 2019; Tandonnet \& Torres, 2017).
Furthermore, whole mRNA libraries and sequencing methods often result in
fragment length bias because longer transcripts are sheared into more
fragments so that a higher number of reads will be assigned to them than
shorter transcripts, causing an overrepresentation of larger transcripts
(Ma et al., 2019; Oshlack \& Wakefield, 2009; Roberts et al., 2011). On
the other hand, 3' RNA-Seq generates an essentially uniform distribution
of fragments with respect to original RNA length (Lohman et al., 2016;
Ma et al., 2019). Although there are methods to correct for the bias in
gene expression due to differences in transcript length, the detection
and sampling of transcripts is still higher -- especially for longer
transcripts -- when using classical mRNA-Seq approaches (Crow et al.,
2022; Ma et al., 2019; Mandelboum et al., 2019; Tandonnet \& Torres,
2017). Finally, whole mRNA-Seq libraries permits identification of
alternative splicing at a single gene, as library and sequencing with
this method capture different fragments and transcripts for the same
locus (Crow et al., 2022).

In many species including fish, RNA-Seq data are commonly used to
investigate the effects of environmental variables (e.g., temperature,
hypoxia) on gene expression (e.g., Krishnan et al., 2020; Long et al.,
2015; Meyer et al., 2011; Smith et al., 2013; Wang et al., 2015;
Jeffries et al., 2021). However, there is little known about the
influence of different methods used to sample individuals , under field
conditions on gene expression. Field conditions may limit the use of
optimal sampling protocols or storage methods(Mutter et al., 2004;
P\selectlanguage{ngerman}érez\selectlanguage{english}-Portela \& Riesgo, 2013). Handling time of individuals before
tissue sampling may also be longer than in the lab and affect gene
expression differently depending on the field sampling technique and
tissue used.

The impacts of handling stress on fish physiology are well understood
(Sopinka et al., 2016). Although most studies focus on glucocorticoid
and blood chemistry responses following capture (Milla et al., 2010;
Wiseman et al., 2007; Wood et al., 1983; Milligan, 1996; Barton, 2002;
Ruane et al., 2001; see also Romero \& Reed, 2005 for influence on
handling time of non-fish species), gene expression responses to
handling stress indicate that the magnitude, intensity, and duration of
changes vary across genes, species, and tissue types (Krasnov et al.,
2005; Lopez et al., 2014). While there is evidence that blood cortisol
and glucose levels are affected by capture method (e.g.,
electrofishing), to our knowledge (Barton \& Dwyer, 1997; Barton \&
Grosh, 1996; Bracewell et al., 2004), it is unclear whether gene
expression is affected by capture method or handling time prior to
sample collection.

Here, we test whether sampling method (electrofishing vs dip netting),
processing time, and RNA-Seq libraries (3' RNA-Seq -- here called
QuantSeq - vs. whole mRNA-Seq -- here called NEB) influence gene
expression data in multiple tissue types from westslope cutthroat trout
(\emph{Oncorhynchus clarkii lewisi} ), a species of conservation concern
native to western North America (Behnke, 2002; Allendorf and Leary,
1988; Shepard et al., 2003). Electrofishing, which consists of a
backpack mounted electrofishing unit that applies an electrical current
in the water to momentarily stun the fish, is one of the most common
fisheries sampling methods. This method may cause the fish to express
genes in response to the electric current, and may affect individual
fish and tissue types differently, increasing variation among biological
replicates. An alternative to electrofishing is dip netting. While nets
may potentially result in a lower effect on gene expression and lower
risk of inadvertently killing both target and non-target organisms, it
is more laborious and time consuming and less effective in the field
where circumstances may not allow for long sampling periods or aquatic
systems may have obstacles that prevent effective capture with nets
(e.g., fallen tree limbs and rocks). Capturing fish by dip netting may
still influence gene expression through stress, as the fish tries to
escape capture.

The results of this study will provide a foundation for improving future
RNA-based study designs for field sampling of wild caught non-model fish
and other species.

\textbf{2 Methods}

\textbf{2.1 Sample collection, group assignment, and tissue harvesting}

All samples of westslope cutthroat trout were collected on a single day
in May 2019 at the Montana Fish, Wildlife, \& Parks Sekokini Springs
hatchery in West Glacier, MT (USA). We collected 30 fish divided in
three treatment groups (10 fish per group), each with 4 tissues for a
total of 120 samples as follows (Tables 1 and 2 and available on Dryad):
group 1 = net-sampling with immediate tissue harvest (n samples = 40);
group 2 = electrofishing with immediate tissue harvest (n samples = 40);
group 3 = electrofishing, tissue harvested from fish 5 minutes after
death by pithing (see below, n samples = 40). All fish were fry (1 year
old; fish were non-sexually mature as they cannot be sexed until 2-4
years old) from the same breeding stock and were offspring (F1) from
wild parents from Danaher Creek (MT). Average size of fish was 108 mm
+/- 11 and average weight was 10g +/- 3.

Fish were gently netted five at the time from the raceway into holding
buckets containing hatchery system water. Fish were then either captured
from the bucket by net or electrofished backpack electrofishing unit set
to 150 volts with a standard pulse for a duration of 3 seconds. Fish
were then euthanized by pithing and processed for tissue harvesting,
except for group 3. Fish from group 3 were sampled in the same way as
fish from group 2, except that after pithing they were placed in a
separate holding bucket of water for 5 minutes before tissue harvesting
to test for the influence of delayed tissue harvesting. Time between
capture and euthanasia and duration of tissue collection were recorded
for each individual. Average time in the bucket was approximately 2
minutes before euthanasia and average time of tissue collection after
pithing was approximately 3 minutes, except for group 3, for which
tissue harvesting began 2-3 min after the 5 min from pithing. Length and
weight data were collected for each fish. Sample information, including
times of tissue harvest after euthanasia for each sample can be found in
Table S1 on Dryad).

Tissue removal was performed using single use scalpels on a nylon
cutting board. Tissue samples from each fish were collected in the
following order: blood, dorsal muscle, liver, and gills. We first
collected the blood immediately before euthanasia as coagulated blood
may affect RNA quality (Chiari and Galtier, 2011). To obtain the blood
sample, the tail was removed by a diagonal cut made through the caudal
peduncle from dorsally anterior of the anal fin to ventrally posterior
of the anal fin to avoid intersecting the gastrointestinal tract. Slight
pressure was applied to the body of the fish and blood was allowed to
drip out of the cut directly into the 2 mL tube. Muscle tissue was
sliced into smaller pieces to allow penetration of the preservative
(Gayral et al., 2011). Sampling tools and the cutting board were
thoroughly cleaned with 10\% bleach first and then purified water
between fish to avoid sample and tissue contamination. Tissue samples
were placed in 2 mL sterile tubes filled with RNAlater (Qiagen) for
preservation. Tubes were left at room temperature overnight and then
stored at -80C (or in dry ice for transportation) until the RNA
extraction was carried out. All sampling was carried out according to
the Institutional Animal Care and Use Committee (IACUC) approved permit
\#AUP 007-19GLFLBS-062819 to GL.

\textbf{2.2 RNA extraction}

RNA extractions and the following laboratory procedures described below
were carried out by a private company (Admera Health). The same
extraction protocol was used for each of the different tissues and
generally followed manufacturer instructions for QIazol (Qiagen)
extraction. Briefly, up to 10 mg of tissue was mechanically homogenized
in 500 \selectlanguage{greek}µ\selectlanguage{english}l of QIazol. After homogenization, QIazol was added to reach 1ml
and then 200 \selectlanguage{greek}µ\selectlanguage{english}l of chloroform were added and mixed. For blood samples,
they were centrifuged at 2000 g for 5 minutes, the supernatant discarded
and 1ml of QIazol (Qiagen) added to the tube. Tubes were then left at
room temperature for 5 minutes and vortexed to ensure homogenization of
the sample. 200 \selectlanguage{greek}µ\selectlanguage{english}l of chloroform was added and mixed. All samples (blood
or other tissues, all containing 1ml of QIazol and 200 \selectlanguage{greek}µ\selectlanguage{english}l of chloroform)
were then incubated at room temperature for 3-5 minutes and centrifuged
at 4 \selectlanguage{ngerman}°C, 12,000 g for 15 minutes. The upper aqueous RNA containing phase
was transferred to a new tube. An equal volume of 70\% ethanol was added
and mixed. The mixture was loaded into a RNeasy mini prep column (Qiagen
RNeasy Mini Plus Kit) and RNA eluted following the manufacturer's
protocol.

The quality of RNA was determined by Qubit HS RNA assay (ThermoFisher),
and the integrity of RNA was evaluated based on RIN (RNA integrity
number, varies between 1 -- 10 with 10 when there is no degradation)
acquired via capillary gel electrophoresis performed using Bioanalyzer
2100 (Agilent Technologies). ANOVA was run in R using the F-test to
compare RIN numbers among the samples belonging to the different groups
(see below) and to compare the RIN numbers among samples belonging to
different tissues in each group. These analyses were run without data
from the liver for which most samples showed signs of degradation
(average RIN 8.0±1.21 , see also data on Dryad ). All the results
presented in the Result section therefore do not include data on the
liver samples.

\textbf{2.4 RNA Library preparation and Sequencing}

Since variation in RNA quality may affect downstream results (Passow et
al. 2019), library construction and sequencing were carried out for 81
tissue samples with a RIN value above 8.8 for QuantSeq and a subset of
14 blood samples (for which we also had QuantSeq data) with
RIN\emph{\textgreater{}} 9.4 for NEB (Tables 1 and 2 and data on Dryad).
None of these samples showed signs of RNA degradation based on the
BioAnalyzer profile(Tables 1 and 2 and data on Dryad ). Whole mRNA
libraries (NEB) were made only for 14 selected blood samples with
similar RIN and concentration among compared groups (Tables 1 and 2).
Library preparation was performed with the NEB Ultra II RNA library prep
kit with NEBNext Poly(A) mRNA Magnetic Isolation Module (New England
Biolabs). For 3'-end RNA Tag-Seq, library preparation was performed with
QuantSeq 3' mRNA-Seq Library Preparation Kit FWD for Illumina (Lexogen).
All procedures were performed according to manufacturer suggested
protocols. The quantity and molecular size of the libraries were
confirmed by Qubit HS DNA assay (ThermoFisher) and Tapestation 2200
system coupled with High Sensitivity D1000 ScreenTapes (Agilent).
Sequencing was performed on Illumina Hiseq X with 150bp pair-end reading
for all the QuantSeq samples (Lexogen) and four NEB samples, while the
remaining 10 NEB samples were sequenced on a NovaSeq machine (see
Results section regarding lack of difference between the NEB samples
sequenced with different machines). Raw reads were deposited on NCBI
(SRA PRJNA691889, available after manuscript acceptance).

\textbf{2.5 RNASeq reads check and genome coverage}

Quality checks of the raw RNA-Seq reads were performed using Fastqc
(Andrews, 2014).

Reads were trimmed with trimmomatic using the default parameters
(version 0.38, Bolger et al. 2014). Raw reads were mapped to
an\emph{Oncorhynchus mykiss} reference genome from NCBI (Omyk\_1.0,
https://www.ncbi.nlm.nih.gov/assembly/GCF\_002163495.1/, Annotation
release ID:100) using STAR (version 2.7.1a; Dobin et al., 2013; Dobin
and Gingeras, 2015) to obtain the number of genes detected by each
technique, QuantSeq vs. NEB.

In order to perform bioinformatic analyses on samples with an equal
number of uniquely mapped reads among the two library types (as also
previously done by others, see Ma et al., 2019), we randomly selected 11
million and 40 million reads per sample for all analyses performed on
QuantSeq and NEB, respectively. Previous work has shown that using
\textasciitilde{} \textgreater{}10M reads does not increase the amount
of uniquely mapped reads, after which the ability to detect differently
expressed genes becomes independent of sequencing depth (Ramsköld, 2009;
Liu et al. 2014; Ma et al., 2019; see also Crow et al., 2022 for an in
depth discussion on the issue of redundancy of reads). Transcripts are
randomly sheared into fragments with NEB but not QuantSeq. Consequently,
the number of reads with NEB are proportional to the number of fragments
not transcripts, whereas the number of reads with QuantSeq is
proportional to the number of transcripts. Because of this, more reads
may be needed for NEB than for QuantSeq to have a similar percentage of
uniquely mapped reads. However, having many fragments corresponding to
the same transcript is redundant and not useful for gene expression
quantification (Crow et al., 2022). At the same time, especially for
longer transcripts, the increased number of fragments per transcript
obtained with NEB library increases transcript detection (Crow et al.,
2022; Ma et al., 2019). The issue for whole mRNA-Seq of over-counting
the same transcript due to multiple fragments corresponding to the same
locus is especially relevant at higher sequencing depths (see Crow et
al., 2022 for a detailed discussion on this issue), so that randomly
selecting a lower number of reads to obtain an equal number of uniquely
mapped reads among the two library types overcome this issue (see also
Crow et al., 2022).

Reads were mapped again to the \emph{Oncorhynchus mykiss} reference
genome. HT-Seq (version 0.11.1; Anders et al. 2015) was then used to
quantify the number of reads uniquely mapped to each gene of the\emph{O.
mykiss} reference genome. Finally, a python script provided with
Stringtie (prepDE.py) was used to generate a gene counts matrix (Pertea
et al., 2016).

\textbf{2.6 Similarity in gene expression among samples}

To assess the variation and direction of variation among samples based
on their gene expression, we calculated the correlation of gene
expression levels among samples and the Euclidean distances among
samples in DESeq 2 (version 1.22.2; Love et al., 2014) following the
program directions. These measures are especially useful to assess the
similarity of biological replicates (e.g., samples belonging to the same
group) (Koch et al. 2018) and therefore to detect anomalies among the
samples. The sample correlation matrix was calculated by performing the
Pearson correlation of the normalized matrix after the variance
stabilizing transformation (\emph{vst} ) was performed on the most
variable 2000 genes based on the HTSeq data produced. \emph{vst} allows
taking into account the sample variability of low counts.

Sample Pearson correlation is calculated in pairwise comparison between
samples and ranges from -1 to 1, where a value of 0 indicates no
correlation (gene expression is completely dissimilar between the two
samples), while values of 1 indicate that the samples have identical
expression level (and -1 corresponds to negative correlation). The
Euclidean distance between sample expression profiles was calculated by
this equation: \emph{dist = sqrt(1- cor\textsuperscript{2})} ,
where\emph{cor} stands for the correlation coefficient of 2 samples. The
smaller the distance, the higher the correlation between samples. These
distances were then used to build the heatmaps of sample distance of
each normalized matrix, which allows the data to be shrunken towards the
genes' average expression across all samples. Gene heatmaps were instead
based on vst transformation to normalize the raw count. After this, the
mean expression in each sample is then normalized to 0. Finally,
differences in gene expression among the studied groups (see below) were
visualized by a PCA plot using the gene count matrix after applying the
variance stabilizing transformation (\emph{vst} ) to normalize the raw
counts. PCA plots are useful to assess the effect of covariates and
batch effects (non-biological variation due to experimental artifacts
(Reese et al. 2013).

\textbf{2.7 Differential gene expression analysis}

Differences in gene expression among groups were identified by
differential expression analysis performed using DESeq2 on raw read
counts (non-normalized, as suggested by DESeq2). To allow comparison
between QuantSeq reads and NEB reads data, we only used one sequencing
direction as suggested for the DESeq2 program for the latter. The false
discovery rate (FDR) was adjusted to 0.05, corresponding to a recovery
at most of 5\% of false positives following the DESeq2 manual. We use
the default options for all the other parameters. We look at differences
in gene expression between sampling methods, harvest tissue time, tissue
type, and QuantSeq vs. NEB in Table 2: (see Tables 1, 2 and data on
Dryad per detailed information about comparisons and sample size for
each comparison; minimum N = 4). The log2 fold changes obtained from
DESeq2 were used as a measure of how many more (or less) genes are
expressed in one group versus the other. We considered genes having
different expression if the adjusted p-value (using the adjusted p-value
results in less false positives) was \textless{} 0.05.

Finally, previous work has indicated an increase in read count for
longer transcripts using NEB than QuantSeq (Ma et al. 2019; but see Crow
et al., 2022). To further address the relationship between gene length
and genes differentially expressed between QuantSeq and NEB, we
conducted an assessment using the known transcript length from
orthologous genes in zebrafish in Ensembl 101 (Yates et al. 2020) based
on gene name for genes that were detected to be differentially expressed
between the two library types. We used the zebrafish instead of the
rainbow trout genome as the former has a more curated annotation (and
thus more precise gene length information) than the latter. We also used
the same approach to specifically assess if transcript length could
influence absence of gene expression or not detection (mean base = 0 in
DESeq2 output) in one but not the other library type, QuantSeq or NEB.

\textbf{3. Results}

\textbf{3.1 RNA and raw sequencing data quality statistics}

Out of the 120 samples for which RNA was extracted, 86 had a RIN value
(a measure of RNA integrity) equal or above 8.8. Little variation in RIN
scores was observed among the sampled tissues and sampling methods
(Table 1 and data on Dryad). Mean and standard deviation for RIN values
for the three tissues were: 9.6±0.22 (blood), 9.2±0.40 (muscle), 9.0±1
(gill). Mean and standard deviation for RIN values for the three
treatment groups: 9.2±0.43 (dip netting), 9.3±0.34 (electrofishing), and
9.2±1.06 (tissue harvesting after 5 minutes). We found no differences in
RIN values among groups (F = 0.299, \emph{df} = 2, \emph{p} = 0.74) and
in RIN values among tissues within each group (F = 0.595, df =
4,\emph{p} = 0.67).

Table S1 on Dryad ). The final number of reads per individual for
QuantSeq libraries ranged from 11 million to 15.6 million (mean = 12.88
million ± 0.67; data on Dryad). On average, of the 11 million reads
randomly selected for each sample, we obtained around 77\% of uniquely
mapped reads on the rainbow trout (\emph{O. mykiss)} genome
independently of the sampling method used (range: 67.7 \selectlanguage{english}- 86.3\%, Table
S1 on Dryad .

RNA sequencing from the 14 NEB samples (blood only) yielded a total of
564 million reads for individuals captured by net (mean = 112.9
million\emph{+} 13.95; N = 5), 563.4 million reads for samples collected
by electrofishing and sampled immediately (mean = 112.7 million
\emph{+}22.4; N = 5), and 350.4 million reads from electrofishing
samples processed after 5 minutes (mean = 87.6 million \emph{+} 7.4; N =
4). The final number of reads per individual ranged from 77.8 to 148.8
million reads (mean = 105.6 million +- 19.1). On average, of the 40
million reads randomly selected for each sample, we obtained on average
75\% of uniquely mapped reads on the \emph{O. mykiss} genome (Table S1
on Dryad).

Mapping reads on the rainbow trout (\emph{O. mykiss)} genome and looking
at the output of DESeq2 between the two library types and for the same
samples, we found that NEB detected 30\% more genes than QuantSeq (Table
S2 on Dryad). A gene was considered to be detected/expressed when
basemean was different from 0. Specifically, for the 14 blood samples
for which the two different types of RNA-Seq libraries were built, we
found that NEB and QuantSeq detected approximately 35K and 26K genes,
respectively, which mapped onto the annotated rainbow trout (\emph{O.
mykiss} ) genome. 25K genes were detected by both library types.
However, 9K genes were detected by NEB but not QuantSeq, and 1K genes
were detected by QuantSeq but not NEB. Presence/absence of genes
detected by one or the other library type is independent of gene
transcript length (Figure 1).

\textbf{3.2 Similarity in gene expression among samples}

Similarity in gene expression among biological replicates - i.e.,
individuals belonging to the same treatment group - gives an idea of
reproducibility of our data and of the overall variation among samples.
Similarity in gene expression within and among groups can be estimated
using the sample correlation or Euclidean distances (see Materials and
Methods for further details). Pearson correlation coefficients (r) for
biological replicates were equal or above 0.9 for 97\% of comparisons
(same tissue and some tissue within a group) (Table S3 on Dryad). This
indicates that although variation in gene expression occurs among
individuals, biological replicates are generally very similar.

Pearson r values between the two sequencing platforms for NEB are all
above 0.9 for the samples belonging to the same group (Supporting
Information Table S2), indicating that different sequencing methods did
not influence the number of uniquely mapped reads. Finally, r among
different tissues (for QuantSeq) and among QuantSeq vs. NEB are
generally \textless{}0.5 and sometimes negative, suggesting different
levels of gene expression among tissues and among the same mapped genes
between the two library types.

Heatmaps of the distance matrices for the different group comparisons
provide hierarchical clustering based on sample distances. When heatmaps
were made combining data from the three different tissues for QuantSeq,
we found three clusters corresponding to the three different tissues
(Figure 2A). However, within each cluster, as also shown by the heatmaps
built with data from each tissue separately, samples belonging to
different groups are clustered together, indicating no clear difference
in gene expression among the tested groups (Supporting Information
Figure S1). Lack of difference in gene expression among the different
groups was also found using NEB data (Figure 2).

Finally, comparison of QuantSeq vs. NEB found differences in gene
expression between the two methods; however, this difference was not
associated with any of the groups (Figure 2). Principal component
analysis (PCA), another way to visualize variation in gene expression
among samples, further supports the lack of differences among sampling
methods and time of tissue harvesting and the differentiation between
QuantSeq versus NEB and among the three sampled tissues (Figures 3 and 4
and Supporting Information Figure S2).

\textbf{3.3 Differential gene expression}

\textbf{3.3.1 Dip Netting versus electrofishing sampling method}

QuantSeq data identified only 3 out of 39,212 genes (0.008\%) that were
significantly (adjusted \emph{p-} value \textless{} 0.05) differentially
expressed between dip netting and electrofishing across all tissue
types. When gene expression between dip netting and electrofishing was
analyzed separately for different tissue types, 0 out of 27,118 genes
were differentially expressed in blood, 16 out of 20,465 genes (0.08\%)
in gills, and 155 out of 9,201 genes (1.7\%) in muscle. NEB data
(available for blood samples only) confirmed QuantSeq data, with no
genes showing differential expression between the two sampling methods
(0 out of 33,236 genes) (Table S2 on Dryad).

\textbf{3.3.2 Immediate versus delayed postmortem tissue harvesting}

We found no significant difference (adjusted \emph{p-} value \textless{}
0.05) in gene expression between samples for which tissues were
harvested immediately versus \textasciitilde{}5 minutes after
euthanasia. QuantSeq data identified only 1 out of 38,864 genes
(0.003\%) that were significantly differentially expressed between
tissue harvesting times across all tissue types. Similar to what was
observed between sampling techniques, when tissues were analyzed
separately for harvesting times, 0 out of 27,401 genes were
differentially expressed in blood, 18 out of 15,310 genes (0.1\%) in
gills, and 3 out of 34,460 genes (0.009\%) in muscle (Table S2 on
Dryad). NEB data (available for blood samples only) was similar to
QuantSeq data, with only 6 out of 23,265 genes (0.03\%) showing
differential expression between the two sampling methods (Table S2 on
Dryad).

\textbf{3.3.3 QuantSeq versus whole mRNA-Seq}

We compared the gene expression of the same 14 blood samples (Ntot=28)
among RNA-Seq libraries built using QuantSeq and NEB. We found 15,328
out of 26,316 genes (58\%) that were significantly (adjusted\emph{p}
-value \textless{}0.05) differentially expressed among library types
(Table S2 on Dryad). Specifically, we found 9852 (64\%) transcripts with
higher expression in NEB versus QuantSeq and 5476 (36\%) transcripts
with higher expression in QuantSeq versus NEB. Within sampling methods,
the percentage of differentially expressed genes between QuantSeq and
NEB was 45\%, 41\%, and 44\% for dip netting, electrofishing sampled
immediately, and electrofishing sampled after 5 min, respectively (Table
S2 on Dryad). The proportion of identical genes with differential
expression favoring one of the two library methods is very similar,
although does favor NEB with increased transcript length (Figure 5). Log
of the basemean expression differences for NEB relative to QuantSeq was
only 0.098 for transcripts \textless{}1000 bp, but was 0.27 for
transcripts between 5000-6000 bp (Figure 5).

\textbf{3.3.4 Comparisons between tissue type within the same sampling
groups}

We found that blood and muscle have the lowest number of genes with
detected expression in both tissues compared to the other pairwise
tissue comparisons (Tables S2 and S4 on Dryad). We observed differences
in gene expression among tissue types (blood, muscle, gills) for each of
the sampling methods. Of all the genes that are significantly
differentially expressed (adjusted \emph{p} - value \textless{}0.05)
among tissue types, only 25-32\% are more highly expressed in blood
compared to gill and muscle (Tables S2 and S4 on Dryad). Conversely, the
proportion of significantly differentially expressed genes in gill and
muscle was much more even and within +-10 of 50\% (Tables S2 and S4 on
Dryad).

\textbf{4 Discussion}

The increasing use of RNA-Seq for ecological, physiological, and
evolutionary studies on wild caught organisms has required appraisal of
the influence of different sampling techniques, storage methods,
processing time, and tissue types on RNA quality and data production
(Camacho-Sanchez et al., 2013; Cheviron et al., 2011; Nakatsuji et al.,
2019; Yu et al., 2013). Among the most important applications of RNA-Seq
currently is testing for rapid adaptation to environmental change (e.g.,
Narum \& Campbell, 2015; Connon et al., 2018) and for its inheritance
(e.g., Christie et al. 2016, Charlesworth et al. 2017, Skvortsova et al.
2018, Navarro-Martin et al. 2020, Savilammi et al. 2020), and for
addressing questions in evo-devo (e.g., Roux et al., 2015; Liu et al.,
2020).

In our work, we tested if different sampling techniques influenced gene
expression across different tissues in westslope cutthroat trout.
Overall, we obtained high RNA quality for all tissues (mean
RIN\emph{\textgreater{}} 9.0 for the different tissues) except liver
(mean RIN = 8.0). Liver has a high rate of protein synthesis and
degradation, and the higher RNA degradation observed for this tissue in
comparison to blood, muscle, and gills is likely the result of higher
enzymatic activity in the liver (Carter et al. 2001, Wiseman et al.
2007). In our experiment, liver was the third tissue sampled after
euthanasia, after blood and muscle, and it took us between 2 and 3
minutes to sample. Because of its importance in detoxification
mechanisms, physiological studies may require target sampling of this
tissue. We therefore suggest sampling of liver first if more than one
tissue is sampled to minimize RNA degradation.

We also found no difference in RNA quality among samples obtained
through dip netting or electrofishing even when tissue was not harvested
until 5 minutes after euthanasia. While opinions on a cutoff threshold
RIN value to obtain reliable gene expression data differ, it has been
shown that partially degraded RNA may still detect the same uniquely
mapped genes as non-degraded RNA, although the coverage of mapped reads
is lower for partially degraded RNA and gene specific (Romero et al.
2014, Wang et al. 2016). However, while RNA degradation may not strongly
affect mapping, it may drastically affect estimates of differential gene
expression (Chen et al. 2014, Romero et al. 2014). Furthermore,
different RNA-Seq techniques may be differentially affected by RNA
degradation (Adiconis et al. 2013), requiring selecting the most
appropriate RNA-Seq library depending on RNA quality (Adiconis et al.
2013).

We found that gene expression among individuals belonging to the same
group and tissue type were very similar for the majority of comparisons
(correlation coefficients \emph{\textgreater{}} 0.9), independent of the
sampling method or harvesting time. However, we observed among-sample
variation in gene expression, reflecting the importance of larger sample
size in RNA-Seq studies to decrease the influence of stochastic effects
on variation in gene expression that could otherwise be interpreted as
biologically relevant (Ching et al. 2020). Furthermore, we also observed
similarity of expression levels among samples obtained with the two
sampling methods, dip netting or electrofishing, or subjected to
different tissue harvest times (immediate or 5 minutes after
euthanasia). Sampling individuals of the same age, in the same
environment and on the same day, with many biological replicates per
treatment and using only samples with highly similar RNA quality most
likely reduced the effects of non-biological variation and of
non-relevant biological variation in our experiments (Fang \& Cui 2010,
Wong et al 2012, Yu et al. 2014).

We detected a higher number of mapped and expressed genes
(\textasciitilde{}30\% higher) for samples processed with NEB than with
QuantSeq, independent of gene transcript length. Others have proposed
traditional whole mRNA to detect more genes than 3' RNA-Seq libraries
(e.g., Crow et al., 2022; Ma et al., 2019; Xiong et al., 2017).
Furthermore, we observed that while NEB detected 30\% more expressed
genes than QuantSeq, QuantSeq also detected a smaller number of
expressed genes that were not found with NEB. Finally, for the same
samples processed with both library types, we found different gene
expression between NEB and QuantSeq, with a higher proportion of genes
(58\%) with greater gene expression for NEB relative to QuantSeq. As we
did not find differentially expressed genes between sampling techniques
and processing time after euthanasia (see below), we could not estimate
if NEB and QuantSeq would detect a different number of differentially
expressed genes, as previously reported by others (e.g., Crow et al.,
2022; Ma et al., 2019; Tandonnet and Torres, 2017). Different detection
of genes and differentially expressed genes between the two library
types has been proposed to depend on the length of the transcript and
how accurate and complete the annotation of the genome of the organism
is. In general, QuantSeq seems to perform better in detecting shorter
transcripts and whole mRNA-Seq in detecting longer ones. Furthermore,
since QuantSeq library data rely on mapping the reads to the 3' UTR on
the species' genome to detect genes, and since UTR regions are generally
more variable in the genome than protein coding regions, better
annotated and complete genomes facilitate mapping and detection of
transcripts/genes. In this study, mapping was carried out on a closely
related salmonid species, since the genome of the westslope cutthroat
trout is currently not available, and this can explain why many more
expressed genes were detected with NEB than QuantSeq. The higher number
of detected expressed genes suggests that researchers should use whole
mRNA-Seq for work on species with limited genomic resources (Crow et
al., 2022). Furthermore, as QuantSeq libraries only allow amplification
of the 3' end of the transcript, different transcripts resulting from
alternative cleavage sites and splicing would only be detected if the 3'
UTR were different. Traditional whole mRNA-Seq should therefore be
preferred if identifying distinct spliced transcripts may be of interest
for the study question. Finally, others (Crow et al., 2022) have
reported how the increased sequencing depth of traditional whole
mRNA-Seq methods may produce redundancy of reads without increasing the
power of detection of differentially expressed genes. As one of the
major advantages that have been reported for QuantSeq (and 3' RNA-Seq in
general) is the reduced sequencing cost, a reduced sequencing depth for
whole mRNA-Seq library can produce the desired data at a reduced cost.

One of the goals of this study was to test if different sampling methods
and processing time would affect gene expression. Although stress levels
associated with dip netting and electrofishing may differ, we found that
sampling technique did not affect gene expression levels. This result
was independent of the RNA-Seq library type (QuantSeq or NEB) and tissue
used. Although whole mRNA-Seq has been reported to be more sensitive to
differentially expressed genes than 3' RNA-Seq methods (Ma et al. 2019),
independent of the RNA-Seq library used, we found no difference in
estimated gene expression between the two field collection methods. As
field conditions often change among sampling locations, researchers
could opt to use electrofishing, where more efficient, and compare with
fish obtained by netting in other localities without introducing
extraneous variation in gene expression.

We also found that harvesting the tissue immediately or 5 minutes after
euthanasia did not produce variation in gene expression, suggesting that
it is safe to euthanize fish in batches and then proceed to tissue
harvesting. In our work, the maximum processing time of the last tissue
harvested after euthanasia was approximately 10 min (for fish processed
5 minutes after euthanasia). Although sampling techniques and tissue
processing time did not influence variation in gene expression, we
observed a large proportion of differentially expressed genes among the
different tissues.

We found fewer expressed genes in blood compared to gill and muscle, and
a smaller proportion of genes with higher expression in blood than in
the other two tissues. Blood and muscle were also the tissues with the
fewest expressed genes in common. Gill tissue had the highest number of
detected expressed genes. This may be due to active cellular processes
occurring in gills further amplified in actively growing animals
(Stolper et al. 2019). Depending on the study question, sampling
different tissues may ensure that multiple genes and multiple biological
processes are considered for studies on differential gene expression.

In summary, our study indicates that differential gene expression
results are likely to be comparable for dip netting and electrofishing.
Additionally, gill, blood, and muscle all produce good quality RNA with
reliable results if sampled within 5-10 minutes from euthanasia. Only
liver samples showed RNA degradation. Finally, although NEB library
detected more expressed genes, this did not lead to different results in
terms of distinct gene expression among the groups tested here. If
detecting alternative splicing is not of interest for the study question
and if working with an organism with good genomic resources (available
genome or a good genome available for a closely related species),
researchers can rely on using either of the library types tested here,
QuantSeq or NEB and sequence them at similar depth, reducing the cost of
NEB library sequencing. However, when it is crucial to detect as many
genes as possible and when working with an organism lacking good genomic
resources, whole RNA-Seq is recommended. These findings advance the use
of RNA-Seq to investigate gene expression variation and its role in
phenomena such as adaptation to environmental variation and climate
change in natural populations

\textbf{Acknowledgments}

MCB, SR, and JD were supported by Bonneville Power Administration grant
199101903 to Montana Fish, Wildlife \& Parks. We are very thankful to
three anonymous reviewers who provided valuable comments and improved
the quality of this manuscript.

\textbf{References:}

Adiconis, X., Borges-Rivera, D., Satija, R., DeLuca, D. S., Busby, M.
A., Berlin, A. M., Sivachenko, A., Thompson, D. A., Wysoker, A.,
Fennell, T. and Gnirke, A. (2013). Comparative analysis of RNA
sequencing methods for degraded or low-input samples. \emph{Nature
Methods} , \emph{10} (7), 623.

Allendorf, F. W., \& Leary, R. F. (1988). Conservation and distribution
of genetic variation in a polytypic species, the cutthroat
trout.\emph{Conservation Biology} , \emph{2} (2), 170-184.

Alvarez, M., Schrey, A. W., \& Richards, C. L. (2015). Ten years of
transcriptomics in wild populations: what have we learned about their
ecology and evolution? \emph{Molecular ecology} , \emph{24} (4),
710-725.

Anders, S., Pyl, P. T., \& Huber, W. (2015). HTSeq---a Python framework
to work with high-throughput sequencing data. \emph{Bioinformatics}
,\emph{31} (2), 166-169.

Andrews, S. (2014). FastQC a quality control tool for high throughput
sequence data. Retrieved from
https://www.bioinformatics.babraham.ac.uk/projects/fastqc/

Auer, P. L., \& Doerge, R. W. (2010). Statistical design and analysis of
RNA sequencing data. \emph{Genetics} , \emph{185} (2), 405-416.

Bahar, B., Monahan, F. J., Moloney, A. P., Schmidt, O., MacHugh, D. E.,
\& Sweeney, T. (2007). Long-term stability of RNA in post-mortem bovine
skeletal muscle, liver and subcutaneous adipose tissues. \emph{BMC
Molecular Biology} , 8(1), 108.

Barton, B. A. (2002). Stress in fishes: a diversity of responses with
particular reference to changes in circulating
corticosteroids.\emph{Integrative and Comparative Biology} , \emph{42}
(3), 517-525.

Barton, B. A., \& Dwyer, W. P. (1997). Physiological stress effects of
continuous\selectlanguage{english}-and pulsed\selectlanguage{english}-DC electroshock on juvenile bull
trout.\emph{Journal of Fish Biology} , \emph{51} (5), 998-1008.

Barton, B. A., \& Grosh, R. S. (1996). Effect of AC electroshock on
blood features in juvenile rainbow trout. \emph{Journal of Fish Biology}
, \emph{49} (6), 1330-1333.

Bracewell, P., Cowx, I. G., \& Uglow, R. F. (2004). Effects of handling
and electrofishing on plasma glucose and whole blood lactate
of\emph{Leuciscus cephalus} . \emph{Journal of Fish Biology} , \emph{64}
(1), 65-71.

Bi, R., \& Liu, P. (2016). Sample size calculation while controlling
false discovery rate for differential expression analysis with
RNA-sequencing experiments. \emph{BMC Bioinformatics} , \emph{17} (1),
146.

Bolger, A. M., Lohse, M., \& Usadel, B. (2014). Trimmomatic: a flexible
trimmer for Illumina sequence data. \emph{Bioinformatics} ,\emph{30}
(15), 2114-2120.

Cahais, V., Gayral, P., Tsagkogeorga, G., Melo\selectlanguage{english}-Ferreira, J.,
Ballenghien, M., Weinert, L., Chiari, Y., Belkhir, K., Ranwez, V. and
Galtier, N. (2012). Reference\selectlanguage{english}-free transcriptome assembly in non\selectlanguage{english}-model
animals from next\selectlanguage{english}-generation sequencing data. \emph{Molecular Ecology
Resources} , \emph{12} (5), 834-845.

Camacho\selectlanguage{english}-Sanchez, M., Burraco, P., Gomez\selectlanguage{english}-Mestre, I., \& Leonard, J. A.
(2013). Preservation of RNA and DNA from mammal samples under field
conditions. \emph{Molecular Ecology Resources} , \emph{13} (4), 663-673.

Carter, C., Houlihan, D., Kiessling, A., Medale, F., \& Jobling, M.
(2001). Physiological effects of feeding. \emph{Food Intake in Fish} ,
297-331.

Chen, E. A., Souaiaia, T., Herstein, J. S., Evgrafov, O. V., Spitsyna,
V. N., Rebolini, D. F., \& Knowles, J. A. (2014). Effect of RNA
integrity on uniquely mapped reads in RNA-Seq. \emph{BMC Research Notes}
, \emph{7} (1), 753.

Cheviron, Z. A., Carling, M. D., \& Brumfield, R. T. (2011). Effects of
postmortem interval and preservation method on RNA isolated from
field-preserved avian tissues. \emph{The Condor} , \emph{113} (3),
483-489.

Chiari, Y., \& Galtier, N. (2011). RNA extraction from sauropsids blood:
evaluation and improvement of methods. \emph{Amphibia-Reptilia}
,\emph{32} (1), 136-139.

Ching, T., Huang, S., \& Garmire, L. X. (2014). Power analysis and
sample size estimation for RNA-Seq differential expression. \emph{Rna}
,\emph{20} (11), 1684-1696.

Christie, M. R., Marine, M. L., Fox, S. E., French, R. A., \& Blouin, M.
S. (2016) A single generation of domestication heritably alters the
expression of hundreds of genes. \emph{Nature Communications} , 7,
10676.

Conesa, A., Madrigal, P., Tarazona, S., Gomez-Cabrero, D., Cervera, A.,
McPherson, A., Szcze\selectlanguage{polish}ś\selectlanguage{english}niak, M. W., Gaffney, D. J., Elo, L. L., Zhang, X.
\& Mortazavi, A. (2016). A survey of best practices for RNA-seq data
analysis. \emph{Genome Biology} , \emph{17} (1), 13.

Connon, R. E., Jeffries, K. M., Komoroske, L. M., Todgham, A. E., \&
Fangue, N. A. (2018). The utility of transcriptomics in fish
conservation. \emph{Journal of Experimental Biology} , 221, jeb148833

Corlett, R. T. (2017) A bigger toolbox: biotechnology in biodiversity
conservation. \emph{Trends in Biotechnology} 35: 55 -- 65.

Crow, T.M., Gill, J. A., Whitehead, A. \& Runcie, D.E. (2022) 3' RNA
sequencing does not increase power or reduce costs for gene expression
analysis. bioRxiv 2022.04.13.488255

Debey, S., Schoenbeck, U., Hellmich, M., Gathof, B. S., Pillai, R.,
Zander, T., \& Schultze, J. L. (2004). Comparison of different isolation
techniques prior gene expression profiling of blood derived cells:
impact on physiological responses, on overall expression and the role of
different cell types. \emph{The Pharmacogenomics Journal} , \emph{4}
(3), 193-207.

Dobin, A., \& Gingeras, T. R. (2015). Mapping RNA-Seq reads with
STAR.\emph{Curr. Protoc. Bioinform.} 51: 11.14.1 -- 11.14.19

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha,
S., Batut, P., Chaisson, M. and Gingeras, T. R. (2013). STAR: ultrafast
universal RNA-seq aligner. \emph{Bioinformatics} , \emph{29} (1), 15-21.

Ekblom, R., \& Galindo, J. (2011). Applications of next generation
sequencing in molecular ecology of non-model organisms. \emph{Heredity}
,\emph{107} (1), 1-15.

Fang, Z., \& Cui, X. (2011). Design and validation issues in RNA-seq
experiments. \emph{Briefings in Bioinformatics} , \emph{12} (3),
280-287.

Gayral, P., Weinert, L., Chiari, Y., Tsagkogeorga, G., Ballenghien, M.,
\& Galtier, N. (2011). Next-generation sequencing of transcriptomes: a
guide to RNA isolation in nonmodel animals. \emph{Molecular Ecology
Resources} , \emph{11} (4), 650-661.

Hansen, K. D., Wu, Z., Irizarry, R. A., \& Leek, J. T. (2011).
Sequencing technology does not eliminate biological
variability.\emph{Nature Biotechnology} , \emph{29} (7), 572-573.

Hitt, N. P., Frissell, C. A., Muhlfeld, C. C., \& Allendorf, F. W.
(2003). Spread of hybridization between native westslope cutthroat
trout, \emph{Oncorhynchus clarkii lewisi} , and nonnative rainbow
trout,\emph{Oncorhynchus mykiss.} \emph{Canadian Journal of Fisheries
and Aquatic Sciences} , \emph{60} (12), 1440-1451.

Huitink, J. M., Heimerikxs, M., Nieuwland, M., Loer, S. A., Brugman, W.,
Velds, A., Sie, D. and Kerkhoven, R. M. (2010). Volatile anesthetics
modulate gene expression in breast and brain tumor
cells.\emph{Anesthesia \& Analgesia} , \emph{111} (6), 1411-1415.

Jarvis, S., Birsa, N., Secrier, M., Fratta, P., \& Plagnol, V. (2020) A
comparison of low read depth QuantSeq 3' sequencing to total RNA-Seq in
FUS mutant mice. \emph{Frontiers in Genetics,} 11, 562445

Jeffries, M. K. S., Kiss, A. J., Smith, A. W., \& Oris, J. T. (2014). A
comparison of commercially-available automated and manual extraction
kits for the isolation of total RNA from small tissue samples. \emph{BMC
Biotechnology} , \emph{14} (1), 94.

Jeffries, M. K. S., Teffer, A., Michaleski, S., Bernier, N. J., Heath,
D. D., \& Miller, K. M. (2021). The use of non-lethal sampling for
transcriptomics to assess the physiological status of wild
fishes.\emph{Comparative Biochemistry and Physiology Part B:
Biochemistry and Molecular Biology} , 256, 110629

Kaern, M., Elston, T. C., Blake, W. J., \& Collins, J. J. (2005).
Stochasticity in gene expression: from theories to
phenotypes.\emph{Nature Reviews Genetics} , \emph{6} (6), 451-464.

Kim, J. K., Kolodziejczyk, A. A., Ilicic, T., Teichmann, S. A., \&
Marioni, J. C. (2015). Characterizing noise structure in single-cell
RNA-seq distinguishes genuine from technical stochastic allelic
expression. \emph{Nature Communications} , \emph{6} (1), 1-9.

Koch, C. M., Chiu, S. F., Akbarpour, M., Bharat, A., Ridge, K. M.,
Bartom, E. T., \& Winter, D. R. (2018). A beginner's guide to analysis
of RNA sequencing data. \emph{American Journal of Respiratory Cell and
Molecular Biology} , \emph{59} (2), 145-157.

Kovach, R. P., Al-Chokhachy, R., \& Stephens, T. (2018). Proactive
rainbow trout suppression reduces threat of hybridization in the Upper
Snake River Basin. \emph{North American Journal of Fisheries Management}
, \emph{38} (4), 811-819.

Krasnov, A., Koskinen, H., Pehkonen, P., Rexroad, C. E., Afanasyev, S.,
\& M\selectlanguage{ngerman}ölsä, H. (2005). Gene expression in the brain and kidney of rainbow
trout in response to handling stress. \emph{BMC Genomics} , \emph{6}
(1), 3.

Krishnan, J., Persons, J. L., Peuß, R., Hassan, H., Kenzior, A., Xiong,
S., Olsen, L., Maldonado, E., Kowalko, J. E. and Rohner, N. (2020).
Comparative transcriptome analysis of wild and lab populations
of\emph{Astyanax mexicanus} uncovers differential effects of environment
and morphotype on gene expression. \emph{Journal of Experimental Zoology
Part B: Molecular and Developmental Evolution} .

Li, C. I., Su, P. F., \& Shyr, Y. (2013). Sample size calculation based
on exact test for assessing differential expression analysis in RNA-seq
data. \emph{BMC Bioinformatics} , \emph{14} (1), 357.

Liao, Y., Wang, J., Jaehnig, E. J., Shi, Z., \& Zhang, B. (2019).
WebGestalt 2019: gene set analysis toolkit with revamped UIs and
APIs.\emph{Nucleic Acids Research} , \emph{47} (W1), W199-W205.

Liu, J., Frochaux, M., Gardeux, V., Deplancke, B., \& Robinson-Rechavi,
M. (2020) Inter-embryo gene expression variability recapitulates the
hourglass pattern of evo-devo. \emph{BMC Biology} , 18, 129.

Liu, Y., Zhou, J., \& White, K. P. (2014). RNA-seq differential
expression studies: more sequence or more
replication?\emph{Bioinformatics} , \emph{30} (3), 301-304.

Lohman, B. K., Weber, J. N., \& Bolnick, D. I. (2016) Evaluation of
TagSeq, a reliable low-cost alternative for RNAseq. \emph{Molecular
Ecology Resources} , 16: 1315 -- 1321.

Long, Y., Yan, J., Song, G., Li, X., Li, X., Li, Q., \& Cui, Z. (2015).
Transcriptional events co-regulated by hypoxia and cold stresses in
Zebrafish larvae. \emph{BMC Genomics} , \emph{16} (1), 385.

Lopez-Patino, M. A., Hernandez-Perez, J., Gesto, M., Libran-Perez, M.,
Miguez, J. M., \& Soengas, J. L. (2014). Short-term time course of liver
metabolic response to acute handling stress in rainbow
trout,\emph{Oncorhynchus mykiss} . \emph{Comparative Biochemistry and
Physiology Part A: Molecular \& Integrative Physiology} , \emph{168} ,
40-49.

Love, M. I., Huber, W., \& Anders, S. (2014). Moderated estimation of
fold change and dispersion for RNA-seq data with DESeq2. \emph{Genome
Biology} , \emph{15} (12), 550.

Ma, F., Fuqua, B. K., Hasin, Y., Yukhtman, C., Vulpe, C. D., Lusis, A.
J., \& Pellegrini, M. (2019). A comparison between whole transcript and
3'RNA sequencing methods using Kapa and Lexogen library preparation
methods. \emph{BMC Genomics} , \emph{20} (1), 9.

Mandelboum S, Manber Z, Elroy-Stein O, Elkon R (2019) Recurrent
functional misinterpretation of RNA-seq data caused by sample-specific
gene length bias. \emph{PLoS Biology} , \emph{17} (11), e3000481.

Marchuk, L., Sciore, P., Reno, C., Frank, C. B., \& Hart, D. A. (1998).
Postmortem stability of total RNA isolated from rabbit ligament, tendon
and cartilage. \emph{Biochimica et Biophysica Acta (BBA)-General
Subjects} , 1379(2), 171-177.

Meyer, E., Aglyamova, G. V., \& Matz, M. V. (2011). Profiling gene
expression responses of coral larvae (\emph{Acropora millepora} ) to
elevated temperature and settlement inducers using a novel RNA\selectlanguage{english}-Seq
procedure. \emph{Molecular Ecology} , \emph{20} (17), 3599-3616.

Milla, S., Mathieu, C., Wang, N., Lambert, S., Nadzialek, S., Massart,
S., Henrotte, E., Douxfils, J., Melard, C., Mandiki, S.N.M. and
Kestemont, P. (2010). Spleen immune status is affected after acute
handling stress but not regulated by cortisol in Eurasian
perch,\emph{Perca fluviatilis} . \emph{Fish \& Shellfish Immunology}
,\emph{28} (5-6), 931-941.

Milligan, C. L. (1996). Metabolic recovery from exhaustive exercise in
rainbow trout. \emph{Comparative Biochemistry and Physiology Part A:
Physiology} , \emph{113} (1), 51-60.

Morrissy, A. S., Morin, R. D., Delaney, A., Zeng, T., McDonald, H.,
Jones, S., Zhao, Y., Hirst, M. and Marra, M. A. (2009). Next-generation
tag sequencing for cancer gene expression profiling. \emph{Genome
Research} , \emph{19} (10), 1825-1835.

Mutch, D. M., Tordjman, J., Pelloux, V., Hanczar, B., Henegar, C.,
Poitou, C., Veyrie, N., Zucker, J. D. and Clement, K. (2008). Needle and
surgical biopsy techniques differentially affect adipose tissue gene
expression profiles. \emph{The American Journal of Clinical Nutrition}
,\emph{89} (1), 51-57.

Mutter, G. L., Zahrieh, D., Liu, C., Neuberg, D., Finkelstein, D.,
Baker, H. E., \& Warrington, J. A. (2004). Comparison of frozen and
RNALater solid tissue storage methods for use in RNA expression
microarrays. \emph{BMC Genomics} , \emph{5} (1), 88.

Nakatsuji, N., Adachi, K., \& Morioka, K. (2019). Long-term stability of
RNA isolated from muscle of red seabream (\emph{Parus major} ) during
ice storage. \emph{Fish Physiology and Biochemistry} , \emph{45} (2),
819-828.

Narum, S. R., \& Campbell, N. R. (2015) Transcriptomic response to heat
stress among ecologically divergent populations of redband
trout.\emph{BMC Genomics} , 16, 103

Olsvik, P. A., Lie, K. K., \& Hevroy, E. M. (2007). Do anesthetics and
sampling strategies affect transcription analysis of fish
tissues?.\emph{BMC Molecular Biology} , \emph{8} (1), 48.

Oshlack, A., \& Wakefield, M. J. (2009). Transcript length bias in
RNA-seq data confounds systems biology. \emph{Biology Direct} ,\emph{4}
(1), 14.

Passow, C. N., Kono, T. J., Stahl, B. A., Jaggard, J. B., Keene, A. C.,
\& McGaugh, S. E. (2019). Nonrandom RNAseq gene expression associated
with RNAlater and flash freezing storage methods. \emph{Molecular
Ecology Resources} , \emph{19} (2), 456-464.

Pearce, D. A., Arthur, L. M., Turnbull, A. K., Renshaw, L., Sabine, V.
S., Thomas, J. S., Bartlett, J. M., Dixon, J. M. and Sims, A. H. (2016).
Tumour sampling method can significantly influence gene expression
profiles derived from neoadjuvant window studies. \emph{Scientific
Reports} , \emph{6} , 29434.

Perez\selectlanguage{english}-Portela, R., \& Riesgo, A. (2013). Optimizing preservation
protocols to extract high\selectlanguage{english}-quality RNA from different tissues of
echinoderms for next\selectlanguage{english}-generation sequencing. \emph{Molecular Ecology
Resources} , \emph{13} (5), 884-889.

Pertea, M., Kim, D., Pertea, G. M., Leek, J. T., \& Salzberg, S. L.
(2016). Transcript-level expression analysis of RNA-seq experiments with
HISAT, StringTie and Ballgown. \emph{Nature Protocols} , \emph{11} (9),
1650.

Qian, X., Ba, Y., Zhuang, Q., \& Zhong, G. (2014). RNA-Seq technology
and its application in fish transcriptomics. \emph{Omics: a journal of
integrative biology} , \emph{18} (2), 98-110.

R Core Team (2019). R: A language and environment for statistical
computing. R Foundation for Statistical Computing, Vienna, Austria.
URL\emph{https://www.R-project.org/}.

Reese, S. E., Archer, K. J., Therneau, T. M., Atkinson, E. J., Vachon,
C. M., De Andrade, M., Kocher, J. P. A. and Eckel-Passow, J. E. (2013).
A new statistic for identifying batch effects in high-throughput genomic
data that uses guided principal component analysis.\emph{Bioinformatics}
, \emph{29} (22), 2877-2883.

Rey, O., Eizaguirre, C., Angers, B., Baltazar-Soares, M., Sagonas, K.,
Prunier, J. G., Blanchet, S. (2020) Linking epigenetics and biological
conservation: towards a conservation epigenetics
perspective.\emph{Functional Ecology} , 34: 414 -- 427

Roberts, A., Trapnell, C., Donaghey, J., Rinn, J. L., \& Pachter, L.
(2011). Improving RNA-Seq expression estimates by correcting for
fragment bias. \emph{Genome Biology} , \emph{12} (3), R22.

Romero, L. M., \& Reed, J. M. (2005). Collecting baseline corticosterone
samples in the field: is under 3 min good enough?. \emph{Comparative
Biochemistry and Physiology Part A: Molecular \& Integrative Physiology}
, \emph{140} (1), 73-79.

Romero, I. G., Pai, A. A., Tung, J., \& Gilad, Y. (2014). RNA-seq:
impact of RNA degradation on transcript quantification. \emph{BMC
Biology} , \emph{12} (1), 42.

Ross, J. (1995). mRNA stability in mammalian cells. \emph{Microbiol.
Mol. Biol. Rev.} , \emph{59} (3), 423-450.

Roux, J., Rosikiewicz, M., \& Robinson-Rechavi, M. (2015). What to
compare and how: Comparative transcriptomics for Evo-Devo. \emph{J. Exp.
Zool. (Mol. Dev. Evol.)} , 324B, 372--382

RStudio Team (2020). RStudio: Integrated Development for R. RStudio,
PBC, Boston, MA URL http://www.rstudio.com/.

Ruane, N. M., Huisman, E. A., \& Komen, J. (2001). Plasma cortisol and
metabolite level profiles in two isogenic strains of common carp during
confinement. \emph{Journal of fish biology} , \emph{59} (1), 1-12.

Seear, P. J., \& Sweeney, G. E. (2008). Stability of RNA isolated from
post-mortem tissues of Atlantic salmon (\emph{Salmo salar L}
.).\emph{Fish Physiology and Biochemistry} , 34(1), 19-24.

Shepard, B. B., Urie, W., \& May, B. E. (2003). Status of westslope
cutthroat trout (\emph{Oncorhynchus clarkii lewisi} ) in the United
States: \emph{2002} . \emph{USDA Forest Service.}

Smith, S., Bernatchez, L., \& Beheregaray, L. B. (2013). RNA-seq
analysis reveals extensive transcriptional plasticity to temperature
stress in a freshwater fish species. \emph{BMC Genomics} , \emph{14}
(1), 375.

Sopinka, N. M., Donaldson, M. R., O'Connor, C. M., Suski, C. D., \&
Cooke, S. J. (2016). Stress indicators in fish. In \emph{Fish
physiology} (Vol. 35, pp. 405-462). Academic Press.

Staton, J. M., Thomson, A. M., \& Leedman, P. J. (2000). Hormonal
regulation of mRNA stability and RNA-protein interactions in the
pituitary. \emph{Journal of Molecular Endocrinology} , \emph{25} (1),
17-34.

Stolper, J., Ambrosio, E. M., Danciu, D. P., Buono, L., Elliott, D. A.,
Naruse, K., Martinez-Morales, J. R., Marciniak-Czochra, A., \& Centanin,
L. (2019). Stem cell topography splits growth and homeostatic functions
in the fish gill. \emph{eLife, 8} , e43747.

Tandonnet, S., Torres, T. T. (2017) Traditional versus 3\selectlanguage{english}' RNA-seq in a
non-model species. \emph{Genomics Data} , 11, 9-16.

Todd, E. V., Black, M. A., \& Gemmell, N. J. (2016). The power and
promise of RNA\selectlanguage{english}-seq in ecology and evolution. \emph{Molecular Ecology}
,\emph{25} (6), 1224-1241.

Wang, Z., Gerstein, M., \& Snyder, M. (2009). RNA-Seq: a revolutionary
tool for transcriptomics. \emph{Nature Reviews Genetics} , \emph{10}
(1), 57-63.

Wang, Y., Yang, L., Wu, B., Song, Z., \& He, S. (2015). Transcriptome
analysis of the plateau fish (\emph{Triplophysa dalaica} ): implications
for adaptation to hypoxia in fishes. \emph{Gene} , \emph{565} (2),
211-220.

Wang, L., Nie, J., Sicotte, H., Li, Y., Eckel-Passow, J. E., Dasari, S.,
Vedell, P. T., Barman, P., Wang, L., Weinshiboum, R. and Jen, J. (2016).
Measure transcript integrity using RNA-seq data. \emph{BMC
Bioinformatics} , \emph{17} (1), 58.

Wiseman, S., Osachoff, H., Bassett, E., Malhotra, J., Bruno, J.,
VanAggelen, G., Mommsen, T.P. and Vijayan, M. M. (2007). Gene expression
pattern in the liver during recovery from an acute stressor in rainbow
trout. \emph{Comparative Biochemistry and Physiology Part D: Genomics
and Proteomics} , \emph{2} (3), 234-244.

Wong, P. B., Wiley, E. O., Johnson, W. E., Ryder, O. A., O'Brien, S. J.,
Haussler, D., Koepfli, K. P., Houck, M. L., Perelman, P., Mastromonaco,
G. and Bentley, A. C. (2012). Tissue sampling methods and standards for
vertebrate genomics. \emph{GigaScience} , \emph{1} (1), 2047-217X.

Wood, C. M., Turner, J. D., \& Graham, M. S. (1983). Why do fish die
after severe exercise?. \emph{Journal of Fish Biology} , \emph{22} (2),
189-201.

Wu, Z.J., Meyer, C.A., Choudhury, S., Shipitsin, M., Maruyama, R.,
Bessarabova, M., Nikolskaya, T., Sukumar, S., Schwartzman, A., Liu, J.S.
and Polyak, K. (2010). Gene expression profiling of human breast tissue
samples using SAGE-Seq. \emph{Genome Research} , \emph{20} (12),
1730-1739.

Xiong, Y., Soumillon, M., Wu, J., Hansen, J., Hu, B., Van Hasselt, J.G.,
Jayaraman, G., Lim, R., Bouhaddou, M., Ornelas, L., Bochicchio, J.,
Lenaeus, L. \emph{et al.} (2017) A comparison of mRNA sequencing with
random primed and 3-directed libraries. \emph{Scientific Reports} , 7,
14626.

Yu, Y., Fuscoe, J. C., Zhao, C., Guo, C., Jia, M., Qing, T., Bannon,
D.I., Lancashire, L., Bao, W., Du, T. and Luo, H. (2014). A rat RNA-Seq
transcriptomic BodyMap across 11 organs and 4 developmental
stages.\emph{Nature Communications} , \emph{5} (1), 1-11.

\textbf{Data Accessibility Statement}

Raw reads produced for this work are uploaded to NCBI (SRA PRJNA691889)
and will be available after manuscript acceptance. All the supplementary
tables mentioned in this study have been deposited on Dryad and can be
found here:
https://datadryad.org/stash/share/0e5VwzuQzjim9fzQppIl8Q6g7eVqaX2NlaQ22kGetX4

\textbf{Data Citation}

Raw reads and supplementary tables and data. Chiari, Y., Howard, L.,
Moreno, N., Relyea, S., Dunnigan, J., Boyer, M.C., Kardos, M.,
Glaberman, S., Luikart, G. (2022) RNA-Seq reads and analysis outputs.
NCBI (SRA PRJNA691889) and Dryad

\textbf{Legend for the file deposited on Dryad}

\textbf{Table S1: Sample, RNA quality, gene counts, and library
information.} Sheet ``\emph{Samples All} '' lists all samples collected
(sample ID and Admera Health ID for QuantSeq and for NEB) with
information about the treatment group they belong to, tissue type,
sampling method, length and weight of the fish, RIN value, and RNA
concentration. Sample size used for each comparison and divided for
tissue type, treatment group, and library preparation is also indicated.
Sheet ``\emph{QuantSeq} '' lists all samples used for the QuantSeq
library with the following information for each sample: treatment group,
sample ID, Admera Health ID's, tissue type, sampling method, RIN value,
concentration, raw read count, read count after mapping the randomly
selected 11 million reads, and percentage of uniquely mapped genes on
the reference genome. Sheet ``\emph{NEB} '' lists all samples used in
the NEB library detailing for each sample the following: treatment
group, sample ID, Admera Health ID's (and new Admera Health ID if
existing), tissue type, sampling method, RIN value, concentration, raw
read count (PE and single), read count after mapping the randomly
selected 40 million reads, and percentage of uniquely mapped genes on
the reference genome.

\textbf{Table S2: Output results of the Differential Expression
Analysis} . Results of Differential Expression Analysis done with DESeq2
for all comparisons, each of them presented on a separate sheet. Groups
1, 2, and 3 refer to sampling by netting, electrofishing, and
electrofishing with processing 5 min after euthanasia, respectively.

\textbf{Table S3: Summary of gene expression patterns for different
sampling methods and tissue types.} The total numbers of genes with
detectable expression for each sampling/tissue comparison are indicated
along with the number and proportion of genes with significantly higher
gene expression in one of the two tissues being compared for each
sampling method.

\textbf{Table S4: Summary of gene expression patterns for different
sampling methods and tissue types.} The total numbers of genes with
detectable expression for each sampling/tissue comparison are indicated
along with the number and proportion of genes with significantly higher
gene expression in one of the two tissues being compared for each
sampling method.

\textbf{Author contribution}

YC, MB, MK, and GL designed the study. LH, SR, JD and GL performed the
sampling. Analyses were performed by SG, YC, and commissioned to Admera
Health. YC, LH, NM, and SG wrote the manuscript. YC and GL provided
funding for this project. All authors have edited and approved the
manuscript.

\textbf{Table 1: Samples and RNA quality.} List of all the samples on
which transcriptomic data were obtained in this study (sample ID and
Admera Health ID for QuantSeq and NEB) with information about the
treatment group they belong to, tissue type, sampling method, RIN value,
and RNA concentration (ng/ul). Full list of samples and relative
information for which RNA was extracted can be found in Table S1 on
Dryad\selectlanguage{english}
\begin{longtable}[]{@{}llllllll@{}}
\toprule
\textbf{Sample ID} & \textbf{QuantSeq} & \textbf{NEB ID} & \textbf{Group
\#} & \textbf{Tissue} & \textbf{Sampling Method} & \textbf{RIN} &
\textbf{Concentration (ng/ul)}\tabularnewline
\midrule
\endhead
1\_1\_1 & 19101XR-01-41 & & 1 & blood & dip net & 9 &
60.3\tabularnewline
1\_1\_3 & 19101XR-01-42 & & 1 & muscle & dip net & 9.2 &
188.5\tabularnewline
1\_1\_5 & 19101XR-01-44 & & 1 & gill & dip net & 9.3 &
68.9\tabularnewline
1\_10\_1 & 19101XR-01-77 & 19101XR-01-77NEB & 1 & blood & dip net & 9.8
& 166.1\tabularnewline
1\_10\_3 & 19101XR-01-78 & & 1 & muscle & dip net & 8.9 &
586.1\tabularnewline
1\_10\_5 & 19101XR-01-80 & & 1 & gill & dip net & 9.2 &
784\tabularnewline
1\_2\_1 & 19101XR-01-45 & & 1 & blood & dip net & 9.4 &
51.8\tabularnewline
1\_2\_3 & 19101XR-01-46 & & 1 & muscle & dip net & 10 &
150.7\tabularnewline
1\_2\_5 & 19101XR-01-48 & & 1 & gill & dip net & 9.3 &
361.9\tabularnewline
1\_3\_1 & 19101XR-01-49 & & 1 & blood & dip net & 9.2 &
75.3\tabularnewline
1\_3\_3 & 19101XR-01-50 & & 1 & muscle & dip net & 9.6 &
307.3\tabularnewline
1\_3\_5 & 19101XR-01-52 & & 1 & gill & dip net & 9.5 &
69.3\tabularnewline
1\_4\_1 & 19101XR-01-53 & 19101XR-01-53NEB & 1 & blood & dip net & 9.5 &
80.8\tabularnewline
1\_4\_3 & 19101XR-01-54 & & 1 & muscle & dip net & 8.9 &
224.4\tabularnewline
1\_4\_5 & 19101XR-01-56 & & 1 & gill & dip net & 9.2 &
201.8\tabularnewline
1\_5\_1 & 19101XR-01-57 & 19101XR-01-57NEB & 1 & blood & dip net & 9.6 &
195.2\tabularnewline
1\_5\_3 & 19101XR-01-58 & & 1 & muscle & dip net & 9 &
90.7\tabularnewline
1\_5\_5 & 19101XR-01-60 & & 1 & gill & dip net & 8.8 &
44.5\tabularnewline
1\_6\_1 & 19101XR-01-61 & 19101XR-01-61NEB & 1 & blood & dip net & 9.6 &
64.8\tabularnewline
1\_6\_3 & 19101XR-01-62 & & 1 & muscle & dip net & 8.8 &
122.5\tabularnewline
1\_6\_5 & 19101XR-01-64 & & 1 & gill & dip net & 9.2 & 91\tabularnewline
1\_7\_1 & 19101XR-01-65 & & 1 & blood & dip net & 9.4 &
88\tabularnewline
1\_7\_3 & 19101XR-01-66 & & 1 & muscle & dip net & 9.4 &
248.8\tabularnewline
1\_8\_1 & 19101XR-01-69 & 19101XR-01-69NEB & 1 & blood & dip net & 9.4 &
154.3\tabularnewline
1\_8\_3 & 19101XR-01-70 & & 1 & muscle & dip net & 8.8 &
128.5\tabularnewline
1\_8\_5 & 19101XR-01-72 & & 1 & gill & dip net & 9.2 &
575.4\tabularnewline
1\_9\_1 & 19101XR-01-73 & & 1 & blood & dip net & 9.7 &
40.5\tabularnewline
1\_9\_5 & 19101XR-01-76 & & 1 & gill & dip net & 9.3 &
483.3\tabularnewline
3\_1\_1 & 19101XR-01-121 & 19101XR-01-121NEB & 2 & blood & efishing &
9.6 & 105.5\tabularnewline
3\_1\_3 & 19101XR-01-122 & & 2 & muscle & efishing & 8.8 &
384.8\tabularnewline
3\_1\_5 & 19101XR-01-124 & & 2 & gill & efishing & 9.5 &
489.9\tabularnewline
3\_10\_1 & 19101XR-01-157 & 19101XR-01-157NEB & 2 & blood & efishing &
9.8 & 92.7\tabularnewline
3\_10\_3 & 19101XR-01-158 & & 2 & muscle & efishing & 9.4 &
234.3\tabularnewline
3\_10\_5 & 19101XR-01-160 & & 2 & gill & efishing & 9.1 &
526.7\tabularnewline
3\_2\_3 & 19101XR-01-126 & & 2 & muscle & efishing & 9.1 &
388.6\tabularnewline
3\_2\_5 & 19101XR-01-128 & & 2 & gill & efishing & 9.2 &
325\tabularnewline
3\_3\_3 & 19101XR-01-130 & & 2 & muscle & efishing & 9.3 &
241.3\tabularnewline
3\_3\_5 & 19101XR-01-132 & & 2 & gill & efishing & 9.2 &
124.6\tabularnewline
3\_4\_1 & 19101XR-01-133 & 19101XR-01-133NEB & 2 & blood & efishing &
9.6 & 151.3\tabularnewline
3\_4\_3 & 19101XR-01-134 & & 2 & muscle & efishing & 8.8 &
294.2\tabularnewline
3\_4\_5 & 19101XR-01-136 & & 2 & gill & efishing & 9.1 &
158.1\tabularnewline
3\_5\_1 & 19101XR-01-137 & 19101XR-01-137NEB & 2 & blood & efishing &
9.8 & 122\tabularnewline
3\_5\_3 & 19101XR-01-138 & & 2 & muscle & efishing & 9.5 &
225.2\tabularnewline
3\_5\_5 & 19101XR-01-140 & & 2 & gill & efishing & 9.2 &
238.1\tabularnewline
3\_6\_3 & 19101XR-01-142 & & 2 & muscle & efishing & 9.4 &
204.8\tabularnewline
3\_6\_5 & 19101XR-01-144 & & 2 & gill & efishing & 8.8 &
172.2\tabularnewline
3\_7\_1 & 19101XR-01-145 & & 2 & blood & efishing & 9.8 &
195.7\tabularnewline
3\_7\_3 & 19101XR-01-146 & & 2 & muscle & efishing & 9.2 &
329.3\tabularnewline
3\_7\_5 & 19101XR-01-148 & & 2 & gill & efishing & 9.2 &
270.4\tabularnewline
3\_8\_1 & 19101XR-01-149 & & 2 & blood & efishing & 9.8 &
139.5\tabularnewline
3\_8\_3 & 19101XR-01-150 & & 2 & muscle & efishing & 9.4 &
206.8\tabularnewline
3\_8\_5 & 19101XR-01-152 & & 2 & gill & efishing & 9.2 &
299.4\tabularnewline
3\_9\_1 & 19101XR-01-153 & 19101XR-01-153NEB & 2 & blood & efishing &
9.6 & 199.7\tabularnewline
3\_9\_3 & 19101XR-01-154 & & 2 & muscle & efishing & 9.5 &
161.6\tabularnewline
5\_1\_1 & 19101XR-01-201 & & 3 & blood & efishing & 9.1 &
20\tabularnewline
5\_1\_3 & 19101XR-01-202 & & 3 & muscle & efishing & 8.9 &
167.3\tabularnewline
5\_1\_5 & 19101XR-01-204 & & 3 & gill & efishing & 9.4 &
129.7\tabularnewline
5\_10\_1 & 19101XR-01-237 & 19101XR-01-237NEB & 3 & blood & efishing &
9.6 & 133.3\tabularnewline
5\_10\_3 & 19101XR-01-238 & & 3 & muscle & efishing & 9.4 &
642.9\tabularnewline
5\_10\_5 & 19101XR-01-240 & & 3 & gill & efishing & 9.4 &
147.7\tabularnewline
5\_2\_1 & 19101XR-01-205 & & 3 & blood & efishing & 9.7 &
66.6\tabularnewline
5\_2\_3 & 19101XR-01-206 & & 3 & muscle & efishing & 9 &
506.6\tabularnewline
5\_3\_1 & 19101XR-01-209 & 19101XR-01-209NEB & 3 & blood & efishing &
9.6 & 222.1\tabularnewline
5\_3\_3 & 19101XR-01-210 & & 3 & muscle & efishing & 9.4 &
275\tabularnewline
5\_3\_5 & 19101XR-01-212 & & 3 & gill & efishing & 8.9 &
935.2\tabularnewline
5\_4\_1 & 19101XR-01-213 & & 3 & blood & efishing & 9.7 &
71\tabularnewline
5\_4\_3 & 19101XR-01-214 & & 3 & muscle & efishing & 9.1 &
280\tabularnewline
5\_4\_5 & 19101XR-01-216 & & 3 & gill & efishing & 9.2 &
570.5\tabularnewline
5\_5\_5 & 19101XR-01-220 & & 3 & gill & efishing & 9.4 &
370.9\tabularnewline
5\_6\_1 & 19101XR-01-221 & 19101XR-01-221NEB & 3 & blood & efishing &
9.7 & 73.7\tabularnewline
5\_6\_3 & 19101XR-01-222 & & 3 & muscle & efishing & 9.3 &
277.7\tabularnewline
5\_6\_5 & 19101XR-01-224 & & 3 & gill & efishing & 9.2 &
233\tabularnewline
5\_7\_1 & 19101XR-01-225 & & 3 & blood & efishing & 9.4 &
112.2\tabularnewline
5\_7\_3 & 19101XR-01-226 & & 3 & muscle & efishing & 9.6 &
31.7\tabularnewline
5\_7\_5 & 19101XR-01-228 & & 3 & gill & efishing & 9.1 &
96.9\tabularnewline
5\_8\_1 & 19101XR-01-229 & & 3 & blood & efishing & 9.7 &
68.9\tabularnewline
5\_8\_3 & 19101XR-01-230 & & 3 & muscle & efishing & 9.6 &
156.1\tabularnewline
5\_8\_5 & 19101XR-01-232 & & 3 & gill & efishing & 9.5 &
68.7\tabularnewline
5\_9\_1 & 19101XR-01-233 & 19101XR-01-233NEB & 3 & blood & efishing &
9.7 & 212.3\tabularnewline
5\_9\_3 & 19101XR-01-234 & & 3 & muscle & efishing & 9.8 &
60.5\tabularnewline
5\_9\_5 & 19101XR-01-236 & & 3 & gill & efishing & 9 &
115.5\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 2:} \textbf{Sample size used for each comparison} . Sample
size for the samples for which transcriptomic data were obtained divided
for tissue type, treatment group, and library preparation is also
indicated. Group number is as in Table 1: Group 1 -- dip net with tissue
harvested immediately after death; Group 2 -- efishing with tissue
harvested immediately after death; Group 3 -- efishing with tissue
harvested 5 min after death.\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
\textbf{Comparisons} & & &\tabularnewline
\midrule
\endhead
\textbf{Group1 vs Group2 QuantSeq} & \textbf{Blood} & \textbf{Gills} &
\textbf{Muscle}\tabularnewline
Ntot=54 & Ntot=17 & Ntot=18 & Ntot=19\tabularnewline
N1=28 & N1=10 & N1=9 & N1=9\tabularnewline
N2=26 & N2=7 & N2=9 & N2=10\tabularnewline
\textbf{Group2 vs Group3 QuantSeq} & \textbf{Blood} & \textbf{Gills} &
\textbf{Muscle}\tabularnewline
Ntot=53 & Ntot=16 & Ntot=18 & Ntot=19\tabularnewline
N2=26 & N2=7 & N2=9 & N2=10\tabularnewline
N3=27 & N3=9 & N3=9 & N3=9\tabularnewline
\textbf{NEB} & \textbf{Group 1vs Group 2} & \textbf{Group 2vs Group 3}
&\tabularnewline
Ntot = 14 & Ntot=10 & Ntot=9 &\tabularnewline
& N1=5 & N2=5 &\tabularnewline
& N2=5 & N3=4 &\tabularnewline
\textbf{QuantSeq vs NEB} & & &\tabularnewline
Ntot = 28 & & &\tabularnewline
NQuantSeq = 14 & & &\tabularnewline
NNEB=14 & & &\tabularnewline
\bottomrule
\end{longtable}

\textbf{Table 3: Comparisons for which changes in gene expression was
tested} .\selectlanguage{english}
\begin{longtable}[]{@{}llll@{}}
\toprule
& \textbf{Sampling Method} & \textbf{Time of Tissue harvesting} &
\textbf{Tissue types}\tabularnewline
\midrule
\endhead
\textbf{Comparisons} & Dip-netting vs electrofishing & Electrofishing
and processed immediately vs 5 minutes after death & Dip-netting vs
electrofishing\tabularnewline
& & & Immediately vs 5 minutes after death\tabularnewline
\textbf{Library data Used} & QuantSeq & QuantSeq & Only
QuantSeq\tabularnewline
& NEB & NEB &\tabularnewline
& QuantSeq vs NEB & QuantSeq vs NEB &\tabularnewline
\bottomrule
\end{longtable}



\textbf{Figure Legends}

\textbf{Figure 1. Bar plot of transcript length versus number of
non-expressed/non-detected genes for each RNA-Seq library
technique.}Data based on the 14 blood samples processed and sequenced
using both library types. The plot only depicts genes that have been
annotated with a known gene length.

\textbf{Figure 2. Sample-to-sample distance heatmap} . Sample-to-sample
distance heatmaps for the comparison between different sampling
techniques, different tissue harvesting time, and different mRNA-Seq
libraries. The rows and columns are arranged based on hierarchical
clustering, so that samples with similar expression profiles are
positioned near to each other. The color scale represents the distance
between samples. A value of distance 0 indicates that two samples have
identical gene expression. The smaller the distance is, the higher is
the correlation between two samples. Treatment groups (called
``\emph{condition} '') compared are indicated in different colors next
to each heatmap. Condition 1 = fish captured by dip netting, condition 2
= fish captured by electrofishing processed immediately, condition 3 =
fish captured by electrofishing processed 5 minutes after
euthanasia.\textbf{A.} QuantSeq dip netting versus electrofishing for
all tissues combined, \textbf{B.} QuantSeq electrofishing with immediate
sampling versus electrofishing with delayed sampling for all tissues
combined,\textbf{C.} NEB dip netting versus electrofishing only for
blood samples, \textbf{D.} NEB electrofishing with immediate sampling
versus electrofishing with delayed sampling only for blood samples,
\textbf{E.}NEB versus QuantSeq comparisons for dip netting versus
electrofishing only for blood samples, \textbf{F.} NEB versus QuantSeq
comparisons for electrofishing with immediate sampling versus
electrofishing with delayed sampling only for blood samples.

\textbf{Figure 3. PCA plots showing PC1 and PC2 for samples that are
differentially expressed among sampling techniques, tissue harvesting
time, and library preparation methods.} Treatment groups compared are
indicated in different colored symbols next to each PCA plot.\textbf{A.}
QuantSeq dip netting versus electrofishing for all tissues combined,
\textbf{B.} QuantSeq electrofishing with immediate sampling versus
electrofishing with delayed sampling for all tissues
combined,\textbf{C.} NEB dip netting versus electrofishing only for
blood samples, \textbf{D.} NEB electrofishing with immediate sampling
versus electrofishing with delayed sampling only for blood samples.

\textbf{Figure 4. PCA plots showing PC1 and PC2 for samples that are
differentially expressed among sampling techniques, tissue harvesting
time, and library preparation methods.} Treatment groups compared are
indicated in different colored symbols next to each PCA plot.\textbf{A.}
NEB versus QuantSeq comparisons for dip netting versus electrofishing
only for blood samples, \textbf{B.} NEB versus QuantSeq comparisons for
electrofishing with immediate sampling versus electrofishing with
delayed sampling only for blood samples.

\textbf{Figure 5. Violin and box plots comparing gene expression versus
gene length for NEB and QuantSeq library types.} Each individual plot
compares the number of genes with significantly different base mean
expression for NEB versus QuantSeq, calculated as logbasemean NEB -
logbasemean QuantSeq. Genes with equal expression fall on the zero line
of the y-axis; genes with higher expression for the whole mRNA
transcriptome versus QuantSeq have positive numeric values above 0,
while genes with higher expression for QuantSeq vs whole mRNA
transcriptome have negative numeric values below 0. The plot only
depicts genes that have been annotated with a known gene length.

\textbf{Supporting Information Legends:}

\textbf{Figure S1:} \textbf{Sample-to-sample distance heatmap} .
Sample-to-sample distance heatmaps for the comparison between different
sampling techniques and different tissue harvesting time for the
different tissues. The rows and columns are arranged based on
hierarchical clustering, so that samples with similar expression
profiles are positioned near to each other. The color scale represents
the distance between samples. A value of distance 0 indicates that two
samples have identical gene expression. The smaller the distance is, the
higher is the correlation between two samples. Treatment groups compared
(called ``\emph{condition} '') are indicated in different colors next to
each heatmap. Group 1 = fish captured by dip netting, Grou{[} 2 = fish
captured by electrofishing processed immediately, Group 3 = fish
captured by electrofishing processed 5 minutes after
euthanasia.\textbf{A.} QuantSeq dip netting versus electrofishing only
for blood samples, \textbf{B.} QuantSeq dip netting versus
electrofishing only for gill samples, \textbf{C.} QuantSeq dip netting
versus electrofishing only for muscle samples, \textbf{D.} QuantSeq
electrofishing with immediate sampling versus electrofishing with
delayed sampling only for blood samples, \textbf{E.} QuantSeq
electrofishing with immediate sampling versus electrofishing with
delayed sampling only for gill samples, \textbf{F.} QuantSeq
electrofishing with immediate sampling versus electrofishing with
delayed sampling only for muscle samples.

\textbf{Figure S2:} \textbf{PCA plots showing PC1 and PC2 for samples
that are differentially expressed among sampling techniques and tissue
harvesting time for the different tissues. A.} QuantSeq dip netting
versus electrofishing only for blood samples, \textbf{B.} QuantSeq dip
netting versus electrofishing only for gill samples, \textbf{C.}QuantSeq
dip netting versus electrofishing only for muscle samples,\textbf{D.}
QuantSeq electrofishing with immediate sampling versus electrofishing
with delayed sampling only for blood samples, \textbf{E.}QuantSeq
electrofishing with immediate sampling versus electrofishing with
delayed sampling only for gill samples, \textbf{F.} QuantSeq
electrofishing with immediate sampling versus electrofishing with
delayed sampling only for muscle samples.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig1/Fig1}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig2/Fig2}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig3/Fig3}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig4/Fig4}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/Fig/Fig}

\end{center}
\end{figure}

\selectlanguage{english}
\FloatBarrier
\end{document}

