In life, genetic and epigenetic networks precisely coordinate the expression of genes—but in death, it is not known if gene expression diminishes gradually or abruptly stops or if specific genes and pathways are involved. We studied this by identifying mRNA transcripts that apparently increase in relative abundance after death, assessing their functions, and comparing their abundance profiles through postmortem time in two species, mouse and zebrafish. We found mRNA transcript profiles of 1063 genes became significantly more abundant after death of healthy adult animals in a time series spanning up to 96 h postmortem. Ordination plots revealed non-random patterns in the profiles by time. While most of these transcript levels increased within 0.5 h postmortem, some increased only at 24 and 48 h postmortem. Functional characterization of the most abundant transcripts revealed the following categories: stress, immunity, inflammation, apoptosis, transport, development, epigenetic regulation and cancer. The data suggest a step-wise shutdown occurs in organismal death that is manifested by the apparent increase of certain transcripts with various abundance maxima and durations.
A healthy adult vertebrate is a complex biological system capable of highly elaborate functions such as the ability to move, communicate and sense the environment—all at the same time. These functions are tightly regulated by genetic and epigenetic networks through multiple feedback loops that precisely coordinate the expression of thousands of genes at the right time, in the right place and in the right level . Together, these networks maintain homeostasis and thus sustain ‘life’ of a biological system.
While much is known about gene expression circuits in life, there is a paucity of information about what happens to these circuits after organismal death. For example, it is not well known whether gene expression diminishes gradually or abruptly stops in death—nor whether specific gene transcripts increase in abundance in death. In organismal ‘death’, defined here as the cessation of the highly elaborate system functions in vertebrates, we conjecture that there is a gradual disengagement and loss of global regulatory networks as well as the activation of regulatory genes involved in survival and stress compensation. To test this, we examined the global postmortem abundances of mRNAs in two model organisms: the zebrafish, Danio rerio, and the house mouse, Mus musculus. The purpose of the research was to investigate the ‘unwinding of the clock’ by identifying mRNA transcripts that increase in abundance with postmortem time and assessing their functions based on the primary literature. The biological systems investigated in this study are different from those examined in other studies, such as individual dead and/or injured cells in live organisms, i.e. apoptosis and necrosis (reviewed in [2–5]). In contrast to previous studies, the abundances of mRNA transcripts from the entire D. rerio body, and the brains and livers of M. musculus were assessed through postmortem time. The mRNA transcripts were measured using the ‘Gene Meter’ approach that precisely reports transcript abundances based on a calibration curve for each microarray probe [6–9].
2. Material and methods
2.1. Induced death and postmortem incubation
Forty-four female Danio rerio were transferred from several flow-through aquaria kept at 28°C to a glass beaker containing 1 l of aquarium water. Four individuals were immediately taken out, snap frozen in liquid nitrogen and stored in Falcon tubes at −80°C (two zebrafish per tube). These samples were designated as the first set of live controls. A second set of live controls was immersed in an open cylinder (described below). Two sets of live controls were used to determine whether putting the zebrafish back into their native environment had any effects on gene expression (we later discovered no significant effects).
The rest of the zebrafish were subjected to sudden death by immersion in a ‘kill’ chamber. The chamber consisted of an 8 l styrofoam container filled with chilled ice water. To synchronize the death of the rest of the zebrafish, they were transferred to an open cylinder with a mesh-covered bottom and the cylinder was immersed into the kill chamber. After 20–30 s of immersion, four zebrafish were retrieved from the chamber, snap frozen in liquid nitrogen and stored at −80°C (two zebrafish per Falcon tube). These samples were designated as the second set of live controls. The remaining zebrafish were kept in the kill chamber for 5 min and then the cylinder was transferred to a flow-through aquarium kept at 28°C so that they were returned to their native environment.
Postmortem sampling of the zebrafish occurred at: time 0, 15 min, 30 min, 1 h, 4 h, 8 h, 12 h, 24 h, 48 h and 96 h. For each sampling time, four expired zebrafish were retrieved from the cylinder, snap frozen in liquid nitrogen and stored at −80°C in Falcon tubes (two zebrafish to a tube). One zebrafish sample was lost, but extraction volumes were adjusted to one individual.
The mouse strain C57BL/6JRj (Janvier SAS, France) was used for our experiments. The mice were 20-week old males of approximately the same weight. The mice were highly inbred and were expected to have a homogeneous genetic background. Prior to euthanasia, the mice were kept at room temperature and were given ad libitum access to food and water. Each mouse was euthanized by cervical dislocation and placed in an individual plastic bag with holes to allow air/gas exchange. The bagged carcasses were kept at room temperature in a large, open polystyrene container. Sampling of the deceased mice began at 0 h (postmortem time zero) and continued at 30 min, 1 h, 6 h, 12 h, 24 h and 48 h postmortem. At each sample time, three mice were sampled (except for 48 h when two mice were sampled) and the entire brain (plus stem) and two portions of the liver were extracted from each mouse. For liver samples, clippings were taken from the foremost and rightmost lobes of the liver. The brain and liver samples were snap frozen in liquid nitrogen and stored individually in Falcon tubes at −80°C.
2.2. RNA extraction, labelling, hybridization and DNA microarrays
The number of individuals was 43 for zebrafish and 20 for mice. Samples from two fish were pooled for analysis, resulting in two replicate measurements at each time point. The number of replicated measurements for mice was three at each of the first six time points and two at 48 h. Thus, the total number of samples analysed was 22 for zebrafish and 20 for mice. For the zebrafish, samples were mixed with 20 ml of Trizol and homogenized using a TissueLyzer (Qiagen). For the mice, 100 mg of brain or liver samples were mixed with 1 ml of Trizol and homogenized. One millilitre of the emulsion from each sample was put into a fresh 1.5 ml centrifuge tube for RNA extraction and the rest was frozen at −80°C.
RNA was extracted by adding 200 µl of chloroform, vortexing the sample and incubating it at 25°C for 3 min. After centrifugation (15 min at 12 000g at 4°C), the supernatant (approx. 350 µl) was transferred to a fresh 1.5 ml tube containing an equal volume of 70% ethanol. The tube was vortexed, centrifuged and purified following the procedures outlined in the PureLink RNA Mini Kit (Life Technologies, USA).
The isolated RNA, 400 ng per sample, was labelled, purified and hybridized according to the One-Color Microarray-based Gene Expression Analysis (Quick Amp Labeling) with Tecan HS Pro Hybridization kit (Agilent Technologies). For the zebrafish, the labelled RNA was hybridized to the Zebrafish (v2) Gene Expression Microarray (Design ID 019161). For the mouse, the labelled RNA was hybridized to the SurePrint G3 Mouse GE 8×60K Microarray Design ID 028005 (Agilent Technologies). The microarrays were loaded with 1.65 µg of labelled cRNA for each postmortem sample.
2.3. Microarray calibration
Oligonucleotide (60 nt) probes on the zebrafish and mouse microarrays were calibrated using pooled labelled cRNA of all zebrafish and all mouse postmortem samples, respectively. The dilution series for the Zebrafish array was created using the following concentrations of labelled cRNA: 0.41, 0.83, 1.66, 1.66, 1.66, 3.29, 6.60 and 8.26 µg. The dilution series for the Mouse arrays was created using the following concentrations of labelled cRNA: 0.17, 0.33, 0.66, 1.32, 2.64, 5.28, 7.92 and 10.40 µg. Calibration involved plotting the signal intensities of the probes against a dilution factor and determining the isotherm model (e.g. Freundlich and/or Langmuir) that best fit the relationship between signal intensities and gene abundances.
Consider zebrafish gene transcripts targeted by A_15_P110618 (which happens to be one of the transcriptional profiles of gene Hsp70.3 shown in figure 1a). External file FishProbesParameters.txt shows that a Freundlich model best fit the dilution curve with R2 = 0.99. The equation for this probe is the following: where SI is the observed average signal intensity for the dilution x. The transcript abundance G was calculated by inverting this equation. For each probe signal intensity at postmortem time, SIt, the gene abundance G = (SIt/exp(7.1081))1/0.67632. Specifically, consider two biological replicates of 15 min postmortem zebrafish, the signal intensities of the probe A_15_P110618 are 770.5 and 576, which translates into the abundances 0.50 and 0.33 arbitrary units (arb. units), respectively. The target abundances were further converted to log10 and are shown in external file Fish_log10_AllProfiles.txt.
2.4. Statistical analysis
Abundance levels were log-transformed for analysis to stabilize the variance. A one-sided Dunnett's T-statistic was applied to test for increase at one or more postmortem times compared to live control (fish) or time 0 (mouse). A bootstrap procedure with 109 simulations was used to determine the critical value for the Dunnett's statistics in order to accommodate departures from parametric assumptions and to account for multiplicity of testing. The transcript profile for each gene was centred by subtracting the mean values at each postmortem time point to create ‘null’ profiles. Bootstrap samples of the null profiles were generated to determine the 95th percentile of the maximum (over all genes) of the Dunnett's statistics. A transcript was considered to have a significantly increased abundance when one or more points had Dunnett's T-values larger than the 95th percentile. The corresponding genes were retained for further analyses.
Orthogonal transformation of the abundances to their principal components (PCs) was conducted, and the results were graphed on a two-dimensional ordination plot. The m × n matrix of abundances (sampling times by number of gene transcripts), which is 10 × 548 for zebrafish and 7 × 515 for mouse, was used to produce an m × m matrix D of Euclidean distances between all pairs of sampling times. Principal component analysis (PCA) was performed on the matrix of distances, D. To investigate and visualize differences between the sampling times, a scatterplot of the first two principal components (PC1 and PC2) was created. To establish relative contributions of the gene transcripts, the projection of each sampling time onto the (PC1 and PC2) plane was calculated and those gene transcripts with high correlations (greater than or equal to 0.70) between abundances and either component (PC1 or PC2) were displayed as a biplot.
2.5. Gene annotation and functional categorization
Microarray probe sequences were individually annotated by performing a BLASTN search of the zebrafish and mouse NCBI databases (February 2015). The gene annotations were retained if the bit score was greater than or equal to 100 and the annotations were in the correct 5′–3′ orientation. Transcription factors, transcriptional regulators and cell signalling components (e.g. receptors, enzymes and messengers) were identified as global regulatory genes. The rest were considered response genes.
Functional categorizations were performed by querying the annotated gene transcripts in the primary literature and using UniProt (www.uniprot.org). Genes not functionally categorized to their native organism (zebrafish or mouse) were categorized to genes of phylogenetically related organisms (e.g. human). Cancer-related genes were identified using a previously constructed database (see Additional File 1: table S1 in ).
Extracting the total mRNA, calibrating the microarray probes, and determining the transcript abundances at each postmortem sampling time produced a fine-grain series of transcriptome data for the zebrafish and the mouse. Approximately 84.3% (36 811 of 43 663) zebrafish probes and 67.1% (37 368 of 55 681) mouse probes were found to provide suitable dose–response curves for calibration (electronic supplementary material, files S1–S7; http://dx.doi.org/10.5061/dryad.hv223).
Figure 2 shows the sum of all transcript abundances calculated from the calibrated probes in dependence of postmortem time. In general, the sum of all abundances decreased with time, which means that less transcript targets hybridized to the microarray probes. In the zebrafish, mRNA decreased abruptly at 12 h postmortem (figure 2a), while for the mouse brain (figure 2b), mRNA increased in the first hour and then gradually decreased. For the mouse liver, mRNA gradually decreased with postmortem time. The fact that total mRNA shown in figure 2a,b mirrors the electrophoresis patterns shown in the electronic supplementary material, figures S1 and S2 (ignoring the 28S and 18S rRNA bands) indicates a general agreement of the Gene Meter approach to the gel-based approach (i.e. Agilent Bioanalyzer). Hence, mRNA abundances depend on the organism (zebrafish, mouse), organ (brain, liver) and postmortem time, which is aligned with previous studies [11–16].
The abundance of a transcript is determined by its rate of synthesis and its rate of degradation . We focused here on transcripts that significantly increased in abundance—relative to live controls—because these genes might be actively transcribed after organismal death despite an overall decrease in total mRNA with time. A transcript was defined as having a significantly increased abundance when at least one time point was statistically higher than that of the control (figure 1a–c). It is important to understand that the entire profiles, i.e. 22 data points for the zebrafish and 20 points for the mouse, were subjected to a statistical test to determine significance (see Material and methods). We found 548 zebrafish profiles and 515 mouse profiles had significantly increased transcript abundances.
Based on GenBank gene annotations, we found that, among the transcripts with significantly increased abundances, for the zebrafish 291 were protein-coding genes (53%) and 257 non-annotated mRNA (47%), and for the mouse 324 were known protein-coding genes (63%), 190 non-annotated mRNA (37%) and one an Agilent control sequence of unknown composition. Hence, in the zebrafish and mouse, about 58% of the total genes with significant transcript abundances are known and the rest (42%) are putatively non-annotated RNA.
Examples of genes yielding transcripts that significantly increased in abundance with postmortem time are: the Heat shock protein (Hsp70.3) gene, the Thymocyte selection-associated high mobility group box 2 (Tox2) gene, and an unknown (NULL) gene (figure 1a–c). While the Hsp70.3 transcript abundance increased after 1 h postmortem to reach a maximum at 12 h, the Tox2 transcript increased after 12 h postmortem to reach a maximum at 24 h, and the NULL transcript consistently increased with postmortem time. These figures provide typical examples of transcript profiles and depict the high reproducibility of the sample replicates as well as the quality of output obtained by the Gene Meter approach.
3.1. Non-random patterns in transcript profiles
Ordination plots of the transcript profiles that had significantly increased abundances revealed prominent differences with postmortem time (figure 1d,e), suggesting the increases in transcript abundances of genes followed a discernable (non-random) pattern in both organisms. The biplots showed that 203 zebrafish transcript profiles and 226 mouse profiles significantly contributed to the ordinations. To identify patterns in the transcript profiles, we assigned them to groups based on their position in the biplots. Six profile groups were assigned for the zebrafish (A to F) and five groups (G to K) were assigned for the mouse. Determination of the average gene transcript abundances by group revealed differences in the shapes of the averaged profiles, particularly the timing and magnitude of peak abundances, which accounted for the positioning of data points in the ordinations.
Genes coding for global regulatory functions were examined separately from others (i.e. response genes). Combined results show that about 33% of the genes in the ordination plots were involved in global regulation, with 14% of these encoding transcription factors/transcriptional regulators and 19% encoding cell signalling proteins such as enzymes, messengers and receptors (electronic supplementary material, table S3). The response genes accounted for 67% of the total.
The genes were assigned to 22 categories (electronic supplementary material, File S8) with some genes having multiple categorizations. For example, the Eukaryotic translation initiation factor 3 Subunit J-B (Eif3j2) gene was assigned to protein synthesis and cancer categories .
Genes in the following functional categories were investigated: stress, immunity, inflammation, apoptosis, solute/ion/protein transport, embryonic development, epigenetic regulation and cancer. We focused on these categories because they were common to both organisms, they contained multiple genes and they might provide possible explanations for the postmortem increases in transcript abundances (e.g. epigenetic gene regulation, embryonic development, cancer). The transcriptional profiles were plotted by category and each profile was ordered by the timing of the increased abundance and peak abundances. This allowed comparisons of transcript dynamics as a function of postmortem time for both organisms. For each category, we provided the name and function of the gene and compared transcript dynamics within and between the organisms.
3.2. Stress response
In organismal death, transcripts from stress response genes were anticipated to significantly increase in abundance because these genes are activated in life to cope with perturbations, recover homeostasis  and stabilize the cytoskeleton . The stress response genes were assigned to three groups: heat shock protein (Hsp), hypoxia-related and ‘other’ responses such as oxidative stress.
In the zebrafish, Hsp gene transcripts that significantly increased in abundance included: Translocated promoter region (Tpr), Hsp70.3 and Hsp90 (figure 3). The Tpr gene encodes a protein that facilitates the export of the Hsp mRNA across the nuclear membrane  and has been implicated in chromatin organization, regulation of transcription, mitosis  and controlling cellular senescence . The Hsp70.3 and Hsp90 genes encode proteins that control the level of intracellular calcium , assist with protein folding and aid in protein degradation .
In the mouse, the Hsp gene transcripts included: Tpr, Hsp-associated methyltransferase (Mettl21) and Heat shock protein 1 (Hspe1) (figure 3). The Mettl21 gene encodes a protein modulating Hsp functions . The Hspe1 gene encodes a chaperonin protein that assists with protein folding in the mitochondria .
The timing and duration of the Hsp transcript abundances varied by organism. In general, the increase in transcript abundance of Hsp genes occurred much later in the zebrafish than the mouse (4 h versus 0.5 h postmortem, respectively). There were also differences in transcript abundance maxima of Hsp genes since they reached maxima at 9–24 h in the zebrafish, while they reached maxima at 12–24 h in the mouse. Previous studies have examined the increase of Hsp70.3 transcripts with time in live serum-stimulated human cell lines . In both the zebrafish and human cell lines (figure 1a), the Hsp70.3 gene transcript reached a maximum abundance at approximately 12 h indicating the same reactions occur in life and death.
In the zebrafish, hypoxia-related gene transcripts that significantly increased in abundance included: Carbonic anhydrase 4 (Ca4c), Nuclear factor (NF) interleukin-3 (Nfil3), Hypoxia-inducible factor 1-alpha (Hiflab) and Arginase-2 (Arg2) (figure 3). The Carbonic anhydrase 4 (Ca4c) gene encodes an enzyme that converts carbon dioxide into bicarbonate in response to anoxic conditions . The Nfil3 gene encodes a protein that suppresses hypoxia-induced apoptosis  and activates immune responses . The Hiflab gene encodes a transcription factor that prepares cells for decreased oxygen . The Arg2 gene encodes an enzyme that catalyses the conversion of arginine to urea under hypoxic conditions . Of note, the accumulation of urea presumably triggered the increase of Slc14a2 gene transcripts at 24 h, as reported in the Transport Section (below).
In the mouse, the hypoxia-related gene transcripts that significantly increased in abundance included: Methyltransferase hypoxia-inducible domain (Methig1) and Sphingolipid delta-desaturase (Degs2) (figure 3). The Methig1 gene encodes methyltransferase that presumably is involved in gene regulation . The Degs2 gene encodes a protein that acts as an oxygen sensor and regulates ceramide metabolism . Ceramides are waxy lipid molecules in cellular membranes that regulate cell-growth, death, senescence, adhesion, migration, inflammation, angiogenesis and intracellular trafficking .
The increased abundance of Ca4c transcripts in the zebrafish putatively indicates a build up of carbon dioxide 0.1–1 h postmortem in the zebrafish presumably due to lack of blood circulation. The increased abundance of the Nfil3 transcripts in the zebrafish and Methig1 transcripts in the mouse suggests hypoxic conditions exist within 0.5 h postmortem in both organisms. The increased abundance of the other hypoxia gene transcripts varied with postmortem time, with increases of Hiflab, Arg2 and Degs2 transcripts at 4 h, 12 h and 24 h, respectively.
3.2.3. Other stress responses
In the zebrafish, gene transcripts that significantly increased in abundance included: Alkaline ceramidase 3 (Acer3), Peroxirodoxin 2 (Prdx2), Immediate early (Ier2), Growth arrest and DNA-damage-inducible protein (Gadd45a), Zinc finger CCH domain-containing 12 (Zcchc12), Corticotropin-releasing hormone receptor 1 (Crhr1) and Zinc finger AN1-type domain 4 (Zfand4) (figure 3). The Acer3 gene encodes a stress sensor protein that mediates cell-growth arrest and apoptosis . The Prdx2 gene encodes an antioxidant enzyme that controls peroxide levels in cells  and triggers production of Tnfa proteins that induce inflammation . The Ier2 gene encodes a transcription factor involved in stress response . The Gadd45a gene encodes a stress protein sensor that stops the cell cycle , modulates cell death and survival, and is part of the signalling networks in immune cells . The Zcchc12 gene encodes a protein involved in stress response in the brain . The Crhr1 and Zfand4 genes encode stress proteins [44,45].
While the Acer3, Prdx2 and Ier2 transcripts increased within 0.3 h postmortem, indicating a changed physiological state, the Gadd45a transcript increased at 9 h and the other transcripts (Zcchc12, Crhr1 and Zfand4) increased at 24 h postmortem.
In the mouse, gene transcripts that significantly increased in abundance included: Membrane-associated RING-CH 4 (March4), Homocysteine-responsive endoplasmic reticulum-resident ubiquitin-like domain member 2 (Herpud2), Prohibitin-2 (Phb2), Gadd45a and Two-oxoglutarate and iron-dependent oxygenase domain-containing 1 (Ogfod1) (figure 3). The March4 gene encodes an immunologically-active stress response protein . The Herpud2 gene encodes a protein that senses the accumulation of unfolded proteins in the endoplasmic reticulum . The Phb2 gene encodes a cell surface receptor that responds to mitochondrial stress . The Ogfod1 gene encodes a stress-sensing protein .
Note that the stress gene transcripts in the mouse all increased within 1 h postmortem and remained at high abundance for 48 h.
3.2.4. Summary of stress response
In both organisms, organismal death increased the abundance of heat shock, hypoxia and ‘other stress’ gene transcripts, which varied in their timing and duration within and between organisms. Consider, for example, the Tpr and Gadd45a genes, which were common to both organisms. While the transcript abundance for the Tpr gene significantly increased within 0.5 h postmortem in both organisms, the transcript abundance for the Gadd45a gene increased at 9 h in the zebrafish and 0.5 h in the mouse. In addition, the transcriptional profile of the Tpr gene was more variable in the zebrafish than that of the mouse since the transcripts increased in abundance at 0.3 h, 9 h and 24 h postmortem, which suggests that they might be regulated through a feedback loop. By contrast, the transcriptional profile of Tpr gene in the mouse increased at 0.5 h and peaked at 12 and 24 h postmortem.
Taken together, the significant increase in transcript abundance of stress genes in both organisms is presumably to compensate for a loss of homeostasis.
3.3. Innate and adaptive immune responses
In organismal death, an increase of immune response gene transcripts was anticipated since vertebrates have evolved ways to protect the host against infection in life, even under absolutely sterile conditions . Inflammation genes were excluded from this section (even though they are innate immune genes) because we examined them in a separate section (below).
In the zebrafish, gene transcripts that significantly increased in abundance included: Early growth response-1 and -2 (Egr1, Egr2), Interleukin-1b (Il1b), l-amino acid oxidase (Laao), Interleukin-17c (Il17c), Membrane-spanning 4-domains subfamily A member 17A.1 (Ms4a17.a1), Mucin-2 (Muc2), Immunoresponsive gene 1 (Irg1), Interleukin-22 (Il22), Ubl carboxyl-terminal hydrolase 18 (Usp18), ATF-like 3 (Batf3), Cytochrome b-245 light chain (Cyba) and Thymocyte selection-associated high mobility group box protein family member 2 (Tox2) (figure 4). The Egr1 and Egr2 genes encode proteins that regulate B- and T-cell functions in adaptive immunity [51,52]. The Il1b gene encodes an interleukin that kills bacterial cells through the recruitment of other antimicrobial molecules . The Laao gene encodes an oxidase involved in innate immunity . The Il17c and Il22 genes encode interleukins that work synergically to produce antibacterial peptides . The Ms4a17.a1 gene encodes a protein involved in adaptive immunity . The Muc2 gene encodes a protein that protects the intestinal epithelium from pathogenic bacteria . The Irg1 gene encodes an enzyme that produces itaconic acid, which has antimicrobial properties . The Usp18 gene encodes a protease that plays a role in adaptive immunity . The Batf3 gene encodes a transcription factor that activates genes involved in adaptive immunity . The Cyba gene encodes an oxidase that is used to kill microorganisms . The Tox2 gene encodes a transcription factor that regulates natural killer (NK) cells of the innate immune system .
Increases of immunity gene transcripts in the zebrafish occurred at different times with varying durations. While transcripts of genes involved in adaptive immunity increased in abundance at 0.1–0.3 h (Egr), 9 h (Ms4a17.a1) and 24 h (Usp18, Batf3) postmortem, transcripts of genes involved in innate immunity increased at 4 h (Il1b), 9 h (Laao, Il17c), 12 h (Muc2, Irg1) and 24 h (Il22, Cyba,Tox2), indicating a multi-pronged and progressive approach to deal with injury and the potential of microbial invasion.
In the mouse, gene transcripts that significantly increased in abundance included: Catalytic polypeptide-like 3G (Apobec3g), CRISPR-associated endonuclease (Cas1), Perforin-1 (Prf1), Immunoglobulin heavy variable 8–11 (Ighv8-11), C4b-binding protein (C4b), Complement component C7 (C7), T-cell receptor alpha and delta chain (Tcra/Tcrd), High affinity immunoglobulin gamma Fc receptor I (Fcgr1a), Defensin (Defb30), Chemokine-4 (Ccr4), Interleukin-5 (Il5), NK cell receptor 2B4 (Cd244), Cluster of differentiation-22 (Cd22), Lymphocyte cytosolic protein 2 (Lcp2), Histocompatibility 2 O region beta locus (H2ob) and Interferon-induced transmembrane protein 1 (Ifitm1) (figure 4). The Apobec3g gene encodes a protein that plays a role in innate anti-viral immunity . The Cas1 gene encodes a protein involved in regulating the activation of immune systems [64–67]. The Prf1, C7 and Defb30 genes encode proteins that kill bacteria by forming pores in plasma membrane of target cells [68–70]. The Ighv8-11 gene encodes an immunoglobulin of uncertain function. The C4b gene encodes a protein involved in the complement system . The Tcra/Tcrd genes encode proteins that play a role in the immune response . The Fcgr1a gene encodes a protein involved in both innate and adaptive immune responses . The Ccr4 gene encodes a cytokine that attracts leucocytes to sites of infection . The Il5 gene encodes an interleukin involved in both innate and adaptive immunity [75,76]. The Cd244 and Cd22 genes encode proteins involved in innate immunity . The Lcp2 gene encodes a signal-transducing adaptor protein involved in T cell development and activation . The H2ob gene encodes a protein involved in adaptive immunity. The Ifitm1 gene encodes a protein that has anti-viral properties .
Most of the transcripts of immune response genes increased in abundance within 1 h postmortem in the mouse (n = 14 out of 16 genes), indicating a more rapid response than that of the zebrafish.
3.3.1. Summary of immune response
The increase in transcript abundance of immune response genes in both organisms included innate and adaptive immunity components. An interesting phenomenon observed in the mouse (but not zebrafish) was that four genes (C7, Tcra/Tcrd, Fcgr1a and Defb30) reached transcript abundance maxima at two different postmortem times (i.e. 1 h and 12 h) while others reached only one. The variability in the gene transcript abundances suggests possible regulation by feedback loops.
3.4. Inflammation response
The increased abundance of transcripts of inflammation genes in organismal death was anticipated because inflammation is an innate immunity response to injury. In the zebrafish, inflammation gene transcripts that increased in abundance included: Egr1, Egr2, Il1b, Tumour necrosis factor receptor (Tnfrsf19), Haem oxygenase 1 (Hmox1), Tumour necrosis factor (Tnf), G-protein receptor (Gpr31), Interleukin-8 (Il8), Tumour necrosis factor alpha (Tnfa), NF kappa B (Nfkbiaa), MAP kinase-interacting serine/threonine kinase 2b (Mknk2b) and Corticotropin-releasing factor receptor 1 (Crhr1) (figure 5). The Egr1 and Egr2 genes encode transcription factors that are pro- and anti-inflammatory, respectively [51,52,80]. The Il1b gene encodes a pro-inflammatory cytokine that plays a key role in sterile inflammation [81,82]. The Tnfrsf19 gene encodes a receptor that has pro-inflammatory functions . The Hmox1 gene encodes an enzyme that has anti-inflammatory functions and is involved in haem catabolism [84,85]. The Tnf and Tnfa genes encode pro-inflammatory proteins. The Gpr31 gene encodes a pro-inflammatory protein that activates the NF-κB signalling pathway . The Il8 gene encodes a cytokine that has pro-inflammatory properties . The Nfkbiaa gene encodes a protein that integrates multiple inflammatory signalling pathways including Tnf genes . The Mknk2b gene encodes a protein kinase that directs cellular responses and is pro-inflammatory . The Crhr1 gene modulates anti-inflammatory responses .
The increased abundance of the pro-inflammatory Egr1 transcript at 0.1 h was followed by an increase of anti-inflammatory Egr2 transcript at 0.2 h, suggesting the increase of one transcript was effecting another (figure 5). Similarly, the increased abundance of pro-inflammatory Il1b transcript at 4 h postmortem was followed by: increased abundance of pro-inflammatory Tnfrsf19, Tnf, Gpr31 and Il8 transcripts and the anti-inflammatory Hmox1 transcript at 9 h, the increased abundance of pro-inflammatory Tnfa, Nfkbiaa and Mknk2b transcripts at 12 h, and the increased abundance of anti-inflammatory Crhr1 transcripts at 24 h. Of note, while none of the pro-inflammatory gene transcripts increased in abundance past 24 h, the anti-inflammatory Crhr1 gene remained at high abundance at 48 h. It should also be noted that the Il1b, Il8 and Tnfa gene transcripts have been reported to be increased in traumatic impact injuries in postmortem tissues from human brains .
In the mouse, inflammation gene transcripts that increased in abundance included: mitogen-activated protein kinase (Map3k2), TNF receptors (Tnfrsf9, Tnfrs14), B-cell lymphoma 6 protein (Bcl6), C-C chemokine receptor-type 4 (Ccr4), Prokineticin-2 (Prok2) and platelet-activating factor receptor (Pafr) (figure 5). The Map3k2 gene encodes a kinase that activates pro-inflammatory NF-κB genes . The Tnfrsf9 and Tnfrs14 genes encode receptor proteins that have pro-inflammatory functions . The Bcl6 gene encodes a transcription factor that has anti-inflammatory functions . The Ccr4 gene encodes a cytokine receptor protein associated with inflammation . The Prok2 gene encodes a cytokine-like molecule, while the Pafr gene encodes a lipid mediator; both have pro-inflammatory functions [93,94].
Most inflammation-associated gene transcripts increased in abundance within 1 h postmortem and continued to be abundant for 12–48 h. The anti-inflammatory Bcl6 gene transcripts increased in abundance at two different times, 0.5–6 h and 24 h, suggesting that their abundances might be regulated by a feedback loop. It should also be noted that pro-inflammatory Map3k2 and Tnfrs14 gene transcripts were not at high abundance after 24 and 12 h, respectively, which also suggests regulation by a putative feedback loop from the Bcl6 transcript product.
3.4.1. Summary of inflammation response
In both organisms, some transcripts that increased in abundance have pro-inflammatory functions while others have anti-inflammatory functions. It is possible that the increases in transcript abundances are regulated by feedback loops involving an initial inflammatory reaction followed by an anti-inflammatory reaction to repress it . The variation in the gene transcript abundances of these inflammatory genes suggests the underlying regulatory network is still active in organismal death.
3.5. Apoptosis and related genes
Since apoptotic processes kill damaged cells for the benefit of the organism as a whole, we anticipated a significant increase in the abundance of apoptosis gene transcripts in organismal death.
In the zebrafish, apoptosis gene transcripts that increased in abundance included: Jun (Jdp2, Jun), Alkaline ceramidase 3 (Acer3), Fos (Fosb, Fosab, Fosl1), IAP-binding mitochondrial protein A (Diabloa), Peroxiredoxin-2 (Prdx2), Potassium voltage-gated channel member 1 (Kcnb1), Caspase apoptosis-related cysteine peptidase 3b (Casp3b), DNA-damage-inducible transcript 3 (Ddit3), BCL2 (B-cell lymphomas 2)-interacting killer (Bik) and Ras association domain family 6 (Rassf6) (figure 6). The Jdp2 gene encodes a protein that represses the activity of the transcription factor activator protein 1 (AP-1) . The Acer3 gene encodes an enzyme that maintains cell membrane integrity/function and promotes apoptosis . The Fos genes encode proteins that dimerize with Jun proteins to form part of the AP-1 that promotes apoptosis [98,99]. The Diabloa gene encodes a protein that neutralizes inhibitors of apoptosis (IAP)-binding protein  and activates caspases . The Prdx2 gene encodes antioxidant enzymes that control cytokine-induced peroxide levels and inhibit apoptosis . Although the Kcnb1 gene encodes a protein used to make ion channels, the accumulation of these proteins in the membrane promotes apoptosis via a cell signalling pathway . The Casp3b encodes a protein that plays a role in the execution phase of apoptosis . The Ddit3 gene encodes a transcription factor that promotes apoptosis. The Bik gene encodes a protein that promotes apoptosis . The Rassf6 gene encodes a protein that promotes apoptosis .
In the zebrafish, transcripts of both anti-apoptosis Jdp2 and pro-apopotosis Acer3 genes increased in abundances within 0.1 h postmortem (figure 6). These increases were followed by increases of five pro-apoptosis gene transcripts and one anti-apoptosis gene transcript within 0.3–0.5 h. The transcriptional dynamics varied among the genes. Specifically, (i) the increased abundance of the Fosb gene transcript stopped after 1 h, (ii) the transcripts of the Diabloa and Fosab genes reached abundance maxima at 0.5–4 h and then their abundances decreased after 9 h for the Diabloa and after 24 h for the Fosab genes, (iii) the Jun gene transcripts reached two maxima (one at 0.5 and another at 4–12 h)—then its abundance decreased after 24 h, and (iv) the transcript of the Prdx2 gene showed a continuous increase in abundance until reaching a maximum at 24 h and then the abundance decreased. The remaining genes were pro-apoptosis and their transcripts increased in abundance after 1–24 h postmortem. The transcripts of the Ddit3 and Rassf6 genes were very different from the other transcripts because they increased in abundance at one sampling time (12 h and 24 h, respectively) and then decreased. Apparently none of the transcripts of apoptosis genes increased in abundance after 24 h, in contrast to genes in other categories (e.g. transcripts of some of stress and immunity genes increased in abundance up to 96 h postmortem).
In the mouse, apoptosis gene transcripts that increased in abundance included: BCL2-like protein 11 (Bcl2L11), Casein kinase IIa (Csnk2a1), Interleukin 15 receptor subunit a (Il15ra), Myocyte enhancer factor 2 (Mef2a), F-box only protein 10 (Fbxo10), Sp110 nuclear body protein (Sp110), TGFB-induced factor homeobox 1 (Tgif1), Intersectin 1 (Itsm1), the Ephrin type-B receptor 3 (Ephb3) and the p21 protein-activated kinase 4 (Pak4) (figure 6). The Bcl2L11 gene encodes a protein that promotes apoptosis . The Csnk2a1 gene encodes an enzyme that phosphorylates substrates and promotes apoptosis . The Il15ra gene encodes an anti-apoptotic protein . The Mef2a gene encodes a transcription factor that prevents apoptosis . The Fbxo10 gene encodes a protein that promotes apoptosis . The Sp110 gene encodes a regulator protein that promotes apoptosis . The Tgif1 gene encodes a transcription factor that blocks signals of the transforming growth factor beta (TGFβ) pathway, and therefore is pro-apoptosis . The Itsn1 gene encodes an adaptor protein that is anti-apoptosis . The Ephb3 gene encodes a protein that binds ligands on adjacent cells for cell signalling and suppresses apoptosis . The Pak4 gene encodes a protein that delays the onset of apoptosis .
In the mouse, transcripts for the pro- and anti-apoptosis genes increased in abundance within 0.5 h postmortem; however, with the exception of Bcl2L11, most reached transcript abundance maxima at 12–48 h postmortem (figure 6). The Bcl2L11 transcripts reached abundance maxima at 1 and 6 h postmortem.
3.5.1. Summary of apoptotic response
In both organisms, transcripts of both pro- and anti-apoptosis genes increased in abundance in organismal death. However, the timings of the increases, the transcript maximum abundance and the duration of the increased abundances varied by organism. The results suggest the apoptotic genes and their regulation are distinctly different in the zebrafish than the mouse, with transcripts of the mouse genes having increased abundance to 48 h postmortem while those of the zebrafish having increased abundance to 24 h. Nonetheless, the pro- and anti-apoptosis genes appear to be inter-regulating each another.
3.6. Transport gene response
Transport processes maintain ion/solute/protein homeostasis and are involved in influx/efflux of carbohydrates, proteins, signalling molecules and nucleic acids across membranes. Transcripts of transport genes should increase in abundance in organismal death in response to dysbiosis.
In the zebrafish, transport-associated gene transcripts that increased in abundance included: Solute carrier family 26 anion exchanger member 4 (Slc26a4), Potassium channel voltage-gated subfamily H (Kcnh2), Transmembrane emp24 domain-containing protein 10 (Tmed10), Leucine-rich repeat-containing 59 (Lrrc59), the Nucleoprotein TPR (Tpr), Importin subunit beta-1 (Kpnb1), Transportin 1 (Tnpo1), Syntaxin 10 (Stx10) and Urea transporter 2 (Slc14a2) (figure 7). Of note, the four Tmed10 transcripts shown in figure 7 each represents a profile targeted by an independent probe. The transcription profiles of this gene were identical indicating high reproducibility of the Gene Meter approach. The Slc26a4 gene encodes prendrin that transports negatively charged ions (i.e. Cl−, bicarbonate) across cellular membranes . The Kcnh2 gene encodes a protein used to make potassium channels and is involved in signalling . The Tmed10 gene encodes a membrane protein involved in vesicular protein trafficking . The Lrrc59, Tpr, Tnpo1 and Kpnb1 genes encode proteins involved in trafficking across nuclear pores [118–121]. The Stx10 gene encodes a protein that facilitates vesicle fusion and intracellular trafficking of proteins to other cellular components . The Slc14a2 gene encodes a protein that transports urea out of the cell .
The transcripts of Slc26a4, Kcnh2, Lrrc59 and Tpr genes initially increased in abundance within 0.3 h postmortem and remained in high abundance for 12–24 h. The transcripts of the Tnpo1 gene increased in abundance twice, at 4 and 12 h, suggesting putative regulation by a feedback loop. The transcripts of the remaining genes increased in abundance at 24 h. The increased abundance of the Slc14a2 gene transcript suggests a build up of urea in zebrafish cells at 24–96 h postmortem, which could be due to the accumulation of urea under hypoxic conditions by the Arg2 gene (see Hsp stress response section).
In the mouse, transport-associated gene transcripts that increased in abundance included: Calcium-binding mitochondrial carrier protein (Aralar2), Sodium-coupled neutral amino acid transporter 4 (Slc38a4), SFT2 domain-containing 1 (Sft2d1), Uap56-interacting factor (Fyttd1), Solute carrier family 5 (sodium/glucose co-transporter) member 10 (Slc5a10), Mitochondrial import receptor subunit (Tom5), Translocated promoter region (Tpr), ATP-binding cassette transporter 12 (Abca12), Multidrug resistant protein 5 (Abc5), LIM and SH3 domain-containing protein (Lasp1), Chromosome 16 open reading frame 62 (C16orf62), Golgi transport 1 homologue A (Golt1a), ATP-binding cassette transporter 17 (Abca17), Nucleotide exchange factor (Sil1), Translocase of inner mitochondrial membrane 8A1 (Timm8a1), Early endosome antigen 1 (Eea1) and Potassium voltage-gated channel subfamily V member2 (Kcnv2) (figure 7). The Aralar2 gene encodes a protein that catalyses calcium-dependent exchange of cytoplasmic glutamate with mitochondrial aspartate across the mitochondrial membrane and may function in the urea cycle . The Slc38a4 gene encodes a symport that mediates transport of neutral amino acids and sodium ions . The Sft2d1 gene encodes a protein involved in transporting vesicles from the endocytic compartment of the Golgi complex . The Fyttd1 gene is responsible for mRNA export from the nucleus to the cytoplasm . The Slc5a10 gene encodes a protein that catalyses carbohydrate transport across cellular membranes . The Tom5 gene encodes a protein that plays a role in importation to proteins destined for mitochondrial sub-compartments . The Abca12, Abca17 and Abc5 genes encode proteins that transport molecules across membranes [130–132]. The Lasp1 gene encodes a protein that regulates ion transport . The C16orf62 gene encodes a protein involved in protein transport from the Golgi apparatus to the cytoplasm . The Golt1a gene encodes a vesicle transport protein . The Sil1 gene encodes a protein involved in protein translocation into the endoplasmic reticulum . The Timm8a1 gene encodes a protein that assists importation of other proteins across inner mitochondrial membranes . The Eea1 gene encodes a protein that acts as a tethering molecule for vesicular transport from the plasma membrane to the early endosomes . The Kcnv2 gene encodes a membrane protein involved in generating action potentials .
Within 0.5 h postmortem, transcripts of genes involved in: (i) ion and urea regulation (Aralar), (ii) amino acid (Slc38a4), carbohydrate (Slc5a10) and protein (Sft2d1, Tom5) transport, (iii) mRNA nuclear export (Fyttd1, Tpr) and (iv) molecular efflux (Abca12, Abc5) increased in abundance in the mouse. The transcription profiles of these genes varied in terms of transcript abundance maxima and duration. While the transcripts of Aralar, Sft2d1, Slc38a4, Fyttd1 and Slc5a10 reached abundance maxima at 1 h, those of Tom5, Tpr, Abca12 and Abc5 reached maxima at 12–24 h postmortem. The duration of the increased abundance also varied for these transcripts since most remained at high abundances for 48 h postmortem, while the Sft2d1, Fyttd1 and Slc5a10 transcripts were at high abundances from 0.5 to 12+ h. The shorter duration of increased abundance suggests prompt gene repression. The transcript abundances of Lasp1, C16orf62, Golt1a and Abca17 increased at 1 h postmortem and remained elevated for 48 h. The transcripts of Sil1, Timm8a1 and Eea1 increased in abundance at 6 h, while those of Kcnv2 increased at 24 h postmortem and remained elevated for 48 h.
3.6.1. Summary of transport genes
The increased abundance of transcripts of transport genes suggests attempts by zebrafish and mice to reestablish homeostasis. Although the transcripts of half of these genes increased in abundance within 0.5 h postmortem, many increased at different times and for varying durations. While most of the transcripts of transport genes in the zebrafish were not abundant after 24 h, most transcripts of transport genes in the mouse remained abundant for 24–48 h postmortem.
3.7. Developmental control genes
An unexpected finding in this study was the increased abundance of transcripts of developmental control genes in organismal death. Developmental control genes are mostly involved in regulating developmental processes from early embryo to adult in the zebrafish and mouse; therefore, we did not anticipate their transcripts to become more abundant in organismal death.
In the zebrafish, development-associated gene transcripts that increased in abundance included: LIM domain-containing protein 2 (Limd2), Disheveled-associated activator of morphogenesis 1 (Daam1b), Meltrin alpha (Adam12), Hatching enzyme 1a (He1a), Midnolin (Midn), Immediate early response 2 (Ier2), Claudin b (Cldnb), Regulator of G-protein signalling 4-like (Rgs4), Proline-rich transmembrane protein 4 (Prrt4), Inhibin (Inhbaa), Wnt inhibitory factor 1 precursor (Wif1), Opioid growth factor receptor (Ogfr), Strawberry notch homolog 2 (Sbno2) and Developing brain homeobox 2 (Dbx2) (figure 8). The Limd2 gene encodes a binding protein that plays a role in zebrafish embryogenesis . The Daam1b gene regulates endocytosis during notochord development . The Adam12 gene encodes a metalloprotease-disintegrin involved in myogenesis . The He1a gene encodes a protein involved in egg envelope digestion . The Midn gene encodes a nucleolar protein expressed in the brain that is involved in the regulation of neurogenesis [143,144]. The Ier2 gene encodes a protein involved in left–right asymmetry patterning in the zebrafish embryo . The Cldnb gene encodes a tight junction protein in larval zebrafish . The Rgs4 gene encodes a protein involved in brain development . The Prrt4 gene encodes a protein that is predominantly expressed in the brain and spinal cord in embryonic and postnatal stages of development. The Inhbaa gene encodes a protein that plays a role in oocyte maturation . The Wif1 gene encodes a WNT inhibitory factor that controls embryonic development . The Ogfr gene plays a role in embryonic development . The Sbno2 gene plays a role in zebrafish embryogenesis . The Dbx2 gene encodes a transcription factor that plays role in spinal cord development .
Although the abundances of Limd2, Daam1, Adam12 and He1a transcripts increased in the zebrafish within 0.1 h postmortem, other gene transcripts in this category increased from 0.3 to 24 h postmortem, reaching abundance maxima at 24 h or more.
In the mouse, development-associated gene transcripts that increased in abundance included: MDS1 and EVI1 complex locus protein EVI1 (Mecom), MAM domain-containing glycosylphosphatidylinositol anchor 2 (Mdga2), FYVE, RhoGEF and PH domain-containing 5 (Fgd5), RNA-binding motif protein 19 (Rbm19), Chicken ovalbumin upstream promoter (Coup), Single-minded homolog 2 (Sim2), Solute carrier family 38, member 4 (Slc38a4), B-cell lymphoma 6 protein (Bcl6), Sema domain transmembrane domain (TM) cytoplasmic domain (semaphorin) 6D (Sema6d), RNA binding motif protein 45 (Rbm45), Transcription factor E2F4 (E2f4), Long chain fatty acid-CoA ligase 4 (Lacs4), Kallikrein 1-related peptidase b3 (Klk1b3), Sema domain, immunoglobulin domain, TM and short cytoplasmic domain (Sema4c), TGFB-induced factor homeobox 1 (Tgif1), Interferon regulatory factor 2-binding protein-like (Irf2bpl), Ephrin type-B receptor 3 (Ephb3), Testis-specific Y-encoded-like protein 3 (Tspyl3), Protein ripply 3 (Ripply3), Src kinase-associated phosphoprotein 2 (Skap2), DNA polymerase zeta catalytic subunit (Rev3l), MKL/myocardin-like 2 (Mkl2) and Protein phosphatase 2 regulatory subunit A (Ppp2r1a) (figure 8). The Mecom gene plays a role in embryogenesis and development . The Mdga2 gene encodes immunoglobins involved in neural development . The Fgd5 gene is needed for embryonic development since it interacts with hematopoietic stem cells . The Rbm19 gene is essential for preimplantation development . The Coup gene encodes a transcription factor that regulates development of the eye  and other tissues . The Sim2 gene encodes a transcription factor that regulates cell fate during midline development . The Slc38a4 gene encodes a regulator of protein synthesis during liver development and plays a crucial role in fetal growth and development [160,161]. The Bcl6 gene encodes a transcription factor that controls neurogenesis . The Sema6d gene encodes a protein involved in retinal development . The Rbm45 gene encodes a protein that has preferential binding to poly(C) RNA and is expressed during brain development . The E2f4 gene is involved in maturation of cells in tissues . The Lacs4 gene plays a role in patterning in embryos . The Klk1b3 gene encodes a protein that plays a role in the developing embryos . The Sema4c gene encodes a protein that has diverse function in neuronal development and heart morphogenesis [168,169]. The Tgif1 gene encodes a transcription factor that plays a role in trophoblast differentiation . The Irf2bpl gene encodes a transcriptional regulator that plays a role in female neuroendocrine reproduction . The Ephb3 gene encodes a kinase that plays a role in neural development . The Tspyl3 gene plays a role in testis development . The Ripply3 gene encodes a transcription factor involved in development of the ectoderm . The Skap2 gene encodes a protein involved in actin reorganization in lens development . The Rev3l gene encodes a polymerase that can replicate past certain types of DNA lesions and is necessary for embryonic development . The Mkl2 gene encodes a transcriptional co-activator that is involved in the formation of muscular tissue during embryonic development . The Ppp2r1a gene plays a role in embryonic epidermal development .
The transcripts of Mecom, Mdga2, Fgd5, Rbm19, Coup, Sim2, Slc38a4, Bcl6, Sema6d, Rbm45, E2f4 and Lacs4 genes in the mouse significantly increased in abundance within 0.5 h postmortem but the other transcripts increased from 1 h to 48 h reaching abundance maxima at 12 h or more.
3.7.1. Summary of developmental control genes
In organismal death, there is progressive increase in transcript abundances of some developmental control genes suggesting that they are no longer silenced. A possible reason for these increased abundances is that the postmortem physiological conditions resemble those of earlier developmental stages.
3.8. Cancer genes
There are a number of databases devoted to cancer and cancer-related genes. Upon cross-referencing the genes found in this study, we discovered a significant overlap. The genes found in this search are presented below.
In the zebrafish, transcripts of the following cancer genes significantly increased in abundance: Jdp2, Xanthine dehydrogenase (Xdh), Egr1, Adam12, Myosin-IIIa (Myo3a), Fosb, Jun, Integrin alpha 6b (Itga6), Ier2, Tpr, Dual specificity protein phosphatase 2 (Dusp2), Disintegrin and metallopeptidase domain 28 (Adam28), Tnpo1, Ral guanine nucleotide dissociation stimulator-like (Rgl1), Carcinoembryonic antigen-related cell adhesion molecule 5 (Ceacam1), Fosl1, Il1b, Hif1a, Serine/threonine-protein phosphatase 2A regulatory (Ppp2r5d), DNA replication licensing factor (Mcm5), Gadd45, Myosin-9 (Myh9), Casp3, Tnf, Il8, Cyclic AMP-dependent transcription factor (Atf3), small GTPase (RhoA), Mknk2, Ephrin type-A receptor 7 precursor (Epha7), ETS-related transcription factor (Elf3), Nfkbia, Kpnb1, Wif1, RAS guanyl-releasing protein 1 (Rasgrp), Ras association domain-containing protein 6 (Rassf6), Cyba, DNA-damage-inducible transcript 3 (Ddit3), Serine/threonine-protein kinase (Sbk1) and Tyrosine-protein kinase transmembrane receptor (Ror1) (figure 9).
In the mouse, transcripts of the following cancer genes significantly increased in abundance: Retinoblastoma-like protein 1 (Rbl1), Elongation factor RNA polymerase II (Ell), Bcl-2-like protein 11 (Bcl2l11), Sal-like protein 1 (Sall1), Map3k2, Bcl6, Tnfrsf9, CK2 target protein 2 (Csnk2a1), Transcription factor E2f4 (E2f4), Zinc finger DHHC-type containing 14 (Zdhhc14), Tpr, RAS p21 protein activator 1 (Rasa1), Gadd45, Prohibitin (Phb2), Serine/threonine-protein phosphatase PP1-gamma catalytic (Ppp1cc), Lasp1, G protein-coupled receptor kinase 4 (Grk4), LIM domain transcription factor (Lmo4), Protein phosphatase 1E (Ppm1e), Protein sprouty homolog 1 (Spry1), Multiple PDZ domain protein (Mpdz), Kisspeptin receptor (Kiss1), Receptor-type tyrosine-protein phosphatase delta precursor (Ptprd), Small effector protein 2-like (Cdc42), AT-rich interactive domain-containing protein 1A (Arid1a), Lymphocyte cytosolic protein 2 (Lcp2), DNA polymerase zeta catalytic subunit (Rev3l), Tnfrsf14, Integrin beta-6 precursor (Itgb6), Triple functional domain protein (Trio), ATPase class VI type 11C (Atp11c) and Serine/threonine-protein phosphatase 2A regulatory (Ppp2r1a) (figure 9).
3.8.1. Summary of cancer genes
Genes analysed under this category were classified as ‘cancer genes’ in a Cancer Gene Database  (figure 9). The timing, duration and peak transcript abundances differed within and between organisms. Note that some gene transcripts had two abundance maxima. In the zebrafish, this phenomenon occurred for Adam12, Jun, Tpr, Dusp2, Tnpo1 and Hif1a genes and in the mouse, Bcl6, Tnfrs9, Lasp1, Cdc42 and Lcp2 genes, and is consistent with the notion that the transcript abundances are being regulated through feedback loops.
3.9. Epigenetic regulatory genes
Epigenetic regulation of gene expression involves DNA methylation and histone modifications of chromatin into active and silenced states . These modifications alter the condensation of the chromatin and affect the accessibility of the DNA to the transcriptional machinery. Although epigenetic regulation plays an important role in development, modifications can arise stochastically with age or in response to environmental stimuli . Hence, we anticipated that epigenetic regulatory genes would be involved in organismal death.
In the zebrafish, transcripts of the following epigenetic genes significantly increased in abundance: Jun dimerization protein 2 (Jdp2), Chromatin helicase protein 3 (Chd3), Glutamate-rich WD repeat-containing protein 1 (Grwd1), Histone H1 (Histh1l), Histone cluster 1, H4-like (Hist1h46l3) and Chromobox homolog 7a (Cbx7a) (figure 10). The Jdp2 gene is thought to inhibit the acetylation of histones and repress expression of the c-Jun gene . The Chd3 gene encodes a component of a histone deacetylase complex that participates in the remodelling of chromatin . The Grwd1 gene is thought to be a histone-binding protein that regulates chromatin dynamics at the replication origin . The Histh1l gene encodes a histone protein that binds the nucleosome at the entry and exit sites of the DNA and the Hist1h46l3 gene encodes a histone protein that is part of the nucleosome core . The Cbx7a gene encodes an epigenetic regulator protein that binds non-coding RNA and histones and represses gene expression of a tumor suppressor .
The transcripts of both Jdp2 and Chd3 genes increased in abundance within 0.3 h postmortem, and reached abundance maxima at 0.5 h. Note that two different probes targeted the Jdp2 transcript. The transcript of the Grwd1 gene increased in abundance at 1 h and 24 h postmortem. The transcript of the histone genes increased in abundance at 4 h postmortem, reaching abundance maxima at 24 h. The transcript of the Cbx7a gene increased in abundance at 12 h, reaching an abundance maximum at 24 h. The transcript abundances of these genes decreased after 24 h.
In the mouse, transcripts of the following epigenetic genes significantly increased in abundance: Tubulin tyrosine ligase-like family member 10 (Ttll10), Histone cluster 1 H3f (Hist1h3f), Histone cluster 1 H4c (Hist1h4c), YEATS domain-containing 2 (Yeats2), Histone acetyltransferase (Kat7) and Probable JmjC domain-containing histone demethylation protein 2C (Jmjd1c) (figure 10). The Ttll10 gene encodes a polyglycylase involved in modifying nucleosome assembly protein 1 that affects transcriptional activity, histone replacement and chromatin remodelling . The Hist1h3f and Hist1h4c genes encode histone proteins that are the core of the nucleosomes . The Yeats2 gene encodes a protein that recognizes histone acetylations so that it can regulate gene expression in the chromatin . The Kat7 gene encodes an acetyltransferase that is a component of histone-binding origin-of-replication complex, which acetylates chromatin and therefore regulates DNA replication and gene expression . The Jmjd1c gene encodes an enzyme that specifically demethylates Lys-9 of histone H3 and is implicated in the reactivation of silenced genes .
The transcripts of the Ttll10, Yeats2 and histone protein genes increased in abundance 0.5 h postmortem and reached maxima at different times, with the Ttll10 transcript reaching a maximum at 1 to 6 h, the histone transcripts reaching maxima at 6 and 12 h postmortem, and the Yeats2 transcript reaching maxima at 12–24 h postmortem (figure 10). The transcripts of the Kat7 and Jmjd1c genes increased in abundance at 24 h, reaching abundance maxima at 48 h postmortem. Note that the transcripts of the histone genes were no longer abundant after 24 h postmortem.
3.9.1. Summary of epigenetic regulatory genes
The increased abundance of transcripts of genes encoding histone proteins, histone–chromatin modifying proteins, and proteins involved in regulating DNA replication at the origin were common to the zebrafish and the mouse. These findings indicate that epigenetic regulatory genes are still modifying chromatin structure in organismal death and thus change the accessibility of transcription factors to the promoter or enhancer regions.
3.10. Percentage of gene transcripts with significant abundance by postmortem time
The percentage of gene transcripts was defined as the number of gene transcripts with abundances greater than the control over the total number of transcripts with significant abundance in a category at a specific postmortem time. A comparison of the percentage of gene transcripts by postmortem time of all gene categories revealed similarities between the zebrafish and the mouse. Specifically, most gene transcripts increased in abundance between 0.5 and 24 h postmortem, and after 24 h the transcript abundance drastically dropped (figure 11, ‘All genes’). It should be noted that the same pattern was found in stress, transport and development categories for both organisms. However, in the zebrafish, the immunity, inflammation, apoptosis and cancer categories differed from the mouse. Specifically, the gene transcripts in the immunity, inflammation and cancer categories increased in abundance much later (1–4 h) in the zebrafish than the mouse, and the duration of elevated abundances was much shorter. For example, while 90% of the transcripts for genes in the immunity and inflammation categories increased in abundance in the mouse within 1 h postmortem, less than 30% of the transcripts in the same categories were abundant in the zebrafish (figure 11), indicating a slower initial response. It should be noted that while the number of transcripts of immunity genes reached abundance maxima at 24 h postmortem in both organisms, the number of inflammation genes reaching abundance maxima occurred at 1–4 h in the mouse and 24 h in the zebrafish. The significance of these results is that the inflammation response occurs rapidly and robustly in the mouse while in the zebrafish it takes longer to establish, which could be attributed to phylogenetic differences. There were significant differences in the transcript abundances of apoptosis genes between the zebrafish and the mouse. In the mouse, the percentage of transcripts of apoptosis genes reached 100% at 1 h postmortem and remained sustained for 48 h postmortem, while the percentage of transcript genes with increased abundance in the zebrafish never reached 70% and the abundances abruptly decreased after 12 h.
3.11. Upregulation or differential mRNA stability?
Since equal amounts of RNA were used for all time points (see below), although degradation was ongoing, it is theoretically possible that the apparent increase in the abundance of a subset of transcripts is actually due to a higher stability of these transcripts compared to the background of degrading transcripts. Hence, the question arises whether higher transcript abundances are due to upregulation after organismal death or complex decay profiles leading to relative enrichments. To determine whether the significant increases were due to such an enrichment, the expected profile of a hypothetical stable non-degrading cRNA was calculated for the zebrafish, mouse liver and mouse brain. In theory, the abundance of a stable non-degrading cRNA transcript determined by the Gene Meter approach should positively correlate to the amount of total cRNA delivered to the DNA microarray. Below is a rationale and the approach used to identify stable non-degrading cRNAs in the transcript pool.
As outlined in the Material and methods section, the amount of sample taken from an animal was approximately the same and the homogenization volume was the same. The electronic supplementary material, tables S1 and S2 show the quantity of total RNA extracted from a tissue (x, ng µl−1). Since a fixed amount of RNA was taken into labelling, the volume of the homogenized sample was proportional to 1/x, i.e. the effective quantity of tissue taken into the microarray analysis was proportional to 1/x.
Let us assume there was a subset of stable RNAs, while all other RNA molecules were degrading. Hence, the quantity of the stable gene transcript would be directly proportional to the amount of tissue taken into the experiment, 1/x. In order to provide an expected concentration–time profile for the assumed stable cRNA, one can compare (on the log2 scale) the 1/x values for all time points. To make it relative to the live control, the log2 value of the control will be subtracted from each time point. The obtained profile will be the expected profile (fold change) of a stable non-degrading cRNA.
Taking into account the above considerations, we found that the expected profile of the stable zebrafish cRNA would have an eightfold increase at 96 h postmortem (figure 12). By contrast, at 48 h postmortem, we found that the expected profile of the stable cRNA in mouse liver would have an approximately fourfold decrease, while the stable cRNA in the mouse brain would have an approximately twofold decrease, because less tissue was taken into the microarray analysis (since the total RNA yield increased). It is important to note that these expected stable mouse mRNA profiles (lower two panels of figure 12) would not have been selected by the statistical procedure for identifying transcriptional profiles that significantly increased in abundance relative to the live controls.
In the zebrafish, the potentially enriched transcripts (due to their stability) were identified by correlating their abundances to the expected fold change of the stable cRNA. In theory, potentially enriched transcripts should be positively correlated to the expected fold change of the putatively stable transcript. Alternatively, transcripts that are not enriched due to stability effects should be negatively correlated or not correlated at all to the expected fold change. The correlations of the 548 gene transcripts (i.e. those that were significantly increased in abundance with postmortem time) ranged from highly negative to highly positive (figure 13). Note higher frequency gene transcripts on the left side of the histogram indicate that most transcripts were negatively correlated to the expected fold change, while the lower frequency on the right side of the histogram indicates a small portion of the gene transcripts were putatively enriched than other transcripts. A standard statistical table, using d.f. = n − 2 with direction, revealed that correlations greater than 0.685 were statistically significant at α = 0.01 (three bars on the right side of the histogram). Hence, 45 of the 548 gene transcripts were found to be significantly correlated to the expected fold change of the stable cRNA.
The 45 gene transcripts that were putatively enriched are shown in table 1. Figure 14 shows the transcriptional profiles of three selected gene transcripts compared to expected fold change profile as shown in the top three panels. The transcriptional profile of a negatively correlated gene transcript served as a control. None of the gene transcripts from the mouse samples were enriched because the amount of cRNA in the tissue extract increased or stayed about the same with postmortem time.
The primary motivation for our study was driven by curiosity in the processes involved in the shutting down of a complex biological system—which has received little attention so far. Other fields of research have examined the shutdown of complex systems (e.g. societies , government  and electrical black outs ). Yet, to our knowledge, no study has examined long-term postmortem dynamics of transcripts from vertebrates kept in their native conditions. The secondary motivation for our study was to demonstrate the utility of Gene Meter technology for gene expression studies.
4.1. Why study transcriptional dynamics in death?
While the development of a complex biological system requires time and energy, its shutdown and subsequent disassembly entails the dissipation of energy and the unravelling of complex structures and could provide novel insights into interesting pathways. Since the shutdown of a complex system does not occur instantaneously, not all cells in a body will be immediately affected by pending death. Hence, the transcript pools should contain mRNAs involved in day-to-day survival as well as stress compensation—yet, since death is the unequivocal end of life—one would expect mRNAs will change with postmortem time. How the transcription pools dynamically change with postmortem time was the goal of this study.
As one would expect, a living system is a collection of biochemical reactions linked together by the components participating in them. These reactions depending on one another to a certain extent, we conjecture that the observed increases in transcript abundance are due to thermodynamic and kinetic conditions that are encountered during organismal death and that reflect, at least in part, pathways that are used during life. For example, the increased abundances of epigenetic regulatory transcripts suggest that histone modification (e.g. Histh1l) and chromatin interactions (e.g. Grwd1, Chd3, Yeats, Jmjd1c) could be taking place (figure 10). Products of these transcripts could be responsible for the unravelling of the nucleosomes, which enable transcription factors and RNA polymerases to transcribe the developmental control genes previously silenced since embryogenesis (figure 8). Hence, one plausible explanation for the increase in transcript abundances is that specific genes are transcriptionally upregulated. The energy barrier in this example is the tightly wrapped nucleosomes that previously did not allow access to developmental control genes. Other energy or entropy barriers include the nucleopores that allow the exchange of mRNA and other molecules between the mitochondria and the cytosol (e.g. Tpr, Tnpo1, Lrrc59), or the ion/solute protein channels (e.g. Aralar2, Slc38a4) that control intracellular ions regulating apoptotic pathways [194,195].
4.2. Methodological validity
The Gene Meter approach is pertinent to the quality of the microarray output obtained in this study because conventional DNA microarrays yield noisy data [196,197]. The Gene Meter approach determines the behaviour of every microarray probe by calibration—which is analogous to calibrating a pH meter with buffers. Without calibration, the precision and accuracy of a meter is not known, nor can one know how well the experimental data fits to the calibration (i.e. R2). In the Gene Meter approach, the response of a probe (i.e. its behaviour in a dilution series) is fitted to either a Freundlich or Langmuir adsorption model, and probe-specific parameters are calculated. The ‘noisy’ or ‘insensitive’ probes are identified and removed from further analyses. Probes that sufficiently fit the model are retained and later used to calculate the abundance of a specific gene or gene transcript in a biological sample. The models take into consideration the nonlinearity of the microarray signal and the calibrated probes do not require normalization procedures to compare biological samples. By contrast, conventional DNA microarray approaches are biased because different normalizations can yield up to 20–30% differences in the up- or downregulation depending on the procedure selected [198–201]. Another issue with normalization is that it will artificially increase some of the transcripts due to the denominator. For example, if one normalizes gene abundances to the sum of all transcripts (as often done in conventional microarray studies) and the majority of transcripts degrades with postmortem time, the normalized values of the stable transcripts will be artificially increased due to the decreasing denominator.
We recognize that next-generation sequencing (NGS) approaches could have been used to monitor the increase in transcript abundances in this study. However, the same problems of normalization and reproducibility (mentioned above) are pertinent to NGS technology . Hence, the Gene Meter approach is currently the most advantageous to study postmortem gene transcript abundances in a high-throughput manner. Moreover, a recent publication by a group from the US National Institute of Standards and Technology used the same dilution series approach (as we did in this study) to evaluate and calibrate RNASeq . They found RNASeq comparable to microarrays in terms of target quantification, but not superior, as may be perceived by the community.
4.3. Transcription versus decay
Other plausible explanations for the observed increase in transcript abundances include enrichment of stable non-degrading RNA in the transcription pool  and/or the changing of cell types in the samples with postmortem time. In this study, we specifically addressed the enrichment issue by identifying 45 gene transcripts in the zebrafish that could have been artificially ‘enriched’ after 12 h postmortem (table 1 and figure 14). In contrast to the zebrafish, we found none of the gene transcripts were enriched in the mouse. The significance of this finding is that a very low percentage (less than 4%) of the 1063 transcriptional profiles examined in our study could be classified as ‘artificially’ enriched due to differential stability. This inference is based on assuming constant decay rates, but formally we cannot exclude the possibility that complex differential stability effects might also occur, which by themselves might be parts of regulatory loops. Hence, it will be of interest in future experiments to assess with more direct means the fraction of genes that get actively transcribed after death versus those that show complex decay patterns.
It should be noted that a small number of genes has been previously reported to be upregulated in cadavers. Using reverse transcription real-time quantitative PCR (RT-RTqPCR), a study showed significant increases in expression of Myosin light chain 3 (Myl3), Matrix metalloprotease 9 (Mmp9), and Vascular endothelial growth factor A (Vegfa) genes in body fluids after 12 h postmortem . Interestingly, we found an increased abundance of myosin-related and matrix metalloprotease transcripts in our study. Specifically, the myosin-related genes included: Myosin-Ig (Myo1g) in the mouse, and Myosin-IIIa (Myo3a) and Myosin-9 (Myh9) in the zebrafish. The matrix metalloproteinase genes included the Metalloproteinase-14 (Mmp14b) gene in the zebrafish. The Myo1g gene encodes a protein regulating immune response , the Myo3a gene encodes an uncharacterized protein, the Myh9 gene encodes a protein involved in embryonic development , and the Mmp14b gene encodes an enzyme regulating cell migration during zebrafish gastrulation . The Myo1g, Myh9 and Mmp14b transcripts began to increase right after death and reached abundance maxima at 24 h postmortem, while the Myo3a transcript reached an abundance maximum at 12 h postmortem. The significance of these results is twofold: (i) two different technologies (RT-RTqPCR and Gene Meter) have now demonstrated increased transcript abundances and these increases have now been reported in three organisms (human, zebrafish and mouse); and (ii) there might be significant overlap in gene transcripts that increase in abundance in death as we have showed with myosin- and matrix metalloprotease genes, which warrants further studies using other vertebrates. The purpose of such studies would be to understand common mechanisms involved in the shutdown of highly ordered biological systems.
4.4. Stability of cell types
Changes of the cell types in the samples with postmortem times could also account for the increases in the transcript abundances with postmortem time because of the differential survival of various cell types. In the case of human blood cells, for example, eosinophils, monocytes, neutrophils and lymphocytes were present just after death, but at 60 h postmortem eosinophils and monocytes were not found, at 66 h postmortem neutrophils were not found and at 86 h postmortem lymphocytes were not found . Skeletal muscle stem cells from mouse, for example, adopt a dormant cell state and retain regenerative capacity for 14–17 days after organismal death . Fibroblast cells from sheep can be cultured for 56 h , and fibroblast cells from goats can be cultured for 41–160 days postmortem [212,213]. Similarly, inner ear stem cells from mice can be cultured after 5–10 days postmortem . Taken together, some cells are more resilient than others and are the last ones to die. If stem cells are the last ones to die, then the global transcriptome will become more ‘stem-like’ with postmortem time.
If the cellular composition of the organ/tissues did change between sampling times, then this could contribute to the increases in transcript abundances. Our study analysed whole zebrafish and dissected organs of the mice that contain multiple cell types, which could dilute the transcriptional contribution of any one specific cell type. Current limitations prevent us from using cell-type-specific approaches in systematically analysing the postmortem transcriptome. Furthermore, the number of zebrafish and mouse samples analysed so far is not sufficient to investigate the full magnitude of transcriptional changes that could be occurring.
4.5. What do the increases in postmortem transcript abundances mean in the context of life?
Since increases in postmortem transcript abundances occurred in both the zebrafish and the mouse in our study, it is reasonable to suggest that other multicellular eukaryotes will display a similar phenomenon. What does this phenomenon mean in the context of organismal life? We conjecture that the highly ordered structure of an organism—evolved and refined through natural selection and self-organizing processes —undergoes a thermodynamically driven process of spontaneous disintegration through complex pathways, which apparently involve the increased abundance of specific gene transcripts and putative feedback loops. While evolution played a role in pre-patterning of these pathways, it probably does not play any role in its disintegration fate. However, one could argue that some of these pathways have evolved to favour healing or ‘resuscitation’ after severe injury, which would be a possible adaptive advantage. The increased abundance of inflammation response transcripts, for example, putatively indicates that a signal of infection or injury is sensed by the still alive cells after death of the body. Alternatively, these increases could be due to fast decay of some repressors of genes or whole pathways leading to the transcription of genes. Hence, it will be of interest to study this in more detail, since this could, for example, provide insights into how to better preserve organs retrieved for transplantation.
We observed clear qualitative and quantitative differences between two organs (liver and brain) in the mouse in their degradation profiles (figure 2). We also showed an increase in transcript abundances for immunity, inflammation and cancer genes within 1 h of death (figure 11). It would be interesting to explore if these differences are comparable to what occurs in humans, and we wonder how much of the transplant success could be attributed to differences in the synchronicity of postmortem expression profiles rather than immunosuppression agents [216,217]. Our study provides an alternative perspective to the fate of transplant recipients due to the increase of transcripts of regulatory and response genes after the sample has been harvested from the donor.
This is the first comprehensive study to assess changes in transcriptomic profiles after organismal death and raises interesting questions relative to transplantology, inflammation, cancer, evolution and molecular biology.
The euthanasia methods outlined above are approved by the American Veterinary Medical Association (AVMA) Guidelines for Euthanasia (www.avma.org) and carried out by personnel of the Max Planck Institute for Evolutionary Biology (Ploen, Germany). All animal work followed the legal requirements, was registered under no. V312-72241.123-34 (97-8/07) and approved by the ethics commission of the Ministerium für Landwirtschaft, Umwelt und ländliche Räume, Kiel (Germany) on 27 December 2007.
The datasets supporting this article have been uploaded as part of the electronic supplementary material (see additional files) and are available at the Dryad Digital Repository, http://dx.doi.org/10.5061/dryad.hv223 .
D.T., R.N., T.D.L. and A.E.P. designed the study. A.E.P., R.N. and T.D.L. carried out the molecular genetic studies. B.G.L. and A.E.P. identified gene transcript profiles with statistically significant increased abundances. S.S. and P.A.N. annotated the genes. A.E.P. and P.A.N. conducted the statistical analyses. A.E.P., D.T. and P.A.N. wrote the manuscript. All authors read and approved the final manuscript.
The authors declare no competing interests.
The work was supported by funds from the National Cancer Institute P20 Partnership grant nos (P20 CA192973 and P20 CA192976) for S.S. and the Max Planck Society for D.T., R.N., T.D.L. and A.E.P.
We thank Till Sckerl for help with sacrificing the mice, Nicole Thomsen for help with tissue and RNA processing, and Elke Blohm-Sievers for microarray work. We thank Russell Bush for his helpful advice on thermodynamics and entropic barriers.
Electronic supplementary material is available online at https://dx.doi.org/10.6084/m9.figshare.c.3653192.
- Received September 24, 2016.
- Accepted December 12, 2016.
- © 2017 The Authors.
Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.