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 (Phillips et al., 2000). 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 (Phillips et al., 2000).

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.

Materials and Methods

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₂, 1.5 mM MgCl₂, 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²) 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 (Wang et al., 2010).

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⁻⁵, 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™ 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^® 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 two-step 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 2013–2015; Zaiontz, 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 Ct-values 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 2013–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.

Results

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.

TABLE 1

Table 1. Read percentages for the 15 most abundant genes in the erythrocyte samples.

For 13 gene IDs GO terms were found by blast search against the NR database. Seven transcripts were associated with the various GO terms related to hemoglobin function (GO:0005833 hemoglobin complex; GO:0005344 oxygen transporter activity; GO:0005506 iron ion binding; GO:0020037 heme binding; GO:0015671 oxygen transport; GO:0019825 oxygen binding). The EBV-induced G-protein coupled receptor 2 (NM_001140225) was assigned to GO terms GO:0016021 (integral component of membrane), GO:0007186 (G-protein coupled receptor signaling pathway), and GO:0002250 (adaptive immune response). The N-terminal EF-hand calcium-binding protein (XM_014146046) was only associated with GO term GO:0005509 (calcium ion binding), while for the transcript of the TBC1 domain family member 15 (XM_014201444) the GO terms GO:0090630 (activation of GTPase activity), GO:0006886 (intracellular protein transport), and GO:0031338 (regulation of vesicle fusion). Stathmin-3 (XM_014168090) was associated with neuron projection development (GO:0031175), cytoplasmic microtubule organization (GO:0031122), regulation of microtubule polymerization or depolymerization (GO:0031110), regulation of GTPase activity (GO:0043087), negative regulation of Rac protein signal transduction (GO:0035021), and microtubule depolymerization (GO:0007019). Serine protease 27 (XM_014212639) is assigned to proteolysis (GO:0006508), and ribosomal RNA small subunit methyltransferase NEP1 (XM_014152620) with methylation (GO:0032259).

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₂-fold changes of transcripts in old erythrocytes did not exceed 2.05 (Table 2). The three genes with the highest log₂-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₂-fold changes of five genes with unknown products were up to −6.5 in old erythrocytes (Table S2).

FIGURE 1

Figure 1. Heat map of age-related changes in gene expression of erythrocytes (RBCs) with hierarchical clustering (complete linkage method with Euclidean distance measure) of all biological replicates. Red indicates highest expression and blue indicates lowest expression.

TABLE 2

Table 2. List of genes with significantly higher transcript levels in old vs. young erythrocyte fractions.

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 ion binding, and in cellular component category, 4 enriched GO terms were found: intracellular (200 IDs), intracellular part (167 IDs), membrane-bounded organelle (105 IDs) and intracellular organelle (191 IDs) were enriched terms in the cellular component.

FIGURE 2

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 (Table 2; pie charts D–F) as input.

The input list of the transcripts with higher levels in old erythrocytes contained 34 IDs (Table 2), for 18 IDs a total of 25 Gene ontology (GO) annotation terms were identified where 7 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₂ 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).

FIGURE 3

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.

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₂-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 half-life 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 non-functional 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 calcium-binding 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 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₃-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.

FIGURE 4

Figure 4. Most frequent GO terms associated with the 742 DEGs with lower transcript levels in old erythrocytes.

Only a small proportion of the DEGs (N = 41; 4.4%) showed higher levels in old erythrocytes. Among them, 34 are protein-coding and only for three the products are currently not known, although their log₂-fold changes are comparably high. Two of the three genes with the highest log₂-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).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

References

Altschul, S. F., Gish, W., Miller, W., Myers, E. W., and Lipman, D. J. (1990). Basic local alignment search tool. J. Mol. Biol. 215, 403–410. doi: 10.1016/S0022-2836(05)80360-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Alvarez, F., Flaño, E., Villena, A. J., Zapata, A., and Razquin, B. E. (1994). Seasonal intrathymic erythropoietic activity in trout. Dev. Comp. Immunol. 18, 409–420.

PubMed Abstract | Google Scholar

Avery, E. H., Lee, B. L., Freedland, R. A., and Cornelius, C. E. (1992). Bile pigments in gallbladder and freshly-secreted hepatic duct bile from fed and fasted rainbow trout, Oncorhynchus mykiss. Comp. Biochem. Physiol. Comp. Physiol. 101, 857–861.

PubMed Abstract | Google Scholar

Bastos, R. N., Volloch, Z., and Aviv, H. (1977). Messenger RNA population analysis during erythroid differentiation: a kinetical approach. J. Mol. Biol. 110, 191–203.

PubMed Abstract | Google Scholar

Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, M. R. (1988). Senescence of red blood cells: progress and problems. Physiol. Rev. 68, 503–554.

PubMed Abstract

Cossins, A. R., and Richardson, P. A. (1985). Adrenalin-induced Na+/H+ exchange in' trout erythrocytes and its effects upon oxygen-carrying capacity. J. Exp. Biol. 118, 229–246.

