Transcriptomic Analysis of Young and Old Erythrocytes of Fish

Understanding gene expression changes over the lifespan of cells is of fundamental interest and gives important insights into processes related to maturation and aging. This study was undertaken to understand the global transcriptome changes associated with aging in fish erythrocytes. Fish erythrocytes retain their nuclei throughout their lifetime and they are transcriptionally and translationally active. However, they lose important functions during their lifespan in the circulation. We separated rainbow trout (Oncorhynchus mykiss) erythrocytes into young and old fractions using fixed angle-centrifugation and analyzed transcriptome changes using RNA sequencing (RNA-seq) technology and quantitative real-time PCR. We found 930 differentially expressed between young and old erythrocyte fractions; 889 of these showed higher transcript levels in young, while only 34 protein-coding genes had higher transcript levels in old erythrocytes. In particular genes involved in ion binding, signal transduction, membrane transport, and those encoding various enzyme classes are affected in old erythrocytes. The transcripts with higher levels in old erythrocytes were associated with seven different GO terms within biological processes and nine within molecular functions and cellular components, respectively. Our study furthermore found several highly abundant transcripts as well as a number of differentially expressed genes (DEGs) for which the protein products are currently not known revealing the gaps of knowledge in most non-mammalian vertebrates. Our data provide the first insight into changes involved in aging on the transcriptional level and thus opens new perspectives for the study of maturation processes in fish erythrocytes.


INTRODUCTION
As in all vertebrates apart from mammals, the erythrocytes in fish are nucleated (Nikinmaa, 1990). Throughout their lifespan, which is between 80 and 500 days (Avery et al., 1992), they appear to be capable of aerobic metabolism, although the respiratory rate decreases with age . They are also capable of transcription and translation throughout their lifespan, although, again, the capability decreases with age. The RNA content was significantly lowered by 90% in old erythrocyte fractions while DNA and protein content were unaffected by age (Lund et al., 2000). Young erythrocytes respond much more readily to external stimulation by, e.g., adrenergic drugs than old erythrocytes (Lecklin et al., 2000), and their proportion, which can make up to 50% of the circulating erythrocyte population, increases as a result of temperature stress (Lewis et al., 2012) and changes seasonally (Härdig and Hoglund, 1984). In fish that have recently suffered anemia, the proportion of young RBCs is larger compared to unstressed fish (Lane, 1979).
Associated with the decreased capability to respond to external stimulation, several properties of the erythrocytes are different in young and old erythrocytes. For example, the shape changes from circular to elliptical during maturation (Tavares-Dias, 2006), the membrane of young erythrocytes is more fluid (Lecklin et al., 2000), and old erythrocytes have fewer organelles such as mitochondria (Moyes et al., 2002) and free ribosomes .
However, despite these differences in the function and properties of young and old erythrocytes, and although there is a clear difference between the transcriptome of control (11 • C) and heat-shocked (1 h at 25 • C) erythrocytes (Lewis et al., 2010), quantitative changes in transcription between young and old fish erythrocytes have not been studied earlier. Here we have collected blood from rainbow trout in spring (May) when the proportion of young erythrocytes in the circulation increases with rising ambient water temperature (Härdig and Hoglund, 1984;Alvarez et al., 1994;Houston et al., 1996;Lecklin et al., 2000). We have divided the erythrocytes in young and old cohorts on the basis of the old erythrocytes having higher density than young ones (Lane et al., 1982;Tiano et al., 2000). Since old erythrocytes have higher mean cellular hemoglobin concentration (MCHC) than young ones (Lane et al., 1982), the success of age separation was checked by evaluating MCHC. Thereafter RNA sequencing was carried out and differences in transcriptomes in the young and old cohort evaluated.