Götting, M., and Nikinmaa, M. (2017). In vitro study on the regulation of cellular mRNA levels by changes in transcription rate and transcript stability in fish red blood cells. Comp. Biochem. Physiol. Part B Biochem. Mol. Biol. 213, 35–44. doi: 10.1016/j.cbpb.2017.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Götz, S., García-Gómez, J. M., Terol, J., Williams, T. D., Nagaraj, S. H., Nueda, M. J., et al. (2008). High-throughput functional annotation and data mining with the Blast2GO suite. Nucleic Acids Res. 36, 3420–3435. doi: 10.1093/nar/gkn176

PubMed Abstract | CrossRef Full Text | Google Scholar

Grasso, J. A., Chromey, N. C., and Moxey, C. F. (1977). Biochemical characterization of RNA and protein synthesis in erythrocyte development. J. Cell Biol. 73, 206–222.

PubMed Abstract | Google Scholar

Härdig, J., and Hoglund, L. B. (1984). Seasonal variation in blood components of reared Baltic salmon, Salmo salar L. J. Fish Biol. 24, 565–579. doi: 10.1111/j.1095-8649.1984.tb04827.x

CrossRef Full Text | Google Scholar

Hayles, B., Yellaboina, S., and Wang, D. (2010). Comparing transcription rate and mRNA abundance as parameters for biochemical pathway and network analysis. PLoS ONE 5:e9908. doi: 10.1371/journal.pone.0009908

PubMed Abstract | CrossRef Full Text | Google Scholar

Houston, A. H., Dobric, N., and Kahurananga, R. (1996). The nature of hematological response in fish. Fish Physiol. Biochem. 15, 339–347. doi: 10.1007/BF02112361

PubMed Abstract | CrossRef Full Text | Google Scholar

Kim, D., Pertea, G., Trapnell, C., Pimentel, H., Kelley, R., and Salzberg, S. L. (2013). TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions. Genome Biol. 14:R36. doi: 10.1186/gb-2013-14-4-r36

PubMed Abstract | CrossRef Full Text | Google Scholar

Koldkjaer, P., Pottinger, T. G., Perry, S. F., and Cossins, A. R. (2004). Seasonality of the red blood cell stress response in rainbow trout (Oncorhynchus mykiss). J. Exp. Biol. 207, 357–367. doi: 10.1242/jeb.00747

PubMed Abstract | CrossRef Full Text | Google Scholar

Koressaar, T., and Remm, M. (2007). Enhancements and modifications of primer design program Primer3. Bioinformatics 23, 1289–1291. doi: 10.1093/bioinformatics/btm091

PubMed Abstract | CrossRef Full Text | Google Scholar

Lane, H. C. (1979). Some haematological responses of normal and splenectomized rainbow trout Salmo gairdneri to a 12% blood loss. J. Fish Biol. 14, 159–164. doi: 10.1111/j.1095-8649.1979.tb03505.x

CrossRef Full Text | Google Scholar

Lane, H. C., Weaver, J. W., Benson, J. A., and Nichols, H. A. (1982). Some age related changes of adult rainbow trout, Salmo gairdneri Rich., peripheral erythrocytes separated by velocity sedimentation at unit gravity. J. Fish Biol. 21, 1–14. doi: 10.1111/j.1095-8649.1982.tb02818.x

CrossRef Full Text | Google Scholar

Lecklin, T., Tuominen, A., and Nikinmaa, M. (2000). The adrenergic volume changes of immature and mature rainbow trout (Oncorhynchus mykiss) erythrocytes. J. Exp. Biol. 203, 3025–3031.

PubMed Abstract | Google Scholar

Lewis, J. M., Hori, T. S., Rise, M. L., Walsh, P. J., and Currie, S. (2010). Transcriptome responses to heat stress in the nucleated red blood cells of the rainbow trout (Oncorhynchus mykiss). Physiol. Genomics 42, 361–373. doi: 10.1152/physiolgenomics.00067.2010

PubMed Abstract | CrossRef Full Text

Lewis, J. M., Klein, G., Walsh, P. J., and Currie, S. (2012). Rainbow trout (Oncorhynchus mykiss) shift the age composition of circulating red blood cells towards a younger cohort when exposed to thermal stress. J. Comp. Physiol. B. Biochem. Syst. Environ. Physiol. 182, 663–671. doi: 10.1007/s00360-012-0650-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Liao, Y., Smyth, G. K., and Shi, W. (2014). featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics 30, 923–930. doi: 10.1093/bioinformatics/btt656

PubMed Abstract | CrossRef Full Text | Google Scholar

Lien, S., Koop, B. F., Sandve, S. R., Miller, J. R., Kent, M. P., Nome, T., et al. (2016). The Atlantic salmon genome provides insights into rediploidization. Nature 533, 200–205. doi: 10.1038/nature17164

PubMed Abstract | CrossRef Full Text | Google Scholar

Livak, K. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2^−ΔΔCT Method. Methods 25, 402–408. doi: 10.1006/meth.2001.1262

CrossRef Full Text | Google Scholar

Lund, S. G., Phillips, M. C., Moyes, C. D., and Tufts, B. L. (2000). The effects of cell ageing on protein synthesis in rainbow trout (Oncorhynchus mykiss) red blood cells. J. Exp. Biol. 203, 2219–2228.

PubMed Abstract | Google Scholar

Morera, D., Roher, N., Ribas, L., Balasch, J. C., Doñate, C., Callol, A., et al. (2011). RNA-seq reveals an integrated immune response in nucleated erythrocytes. PLoS ONE 6:e26998. doi: 10.1371/journal.pone.0026998

PubMed Abstract | CrossRef Full Text | Google Scholar

Moyes, C. D., Sharma, M. L., Lyons, C., Leary, S. C., Leon, M., Petrie, A., et al. (2002). Origins and consequences of mitochondrial decline in nucleated erythrocytes. Biochim. Biophys. Acta. Mol. Cell Res. 1591, 11–20. doi: 10.1016/S0167-4889(02)00224-0

PubMed Abstract | CrossRef Full Text | Google Scholar

Murphy, J. R. (1973). Influence of temperature and method of centrifugation on the separation of erythrocytes. J. Lab. Clin. Med. 82, 334–41.

PubMed Abstract | Google Scholar

Nikinmaa, M. (1990). Vertebrate Red Blood Cells : Adaptations of Function to Respiratory Requirements. Berlin; New York, NY: Springer-Verlag.

Google Scholar

Nikinmaa, M., and Jensen, F. (1992). Inhibition of adrenergic proton extrusion in rainbow trout red cells by nitrite-induced methaemoglobinaemia. J. Comp. Physiol. B. Biochem. Syst. Environ. Physiol. 162, 424–429. doi: 10.1007/BF00258964

PubMed Abstract | CrossRef Full Text | Google Scholar

Park, K. U., Randazzo, G., Jones, K. L., and Brzezinski, J. A. IV. (2017). Gsg1, Trnp1, and Tmem215 Mark Subpopulations of Bipolar Interneurons in the Mouse Retina. Invest. Ophthalmol. Vis. Sci. 58, 1137–1150. doi: 10.1167/iovs.16-19767

PubMed Abstract | CrossRef Full Text | Google Scholar

Phillips, M. C., Moyes, C. D., and Tufts, B. L. (2000). The effects of cell ageing on metabolism in rainbow trout (Oncorhynchus mykiss) red blood cells. J. Exp. Biol. 203, 1039–1045.

PubMed Abstract | Google Scholar

Speckner, W., Schindler, J. F., and Albers, C. (1989). Age-dependent changes in volume and haemoglobin content of erythrocytes in the carp (Cyprinus carpio L.). J. Exp. Biol. 141, 133–149.

PubMed Abstract | Google Scholar

Tavares-Dias, M. (2006). A morphological and cytochemical study of erythrocytes, thrombocytes and leukocytes in four freshwater teleosts. J. Fish Biol. 68, 1822–1833. doi: 10.1111/j.1095-8649.2006.01089.x

CrossRef Full Text | Google Scholar

Tiano, L., Ballarini, P., Santoni, G., Wozniak, M., and Falcioni, G. (2000). Morphological and functional changes of mitochondria from density separated trout erythrocytes. Biochim. Biophys. Acta 1457, 118–128. doi: 10.1016/S0005-2728(00)00071-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Tonner, P., Srinivasasainagendra, V., Zhang, S., and Zhi, D. (2012). Detecting transcription of ribosomal protein pseudogenes in diverse human tissues from RNA-seq data. BMC Genomics 13:412. doi: 10.1186/1471-2164-13-412

PubMed Abstract | CrossRef Full Text | Google Scholar

Untergasser, A., Cutcutache, I., Koressaar, T., Ye, J., Faircloth, B. C., Remm, M., et al. (2012). Primer3–new capabilities and interfaces. Nucleic Acids Res. 40:e115. doi: 10.1093/nar/gks596

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, L., Feng, Z., Wang, X., Wang, X., and Zhang, X. (2010). DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioinformatics 26, 136–138. doi: 10.1093/bioinformatics/btp612

PubMed Abstract | CrossRef Full Text | Google Scholar

Weber, R. (1990). Functional significance and structural basis of multiple hemoglobins with special reference to ectothermic vertebrates, in Animal Nutrition and Transport Processes Comparative and Environmental Aspects. Comp. Physiol. 6. Transport, Respiration and Excretion, eds. J. Truchot and B. Lahlou (Basel: Karger) 58–75.

Google Scholar

Wiersma, P. A., and Cox, G. S. (1985). Synthesis of messenger-like RNA in avian erythrocyte nuclei. Arch. Biochem. Biophys. 242, 90–103.

PubMed Abstract | Google Scholar

Zaiontz, C. (2015). Real Statistics Using Excel. Available online at: www.realstatistics.com