Animals and Blood Sampling Procedure
Rainbow trout (Oncorhynchus mykiss, N = 5, weight 548.0 ± 70.1 g) were obtained from a commercial hatchery (Finnish Institute for Fisheries and Environment, Parainen, Finland) in May 2016. All procedures were approved by the Finnish Animal Experiment Board (ESAVI/3705/04.10.07/2015). Since the water in the fish tanks of the hatchery is pumped from a nearby bay of the Finnish Archipelago, the water temperature in the tanks follows natural rhythms and is on the rise at this time of the year (K. Malmberg, personal communication). At the date of sampling the water temperature was 14 • C.
Fish were netted from the tanks, quickly killed by a blow on the head and blood was sampled by caudal puncture into heparinized syringes, transferred into sterile falcon tubes, and stored on ice. Blood samples were washed three times by repeated re-suspension in rainbow trout saline (128 mM NaCl, 3 mM KCl, 1.5 mM CaCl 2 , 1.5 mM MgCl 2 , 20 mM Tris-HCl, pH 7.6 (Nikinmaa and Jensen, 1992) and centrifugation at 800 × g and 10 • C for 3 min to remove buffy coat and white blood cells and stored on ice. The erythrocytes were re-suspended in fresh saline at a hematocrit (Hct) of 18-20% and then stored well-aerated over night at 14 • C in cell culture flasks (75 cm 2 ) with open caps to allow equilibration. Fish were weighed and their length was measured.

Erythrocyte Age Class Separation by Density Centrifugation
Density centrifugation of erythrocyte samples followed essentially the procedure described previously (Murphy, 1973;Speckner et al., 1989;Koldkjaer et al., 2004) with minor modifications. A subsample (referred to as "original blood sample") was taken and immediately frozen at −80 • C. The rest of the erythrocyte sample was once more washed in saline and then adjusted to a hematocrit of ∼80% and transferred into polypropylene tubes (length 47 mm, diameter 4 mm, volume 0.5 ml). Tubes were centrifuged in a fixed-angle rotor (45 • ) at 16,000 × g and 10 • C for 30 min. The tubes were cut into three equally sized parts: the top part, containing the youngest and least dense erythrocyte class fraction, the middle part which was discarded, and the lower part, containing the oldest and most dense erythrocytes. Erythrocyte fractions (50 µl) were transferred into new tubes, diluted with 100 µl rainbow trout saline and frozen at −80 • C.
The success of age separation was evaluated by determining the mean cellular hemoglobin concentration (MCHC) of both fractions using conventional methods (centrifugation for haematocrit and cyanmethaemoglobin method for hemoglobin concentration; Speckner et al., 1989;Lund et al., 2000). The MCHC of the putatively old erythrocytes was significantly higher than that of the putatively young erythrocytes (p = 0.015) indicating that the age separation was successful.

RNA-Seq of Young and Old Erythrocyte Fractions
Erythrocyte samples of different age were then homogenized in TissueLyser (Qiagen, Austin, USA) with two stainless steel beads for 2 × 60 s at 30 Hz and RNA was isolated using NucleoZol (Macherey-Nagel) according to the manufacturer's instructions. RNA was quantified using a NanoDrop 2000 (Thermo Scientific, Bonn, Germany) and quality was checked with a Fragment Analyzer TM . Only samples with OD 260/280 and OD 260/230 > 1.8 and RIN values higher than 8.7 (range 8.7-9.8) were used in the analyses.
Sequencing libraries were prepared using the Illumina TruSeq Stranded mRNA Sample Preparation Kit and sequenced in two lanes on a HiSeq 3000 instrument at the Finnish Functional Genomics Centre (Turku, Finland) and single-end sequencing chemistry with 50 bp read length.

Bioinformatic Analysis
Base calling on the reads was done using the Bcl2fastq2 software (version 1.8.4). The quality control of raw sequencing reads was performed with FastQC (www.bioinformatics.babraham.ac. uk/projects/fastqc/), and adapters and low quality bases were trimmed by Trimmomatic (Bolger et al., 2014). The read alignment was performed against the reference genome using TopHat v2.1.0 (Kim et al., 2013). We used the Atlantic Salmon genome [ICSASG v2 genome; provided by the International Cooperation to Sequence the Atlantic SalmonGenome (ICSASG; Lien et al., 2016)] as a reference genome because of the status of the higher coverage of the assembly compared to that of the rainbow trout. The number of uniquely aligned reads was between 2.6 and 6.8 M per sample.
The sample correlation values (Spearman's metrics) were between 0.895 and 0.932 (mean 0.910 for young age class, 0.927 for old age class) for all samples.
The gene-wise read counts were normalized using the TMM normalization method of the edgeR R/Bioconductor package before further statistical testing using the Limma R/Bioconductor package. Differentially expressed genes (DEGs) were selected based on a fold change >2 and an FDR (false-discovery-rate) <0.01. Sample pairing (young and old erythrocyte age classes of one individual) was taken into account when building the linear model for statistical testing.
The association between reads and known genes and the number of reads associated with each gene was assessed using the subreads package v1.5.1 (Liao et al., 2014). Differential expression between age classes was analyzed using MA-plots as implemented in DEGseq R package .
The BLAST2GO software (Götz et al., 2008) was used to predict Gene Ontology (GO) terms for all statistically significant DE genes (FDR < 0.01) and the most abundant genes. Gos are organized hierarchically in terms of biological processes, cellular components and molecular functions. All pseudogenes (labeled with NA) and all tRNAs were removed from the list of DEGs, and the remaining list was subsequently searched against the non-redundant NCBI protein database (NR database) using the BLASTx algorithm (Altschul et al., 1990) with an E-value threshold of 10 −5 , and a maximum of 20 "hit" sequences per query was retained. GO annotations were simplified to a smaller set of high-level GO terms using GO slims as implemented in BLAST2GO.

Differential Expression of Selected Genes
Using Quantitative Real-Time PCR (qPCR) qPCR was used to analyze the expression of selected genes on preparations of young and old erythrocytes relative to the original erythrocyte sample (taken before the fractionation procedure) of five individuals. We selected eight genes which on the basis of earlier research have important functions in fish erythrocytes and designed specific primers (Table S1) using the Primer3 software (Koressaar and Remm, 2007;Untergasser et al., 2012) and checked for secondary structures using Beacon Designer Software TM and UNAFold tool (https://eu.idtdna.com/ UNAFold). qPCR primers were blasted against the rainbow trout database (NCBI) to ensure specificity.
RNA for qPCR was digested using DNase I (Promega) and 500 ng were reverse transcribed using the RevertAid First Strand Synthesis Kit (Thermo Scientific, Bonn, Germany) using random hexamer primers according to the manufacturer's instructions. cDNA products were amplified in triplicates using the KAPA SYBR R Fast qPCR kit (KapaBiosystems) on a QuantStudio 12K Flex Real-Time PCR System (Thermo Fisher Scientific). Each 10 µL reaction mixture contained 2 µl cDNA template (1:10 or 1:20 dilution) and 0.5 µM of each primer. A twostep cycling protocol was applied: 10 min at 95 • C followed by 40 cycles of 15 s at 95 • C and 30 s at 60 • C and 30 s at 72 • C. In a final step, specificity of primers and amplification was evaluated using dissociation curves with a temperature range from 60 to 95 • C. All primer pairs generated a single peak in the dissociation curve and PCR efficiency estimated for each primer pair was within the range of 92-105%. Each qPCR plate contained non-template controls to detect potential contamination in reaction mixes. Data were analyzed with the QuantStudio software. Reaction efficiency for each gene was calculated using a standard curve generated from a 1:2 serial dilution of the pooled samples. Standard curve reactions were performed in duplicate. Calculations of relative expression levels with the 2 − CT method (Livak and Schmittgen, 2001) were done in Microsoft Office Excel and the expression of target genes was normalized against β-actin (act) expression levels. Stability of β-actin expression between young, old and original blood sample was assessed using a one-way repeated measures ANOVA (p = 0.106) in Real Statistics Resource Pack software (Release 4.3;www.realstatistics.com;Copyright 2013Copyright -2015Zaiontz, 2015). Data are shown relative to the respective mixed blood samples and log2 transformed. All data are expressed as means ± SD. Significant differences were assessed first between raw Ctvalues of the young, old and originalsamples, and second between the raw Ct-values of the young and old erythrocyte samples using a one-way repeated measures ANOVA after equal variances (Brown-Forsythe) and normality (Shapiro-Wilk test) of data were tested using the Real Statistics Resource Pack software (Release 4.3;www.realstatistics.com;Copyright 2013Copyright -2015 in Excel 2010 (Zaiontz, 2015). In cases where tests for equal variances and normality failed we used a non-parametric test (Friedman's test). Post-hoc analyses were done with the Holm-Sidak test. Significance was accepted at the P < 0.05 level.

Sequencing Statistics
In this study, an average number of 12,197,337 reads (mean young erythrocytes: 11,692,319; mean old erythrocytes 12,702,356) were generated. From the 8,273,888 mapped reads (mean young erythrocytes: 8,192,383; mean old erythrocytes 8,355,393), 54% (young erythrocytes: 49%; mean old erythrocytes 59%) could be uniquely mapped to the reference genome. Only uniquely mapped reads (average 4,469,945; mean young erythrocytes 3,956,287; mean old erythrocytes 4,983,602) were used for the downstream analysis, resulting in 51,226 unigenes.

Most Abundant Genes in Erythrocytes
The 15 most abundant genes (those among the top 10 abundant genes in any sample) account for 9.81-29.35% of all reads (mean young erythrocytes = 23.01 ± 5.76%; mean old erythrocytes = 12.39 ± 1.58%; Table 1). Naturally, various hemoglobin subunit transcripts are among them; hbα transcripts accounted for 3.46% and hbß for 4.56% of all reads. The products of two abundant transcripts in the list are currently not known.
For 13 gene IDs GO terms were found by blast search against the NR database. Seven transcripts were associated

Aging Results in Lower Transcript Levels
The DEGs between young and old RBCs were identified by selecting a threshold of > two-fold change in expression level (FDR < 0.01). There was a large shift in transcript abundance between young and old erythrocyte samples (Figure 1). The total number of DEGs was 930; 889 (95.6%) of which show higher transcript levels in young erythrocytes (Figure 1; Table S2), while only 41 (4.4%) have higher transcript levels in old erythrocytes (Figure 1; Table 2). Positive log 2 -fold changes of transcripts in old erythrocytes did not exceed 2.05 ( Table 2). The three genes with the highest log 2 -fold changes were the transmembrane protein 215 (XM_014194795), the cancer susceptibility candidate 1, transcript variant X2 (XM_014153053), and the corin, serine peptidase (XM_014206310). Negative log 2 -fold changes of five genes with unknown products were up to −6.5 in old erythrocytes (Table S2).

Functional Classification of DEGs
The gene annotation tool BLAST2GO (Götz et al., 2008), was used to annotate identified DEGs. The input list for the transcripts with higher levels in young erythrocytes contained 742 IDs (Table S2) and for 478 IDs GO annotation terms were found (Figures 2A-C). Eleven terms belonged to the biological process group (Figure 2A), 6 GO terms to the cellular component group ( Figure 2C) and 5 GO terms to the molecular function group (Figure 2B). The largest three subcategories in the biological group were cellular processes (209 IDs), metabolic processes (151 IDs) and single-organism processes (138 IDs). In the molecular function category 140 IDs were associated with FIGURE 2 | Gene ontology (GO) terms associated with young (A-C) and old (D-F) erythrocyte age classes. Pie charts were generated at level 2 for biological process (A,D), and level 3 for molecular function (B,E) and cellular component (C,F) by Blast2GO using 742 transcripts with higher levels in young erythrocytes (Table S2; pie charts A-C) and 34 transcripts with higher levels in old erythrocytes ( were within biological process (Figure 2D), and 9 each in cellular component ( Figure 2F) and molecular function ( Figure 2E). For 16 IDs no GO terms could be found. GO terms were enriched for cellular process (8 IDs), single-organism process (5 IDs) and metabolic process (4 IDs) in the biological process category. In the metabolic function category, no enriched GO terms were found and with each of the GO terms 1-3 IDs were associated. Intracellular components (7 IDs), membrane-bounded organelle (4 IDs) and intracellular organelle (6 IDs) were enriched terms in the cellular component.

Expression Differences between Young and Old Erythrocyte Fractions in Relation to the Original Erythrocyte Sample
We selected eight genes whose products have important functions in fish erythrocytes for a comparison of expression levels between young and old erythrocyte fractions with the original erythrocyte sample with qPCR. Raw data were first normalized to the expression levels of β-actin, and then related to the values of the original erythrocyte sample and finally log 2 transformed. The levels of expression of the genes are shown in Figure 3. The statistical analyses between young, old and original erythrocyte samples revealed four genes with significant differences: gbx (p = 0.0408), β 3b -ar (p = 0.0111), gapdh (p = 0.0150), and hsp70 (p = 0.007). When only young and old erythrocyte fractions were compared, three of the genes remained significantly different: β 3b -ar (p = 0.0159), gapdh (p = 0.0253), and hsp70 (p = 0.0351; Figure 3).

DISCUSSION
In fish, circulating erythrocytes are a population of different age classes, and their proportion changes with season and as a result of environmental stress (Härdig and Hoglund, 1984;Lewis et al., 2012). Maturation of fish erythrocytes in the circulation is associated with marked changes in many aspects of their physiology. We analyzed these changes on the transcriptomic level in young and old erythrocytes in order to get insights into the yet unknown pathways involved in aging. Changes in cellular transcript levels are the result of a number of tightly controlled processes, such as alteration of transcription rate, transcript stability, or mRNA degradation, and the importance of each of these regulatory steps for the observed changes in transcript levels is not yet known (Götting and Nikinmaa, 2017). Measuring cellular steady-state levels as done by RNA-seq does not take into account alterations in either of these processes (Hayles et al., 2010). Thus, we avoid the terms "up-and downregulated" for changes in transcript levels throughout the manuscript, as it is not known if any of the steps is regulated or just changes.
Our results revealed that aging in fish erythrocytes is accompanied by a decrease in transcript levels of a number of genes. It was shown earlier that the rate of transcription and the cellular RNA content decreases during maturation of nucleated non-mammalian erythrocytes (Grasso et al., 1977;Wiersma and Cox, 1985;Lund et al., 2000). High transcriptional and translational activity is required in young, immature erythrocyte because they gradually undergo changes in cell shape, membrane rigidity and other properties until maturation (e.g., Lecklin et al., 2000;Tavares-Dias, 2006). In contrast, in old, mature erythrocytes the need for gene transcription and protein synthesis is much lower and the necessary machinery, such as ribosomes is very low (Lane et al., 1982;Phillips et al., 2000). A total of 51,226 unigenes could be identified in the present study but we found only 930 (1.8%) DEGs.
Most of the DEGs (95.6%) had significantly decreased levels in old erythrocytes and the overall gene expression is much more variable between young erythrocyte samples than between old erythrocyte samples, which appear more homogeneous (Figure 1). The higher transcriptional variability in young erythrocytes is associated with their more active roles in various physiological processes such as in the perception of external stimuli, adaptation and the immune response compared to the terminally differentiated, mature erythrocyte (Lecklin et al., 2000;Koldkjaer et al., 2004;Morera et al., 2011).
The most abundant transcripts in terms of number and expression level belonged to hemoglobin (Hb). Fish exhibit a remarkable multiplicity of Hbs with different functional properties, such as differences in O 2 -affinity and sensitivity to allosteric regulators, which provide an important molecular strategy for adapting to a changing environment (Weber, 1990). During aging some hemoglobin subunit transcripts decrease significantly, while others remain at the same level. It has been reported earlier that Hb accumulates in erythrocytes during maturation and older erythrocytes are tightly packed with Hb (Bastos et al., 1977;Clark, 1988). Older erythrocytes thus serve merely as oxygen transporters. Even if the overall hb transcript levels are decreased in older erythrocytes, a long transcript halflife and a high stability of Hb protein ensures adequate synthesis for required function throughout the lifespan in the circulation.
Surprisingly, our study revealed that for two of the 15 most abundant transcripts in the samples the products and thus the function is not known. A blast search revealed a total of 7 transcript variants for the locus LOC106601081, and a blastx search using the unknown protein product of transcript variant 1 (XM_014192985) as query, found no hits. Further studies are required to identify and characterize the function of this transcript.
The second unknown abundant transcript is a pseudogene. A number of pseudogenes (N = 1993) were found to be present in the erythrocyte samples but more than half of them were only little expressed. Pseudogene transcripts are very often nonfunctional copies of real genes which have lost some or all of their functionality by accumulating mutations that disable, e.g., translation. Although several studies discussed the existence of pseuogenes in various organisms, another source of pseudogenes is the read alignment step during RNA-seq data analysis (Tonner et al., 2012).
Other abundant transcripts have functions in the intracellular protein transport and GTPase activation (TBC1 domain family member 15), calcium ion binding (N-terminal EF-hand calciumbinding protein) or the cytoplasmic microtubule organization (Stathmin-3). Furthermore, highly abundant transcripts of the serine protease 27, responsible for proteolysis, and of the ribosomal RNA small subunit methyltransferase NEP1, involved in methylation, were found. Among the genes with abundant levels only the transcripts of genes encoding EBV-induced G-protein coupled receptor 2 are significantly reduced by 88% in old erythrocytes, suggesting an impairment of the G-protein coupled receptor signaling pathway and a massively reduced adaptive immune response in old erythrocytes. An earlier study revealed already that nucleated erythrocytes play an active role in FIGURE 3 | Comparative gene expression levels between young and old erythrocytes (RBCs). Expression of eight target genes was normalized to β-actin as a reference gene and then related to the expression of the respective gene in the original erythrocyte sample taken before the fractionation procedure. Gene abbreviations are hbα-1, hemoglobin α-1; hbα-4, hemoglobin α-4; gbx, globin x1; β 3b -ar, β 3b -adrenergic receptor; β-nhe, β-Na + /H + -exchanger; hif1a, hypoxia inducible factor 1a; gapdh, glyceraldehyde-3-phosphate dehydrogenase; hsp70, inducible heat shock protein 70 kDa. Statistically significant differences between means of young and old erythrocyte fractions is shown by *p < 0.05. the immune response (Morera et al., 2011), but according to our results this is probably mostly restricted to young erythrocytes.
Age-dependent decreases in transcript levels of proteins involved in many molecular pathways were found in the present study (Figure 4). In particular genes whose products have functions in ion binding and DNA-binding or in signaling and membrane transport have significantly reduced transcript levels in old erythrocytes. Several enzyme classes (e.g., kinases, oxidoreductases, peptidases) seem to be affected as well. Earlier studies found that mRNA levels of carbonic anhydrase and Cl − /HCO 3 -exchanger (AE1) decrease with age (Lund et al., 2000), and also enzyme activity levels of carbonic anhydrase, citrate synthase, cytochrome oxidase and lactate dehydrogenase in rainbow trout RBCs (Lund et al., 2000;Phillips et al., 2000) are lower.
Only a small proportion of the DEGs (N = 41; 4.4%) showed higher levels in old erythrocytes. Among them, 34 are proteincoding and only for three the products are currently not known, although their log 2 -fold changes are comparably high. Two of the three genes with the highest log 2 -fold changes are integral components of the membrane. Currently little is known about the transmembrane protein 215 (XM_014194795) and its function. It is conserved among vertebrate species and is expressed in mouse bipolar cells (Park et al., 2017). An attractive speculation is, because of the increase of its transcription in old cells, that it is involved in the age-dependent removal of cells from circulation. For the corin, serine peptidase (XM_014206310 ;  Table 1) the blastx search against the NR database revealed a hit of the atrial natriuretic peptide-converting enzyme of Oreochromis niloticus, which is involved in proteolysis and serine-type endopeptidase activity. For the cancer susceptibility candidate 1, transcript variant X2 (XM_014153053) no blastx hits and thus no annotations were found. For almost half of the DEGs with higher levels in old erythrocytes (N = 16; 47%) no GO terms could be retrieved from the databases revealing the gaps in knowledge for non-mammalian non-model organisms.
Using quantitative real-time PCR we analyzed the changes of eight genes, which have important functions in fish erythrocytes, in young and old fractions and compared them to the original erythrocyte sample. Four transcripts (β 3b -ar, gapdh, gbx, and hsp70) had significantly different levels between the young, old and the original sample. When only young and old fractions were compared, three transcripts were found to be significantly different (β 3b -ar, gapdh, and hsp70; Figure 3), which are not among the DEGs in the RNA-seq results (Table S2). This is due to the different statistical cutoffs applied for the RNA-seq data (FDR adjusted p-value, see Materials and Methods section) vs. a significance level of p < 0.05 used in the qPCR. When using the not-adjusted T-test p-value for the RNA-seq data, β 3b -ar and gapdh also have significant different levels in young and old fractions (p = 0.02 and p = 0.03, respectively), while hsp70 remained non-significant (p = 0.05). The other transcripts (hbα-1, hbα-4, hif1a, and β-nhe) studied showed the same results as the RNA-seq data and were not different between age class fractions.
The higher β 3b -ar transcript levels in young fractions and the decrease with age is in accordance with earlier reported results and explains why the response of older, mature erythrocytes to external stimuli such as catecholamines is impaired (Cossins and Richardson, 1985;Lecklin et al., 2000). β 3b -ar transcript levels in rainbow trout vary with season, showing high levels in summer (May to October) and low levels in winter (December to May) (Koldkjaer et al., 2004). Although β 3b -ar mRNA was not compared between young and old erythrocyte fractions by Koldkjaer et al. (2004), the observed increase in mRNA can be mainly attributed to the increased proportion of young erythrocytes in the in spring.
Older erythrocytes have a lower ability to produce stress proteins during the recovery period after heat shock (Lund et al., 2000). We found significantly higher transcript levels for the heat shock protein 70 kDa in young erythrocytes in both qPCR and RNA-seq data, but transcript levels of HSF-1 and Hsc71 were unaffected by cell age. This is in contrast to earlier study reporting lower transcript levels for the heat shock cognate 71 kDa (Hsc71) and the heat shock factor (HSF) in old erythrocytes compared to younger ones (Lund et al., 2000).

CONCLUSION
The mechanisms and characteristics which finally promote the removal of erythrocytes from the circulation and their destruction are still not known in fish. We found only a small proportion of DEGs between young and old fractions suggesting that even rather subtle changes in gene expression can have significant impact on a wide variety of molecular functions and the cell fate. The diversity of transcripts found in both erythrocyte age classes implies that the physiological roles of nucleated erythrocytes go far beyond the already known functions in oxygen transport and immune response. Hence, our data provide a rich resource for future investigations on the regulation of aging in fish and in particular fish erythrocytes and help to uncover additional physiological functions.

DATA DEPOSITION
The RNA-seq data have been deposited at the NCBI Gene Expression Omnibus (GEO) with the accession number GSE106570.

AUTHOR CONTRIBUTIONS
MG and MN: conceived the study and wrote the paper; MG: performed the experiments and analyzed the data.

FUNDING
Financial support for this project was provided by a Marie Curie Fellowship (FP7; Agreement No 623338) and The Ella and Georg Ehrnrooth Foundation to MG and the Academy of Finland to MN (grant no. 258078).

ACKNOWLEDGMENTS
We thank Kurt Malmberg (Ammattiopisto Livia/Yrkesinstitut Livia) for providing the rainbow trout. Additionally, we thank V. Vainio and J. Ruohomäki for help with lab work. We acknowledge the Finnish Functional Genomics Centre and the Bioinformatics Unit at the Turku Centre for Biotechnology and Biocenter Finland for their help in plate run services, RNA sequencing and in data analysis. The Units are supported by University of Turku, Åbo Akademi University, and Biocenter Finland. We would like to thank the two reviewers for their valuable comments on an earlier version of the article.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys. 2017.01046/full#supplementary-material Table S1 | Primers used in qPCR. Table S2 | List of genes with significantly higher transcript levels in young vs. old erythrocyte fractions.