Skip to main content
Download PDF

The forgotten variable? Does the euthanasia method and sample storage condition influence an organisms transcriptome – a gene expression analysis on multiple tissues in pigs

BMC Genomics volume 24, Article number: 769 (2023) Cite this article

Abstract

Background

Transcriptomic studies often require collection of fresh tissues post euthanasia. The chosen euthanasia method might have the potential to induce variations in gene expressions that are unlinked with the experimental design. The present study compared the suitability of ‘nitrogen gas in foam’ (ANOXIA) in comparison to a non-barbiturate anaesthetic, T-61® (T61), for euthanizing piglets used in transcriptome research. Further, the effect of common tissue storage conditions, RNAlater™ (RL) and snap freezing in liquid nitrogen (LN2), on gene expression profiles were also analysed.

Results

On comparison of the 3’mRNA-Seq data generated from pituitary, hypothalamus, liver and lung tissues, no significant differential expression in the protein coding genes were detected between the euthanasia methods. This implies that the nitrogen anoxia method could be a suitable alternative for euthanasia of piglets used in transcriptomic research. However, small nuclear RNAs (snRNAs) that constitute the eukaryotic spliceosomal machinery were found to be significantly higher (log2fold change ≥ 2.0, and adjusted p value ≤ 0.1) in pituitary samples collected using ANOXIA. Non-protein coding genes like snRNAs that play an important role in pre-mRNA splicing can subsequently modify gene expression. Storage in RL was found to be superior in preserving RNA compared to LN2 storage, as evidenced by the significantly higher RIN values in representative samples. However, storage in RL as opposed to LN2, also influenced differential gene expression in multiple tissues, perhaps as a result of its inability to inhibit biological activity during storage. Hence such external sources of variations should be carefully considered before arriving at research conclusions.

Conclusions

Source of biological variations like euthanasia method and storage condition can confound research findings. Even if we are unable to prevent the effect of these external factors, it will be useful to identify the impact of these variables on the parameter under observation and thereby prevent misinterpretation of our results.

Peer Review reports

Background

Transcriptome studies in pigs have increased over the past decade, aiming to interpret tissue specific processes that control significant traits related to growth, immunity and disease development [1,2,3,4]. Due to the closeness to human physiology and anatomy as well as the genomic similarities, pigs have been considered as a suitable model for translational research [5,6,7,8]. Several studies have been using pigs as animal models for preclinical studies in humans [9, 10].

Transcriptome studies often require collection of tissues post euthanasia. The recommendations by the European Council [11, 12] as well as the guidelines from the American Veterinary Medical Association [13] enlist several criteria on the choice of euthanasia methods in experimental animals. For the selection of euthanasia method, prime consideration must be given to ensure the welfare of the animals. Additionally, the chosen method should not influence subsequent evaluation of the parameter under investigation. The euthanasia method should keep the tissue of interest intact and more importantly, it should not cause a source of variation in the gene expression profile. Alike euthanasia methods, tissue storage condition has also got the potential to induce variations in gene expressions, that are unlinked with the experimental design [14,15,16,17]. This can further introduce a bias in the top gene expression rankings and subsequently mask the experimental effects. Hence, in the context of transcriptomic research, it is important that the (m)RNA levels related to the tissue under investigation remain unlinked to the method of euthanasia and tissue storage conditions used in the experiment.

Recently, we examined the effect of an inhalant anaesthetic ‘nitrogen gas in foam’ (ANOXIA) on RNA yield and quality parameter in piglets and compared it against an injectable anaesthetic T-61® (T61)[18]. We combined it with two different storage conditions—RNAlater™ (RL) vs snap freezing in liquid nitrogen (LN2). It was found that ANOXIA could be used as a suitable alternative to T61 method based on RNA quality parameters, while storage in RL significantly increased RNA integrity when compared to LN2 storage method. Other studies investigated the effect of various euthanasia methods on protein kinases [19] and brain mRNA levels in mice [20, 21]. In comparison to carbon dioxide (CO2) asphyxiation, euthanasia by decapitation without anaesthesia or with ketamine/xylazine or isoflurane anaesthesia induced significant phosphorylation of mitogen activated protein kinases, highlighting a ‘euthanasia-instigated’ difference in the protein kinases activity. Similarly, ‘focused beam microwave irradiation’ method of euthanasia in mice significantly reduced brain mRNA expression levels, compared to euthanasia by CO2 inhalation. Naïve mice and mice with acute traumatic ‘induced’ brain injury showed varying levels of gene expression in hippocampus when evaluating the effect of euthanasia methods like isoflurane inhalation, chloralhydrate injection and a combination of anaesthetics against an effect of no anaesthesia. Also, the impact of anaesthesia and euthanasia on metabolomics of different mammalian tissues, using a mouse model, has been previously documented [22]. Most notably, elevated levels of multiple nucleotide and purine degradation metabolites were present in skeletal muscle, heart and liver tissues of euthanized animals, indicating deteriorating effect of euthanasia on nucleotides, on comparison with tissues from anaesthetized animals. In pigs, the effect of physical and inhaled euthanasia methods on hormonal measures of stress was previously examined [23], where gradual administration of CO2 or 70% N2/30% CO2 produced similar plasma concentrations of stress indicators to that of physical euthanasia methods.

To our knowledge, no previous research was done at transcriptome level to measure gene expression profiles in porcine tissues that compares the effect of different euthanasia methods and storage conditions. Tissues meant for transcriptome studies are often snap frozen in LN2 after collection and eventually stored at -80 °C until RNA extraction. Due to several advantages over the snap freezing method like improved RNA yield, quality and integrity as well as ease of use in field conditions [18, 24,25,26], a commercial preparation of ammonium salt solution, RNAlater™ (RL), is now popularly used for preserving tissues for gene expression studies. However, there are also recent studies that reported the influence of RL on gene expression. Comparison of microarray profiles from adenocarcinoma tissues in humans indicated that around 2.3% (corresponding to 540 genes) of the total expressed genes displayed more than a two-fold difference in the RL stored samples compared to the snap frozen samples. The upregulated genes were those involved in translational regulation and RNA metabolism and the ones that were downregulated were associated with regulation of enzymatic and energy-dependent processes [14]. Another study involving albino rats reported a significant influence of RL storage on gene expression levels in multiple tissues, measured by qPCR analysis[15]. It was reported that despite being able to maintain RNA integrity, RL is unable to inhibit biological activity and thereby provide a representative gene expression, compared to storage in LN2. A similar finding was also described in a RNAseq study in 30 day old fries (young fish with unabsorbed yolk sac) suggesting a potential role of RL in eliciting a response that can bias the gene expression[16]. Further, evidences for differential expression and post translational modifications were observed in RL’ stored Arabidopsis thaliana tissues, using RNAseq[17]. Hence, it is important to investigate if the storage condition and its interaction with the method of euthanasia would introduce a variation that is not connected to the experimental condition under investigation.

Quantseq is a robust mRNA sequencing method used in transcriptomics. The quantseq 3’ mRNA-Seq (also known as 3’ RNA Tag-Seq method) targets the 3’ poly-A tail of mRNA [27,28,29]. It is cost effective as well as devoid of fragment size bias because it generates a uniform read distribution irrespective of the original RNA length. Compared to whole transcriptome sequencing, quantseq is a suitable alternative especially to quantify gene expression and measure the differential gene expression between two biological conditions, at relatively low sequencing depth (~ 2 million reads per sample) [24,25,26].

The aim of the current study is to investigate whether the choice of euthanasia method (ANOXIA method vs T61), together with the effect of the sample storage condition (RL vs LN2) had an influence on the gene expression profile in pigs. For this, we compared the RNAseq data (quantseq 3’ mRNA) obtained from brain, lung and liver tissues and checked whether the gene expression profiles were affected by the experimental conditions. These tissues were considered the most relevant ones for this study because brain and lung will be exposed immediately to these euthanizing agents eliciting a response whereas liver will play a pivotal role in its metabolism. With this research, we want to highlight the importance of these forgotten variables that can introduce biological or non-biological variations in transcriptomic studies.

Results

Descriptive statistics – RNA extraction

All tissues used in this experiment were sampled and stored within a mean sampling time of 7.75 min (Standard Deviation (SD) = 1.44, Interquartile Range (IQR) 6.40–9.10), although there were differences in sampling time between tissue types (Table 1). High quality RNA was extracted from these tissues and had an average concentration of 674.7 ng/ µl (SD 379.6, IQR 361.4–847.5) RNA. RNA quality parameters comprising A260/230 ratio and A260/280 ratio were 1.99 (SD 0.29, IQR 1.93–2.18) and 2.11 (SD = 0.02, IQR 2.08- 2.13) respectively. The integrity of RNA measured as RIN value had an average measurement of 8.81 (SD = 1.15, IQR 8.28–9.70). A detailed description of the experimental data, based on grouping by euthanasia method, storage condition and tissue types is given in Table 1.

Table 1 RNA measurements and sequencing summary on tissues grouped by euthanasia method and storage condition
Full size table

Summary of sequencing data

Raw sequence data from pituitary, hypothalamus, lung and liver were processed for further analysis. After sequencing, adapter and low quality sequences were trimmed off to attain clean reads. On average, 2.89 ± 0.38 million reads per sample were retained after quality control, and used for further analysis. The reads after alignment to the Sscrofa11.1 genome had an overall alignment rate of around 84.9 ± 1.27%, which included 71.5 ± 3.8% uniquely mapped reads and 13.5 ± 1.27% of the reads mapping to multiple locations on the genome. A summary of sequencing on the basis of grouping by experimental conditions and tissue type is described in Table 1. All the parameters described here had comparable values between the euthanasia methods. However, analysis of variance (ANOVA) revealed a significant effect (p < 0.001) of tissue types and storage conditions on the ‘unique alignment’ percentage (Supplementary Table S1, Additional File 1). Post hoc analysis indicated that the hypothalamus samples (67.5 ± 2.1%) (Fig S1) and samples stored in LN2 (69.9 ± 3.3%) (Fig S2) had a significantly lower ‘unique alignment’ percentage. A detailed description of the sequencing data from all the samples used in this study is given in the supplementary Table S2, Additional File 2.

Differential gene expression analysis

Count data generated following sequence analysis were used for estimating differential gene expression in different tissue types. Pituitary, hypothalamus, lung and liver samples, grouped based on euthanasia method and storage condition, were evaluated for identifying differentially expressed genes. Expression levels of the top 20 highly expressed genes across samples specific to different tissue types is represented in the heatmap (Fig. 1 A –  D). The samples (AnoxRL_Hyp1 and T61RL_Hyp1) representing hypothalamus tissues from two animals, each euthanized by ANOXIA and T61 respectively and stored in RL (Fig. 1 B), showed a different level of gene expression especially for the mitochondrial gene mtCO1. This gene encodes for the cytochrome C oxidase subunit 1, predicted to be involved in mitochondrial electron transport. Barring this deviation, the overall gene expression profile suggests that there is no influence of the euthanasia method on the expression of top ranked genes across different tissues types.

Fig. 1
figure 1

Heat maps of gene expression data from pituitary, hypothalamus, liver and lungs samples (Fig A- D) comparing the two euthanasia methods, ANOXIA and T61. Each row represents a gene and each column represents a sample. The top 20 genes ranked based on magnitude of gene expression are depicted and the colours of the tiles represent the measured expression value of a gene in that sample. The data is normalized using the variance stabilization transformation. The samples are named by combining the euthanasia method and storage condition, followed by the tissue type of the sample For example, AnoxRL_Lv is a liver sample, from an animal euthanized by ANOXIA and stored in RNAlater. The numbers marked on the sample name indicates xth animal per condition

Full size image

Effect of euthanasia methods and storage condition on gene expression profile in pituitary and hypothalamus

Samples from pituitary were stored exclusively in RL and hence the gene expression analysis from this tissue is a direct comparison of the two euthanasia methods, ANOXIA and T61. Differential Gene Expression (DGE) analysis on the gene count data from the pituitary did not identify any statistically significant expression differences between the two euthanasia methods. Out of the 16,777 ‘non zero’ genes from the count matrix that were mapped to the pig genome, none of the protein coding genes were differentially expressed. Although there were no protein coding genes differentially expressed across the two conditions, small nuclear RNAs (snRNAs) namely U1, U2 and U4 that constitute the eukaryotic spliceosomal machinery were found to be significantly higher (Log2 Fold Change (LFC) ≥|2.0|, and adjusted p value ≤ 0.1) in pituitary samples belonging to the ANOXIA group in comparison to T61 (Fig. 2 A).

Fig. 2
figure 2

Enhanced volcano plot identifying differentially expressed genes between the euthanasia methods in different tissue types. In pituitary, no significant deviation in expression of protein coding genes was detected. However, small nuclear RNAs (snRNAs) were found to be significantly higher while using the ANOXIA method (Fig. 2A). No significant differential gene expression was observed in hypothalamus (Fig. 2B), lungs (Fig. 2C) and liver (Fig. 2D). Log2 fold change is plotted on the x axis and negative logarithm (to the base 10) of adjusted p-value is plotted on the y axis

Full size image

In the hypothalamus, no significant effect of the euthanasia methods on the gene expression profile was observed (Fig. 2 B). Out of the 16,599 ‘non zero’ genes that mapped to the pig genome, there was no differential expression detected. However, while comparing the storage methods, it was noticed that the transcripts of the gene cholecystokinin (CCK) appeared to be significantly lower (LFC ≥|2.0| and adjusted p value ≤ 0.1) in the hypothalamus samples stored in RL (Fig. 3A, (See supplementary Table S3, Additional File 3 for full list of genes). It is also noteworthy that no interaction effect of the euthanasia method and storage condition was detected on the gene expression profile in the hypothalamus.

Fig. 3
figure 3

Enhanced volcano plot identifying differentially expressed genes between the storage conditions, from the tissue types Hypothalamus(Fig. 3A), Lungs (Fig. 3B) and Liver (Fig. 3C). Log2 fold change is plotted on the x axis and negative logarithm (to the base 10) of the adjusted p value is plotted on the y axis. Several genes were differentially expressed in liver samples stored in RNAlater

Full size image

Effect of euthanasia methods and storage condition on gene expression profile in lungs and liver

The analysis of gene expression data from lungs mapped a total of 16,805 non zero count genes to the Sus scrofa genome. No genes were found to be differentially expressed at a threshold of LFC ≥|2.0| and adjusted p value ≤ 0.1 (Fig. 2C). Interestingly, the storage condition was found to be a significant factor that influences the detection of differentially expressed genes in lungs. Expression of two genes (quiescin sulfhydryl oxidase 1—QSOX1 and ENSSSCG00000037150 (currently in ensemble archive) were found to be significantly higher (LFC ≥|2.0| and adjusted p value ≤ 0.1) in lung tissues stored in RL, as compared to LN2 (Fig. 3B).

For liver samples, no significant effect of the euthanasia method on gene expression was detected (LFC ≥|2.0| and adjusted p value ≤ 0.1) (Fig. 2D). Interestingly, method of tissue storage had a significant effect on gene expression in the liver tissue, similar to that in lungs. It was observed that out of the 15,778 ‘non zero’ genes from the count matrix mapped to the genome, twelve genes appeared to be upregulated (LFC ≥|2.0| and adjusted p value ≤ 0.1) in the liver tissues stored in RL in comparison with tissues that were snap frozen and stored in LN2 (Fig. 3C, See supplementary Table S4, Additional File 4 for the list of differentially expressed genes). Further, we did not detect any interaction effect of the euthanasia method and storage conditions on the gene expression profile in liver tissue. No genes were found to be differentially expressed when the interaction term was included in the model.

Discussion

In transcriptome studies, different euthanasia methods and tissue storage conditions are routinely used and although this could affect gene expression levels that are unrelated to the experimental design, to our knowledge this has not been investigated so far in pigs. Moreover, the suitability of a method using nitrogen anoxia opposed to T61 injection for euthanizing piglets used in transcriptome research, was not previously tested. Due to the increasing use of pigs in translational research, we have to rule out whether this novel experimental condition can induce a biological variation due to the interaction with the individual or the storage condition.

The present study evaluated the effect of two euthanasia methods (ANOXIA and T61) and two storage conditions (RL and LN2) on the differential gene expression in pituitary, hypothalamus, liver and lungs of male piglets. The statistical model in the differential gene analysis also tested whether there was an interaction effect of these variables on the gene expression profile.

No significant differential expression of any of the known protein coding genes were observed in pituitary and hypothalamus samples, on comparison of the euthanasia methods. Also, for the lung and liver samples, there were only a couple of genes differentially expressed, notably at a very low log2fold change (LFC ≥|0.5|). Considering our small sample size per condition (n = 4), we would ideally set a higher LFC threshold as recommended by Schurch et al.[30], to be absolutely sure about the credibility of the differentially expressed genes. When set at a higher LFC threshold (LFC ≥|2.0|), the expression of these genes in hypothalamus and pituitary was found not to be significantly different between the euthanasia methods studied. When there was no or very little evidence against the null hypothesis, the false discovery rate was very close to one for most of the values and this was evident when the results were plotted with adjusted p value (Fig. 2 A-D). These above observations imply that the nitrogen anoxia method might be a suitable alternative for the T61 injection for euthanizing piglets used in transcriptomics research, as there weren’t any differential expression of protein coding genes. However, there was an interesting observation while analysing the transcriptomic data from the pituitary samples. Small nuclear RNAs (snRNAs) such as U1, U2 and U4, that form a key component of the spliceosomal machinery in the eukaryotes [31, 32], were found to be differentially expressed in pituitary, between the two euthanasia methods. These snRNAs were found to be significantly higher in pituitary samples belonging to the ANOXIA group in comparison to the T61 group.

Differential expression of snRNAs observed across the two conditions tested cannot be overlooked, given the influence of this category of non-coding RNAs on varied cellular processes [33]. These small RNA molecules are found in the cell nucleus and are involved in the processing of pre-messenger RNAs. The sm-class of snRNA genes (classified based on the binding site it possesses for small (Sm) proteins of less than 20 kDa in size) include U1, U2 and U4 snRNAs that play a significant role in catalysing the removal of introns from pre-mRNA [34, 35]. They perform this function by forming a complex with snRNA specific partner proteins that are classified as the small nuclear ribonucleoproteins (snRNP) [35].The gene ontology (GO) terms associated with U1 snRNA is mRNA 5'-splice site recognition and pre-mRNA splice site binding. Similarly, GO terms associated with U2 snRNA is mRNA branch site recognition and pre-mRNA branch point binding. U4 snRNA is also an important member of the spliceosomal process as it is involved in spliceosomal tri-snRNP complex assembly formation. It is also associated with formation of quadruple SL/U4/U5/U6 snRNP as well as U6 snRNA binding [36].

An increase in the level of snRNAs as observed in the nitrogen anoxia method indicates an altered cellular process where snRNA activity was increased. It could be interpreted that these genes were less stable or got degraded during euthanasia with T61 injection, probably as a response to the drug. As a consequence of the euthanasia method used, the splice variant of a gene is either formed or masked, depending on the up/down regulation of the snRNA. Either ways, the consequence of this change must be further investigated as it has the potential to modify the gene expression, like the effect of a splice variant, that we originally intent to measure. In humans, there are evidences for regulation of gene expressions by the U1 snRNAs at the pre-mRNA processing stage [37]. Multiple evidences for down regulation of numerous snRNA genes in response to UV light treatment has also been previously documented [31]. It is worth mentioning that only a small proportion of genes were differentially expressed due to the euthanasia method. However, given the extensive use of pigs in translational research, any effect in response to the genomic stress caused by the method of euthanasia must be ruled out before concluding the research findings.

The impact of storage conditions on RNA measurements and transcriptome profile is equally worth addressing. The integrity of RNA extracted from samples stored in RL had higher mean value compared to the RNA samples that were stored in LN2. This finding is in agreement with the previous finding by Bray et al.[14], where it was reported that RNA integrity is best maintained when tissues are stored in RL. Further analysis of the sequence data indicated that the ‘unique alignment’ percentage was significantly lower for samples that were stored in LN2. We also found evidences for differential gene expression in hypothalamus, liver and lungs, depending on the storage condition. Liver samples reacted the most to the storage condition as several genes appeared to be differentially expressed, upon storage in RL. The inability of RL to inhibit biological activity and thereby provide a stable gene expression [14, 15] might explain the differential gene expression in liver samples stored in RL.

Hypothalamus and lung samples that were stored in RL also showed differential gene expression, when compared to tissues stored in LN2. Several genes appeared to be differentially expressed in both these tissues that were stored in RL. We could not associate any meaningful interpretation while we searched further for gene ontology terms associated with these genes or while a pathway analysis was conducted on this limited number of differentially expressed genes. However, these observations clearly warrants that sample storage condition has the potential to cause a source of technical variation, especially while preserving samples for transcriptome studies. Perhaps, setting higher threshold to identify the differentially expressed genes will ensure that only the highly significant genes are detected in the analysis.

Our study was conducted in a small but accepted sample size for RNAseq studies, and it was specifically done in male piglets. This could be a source of bias and hence if future studies are conducted involving more samples representing both sexes, this can be ruled out. Further, quantseq method of analysis, where only the 3’end of the gene is sequenced, was used to detect the differentially expressed genes in our study. It would be interesting to see if our results are reproducible with other sequencing methodologies as well. However, we could already suggest that ANOXIA method of euthanasia might be a suitable alternative in research animals. We would also like to caution researchers about the impact of these often ‘forgotten’ variables, especially in gene expression studies. Storage medium like RNAlater can influence the differential expression of genes.

Conclusions

The ANOXIA method was found to be comparable to the T61 euthanasia method with regard to gene expression profiles in different tissues. This suggests that the nitrogen anoxia method could serve as a suitable alternative for euthanasia of piglets used in transcriptomic research. However, there could be some regulatory elements that can cause disruptions in the cellular processes and in turn alter the gene expression as was observed in the case of pituitary. Storage of samples in RL proved to better preserve RNA compared to storage in LN2, as indicated by higher RIN values. However, it’s worth noting that RL may be less effective in inhibiting biological activity. Consequently, the choice of storage conditions could potentially introduce confounding variables into our research findings. Researchers must be cautious of these potential sources of bias like euthanasia method and storage conditions, that could influence the research outcomes. Even if we are unable to prevent the effect of these external factors, it will be useful to identify the impact of these factors on the parameter under observation and thereby prevent misinterpretation of our research findings.

Material and methods

Sampling and RNA extraction

The experimental set up, method of euthanasia, sampling and storage conditions and RNA extraction protocol are described in Chakkingal Bhaskaran et al. [18]. In short, 12 one-week-old male piglets were used in the experiment, and 6 animals were euthanized at random either by ANOXIA or by T61. Samples were collected immediately, following euthanasia and confirmation of death. The study was conducted in male piglets since this trial also served as a pilot experiment on transcriptome analysis related to cryptorchidism. Tissue samples from pituitary, hypothalamus, lungs and liver were collected and later processed for transcriptome analysis. For pituitary, samples were only stored in RL, as the gland was too small to be split for the two conditions. Four samples from each tissue type per experimental condition (except lung sample (n = 3) from T61: LN2 group) were processed for the gene expression analysis. A detailed description of the samples (n = 59) used in the transcriptome analysis is given in the supplementary Table S5, Additional File 5.

Library preparation and sequencing

Complementary DNA (cDNA) preparation, library preparation and sequencing was done at the Genomics Core facility of KU Leuven. Double stranded DNA (dsDNA) was constructed from the mRNA by random-primed reverse transcription and second-strand cDNA synthesis. Sequence libraries were prepared with the “3' mRNA-Seq library prep kit for Illumina (FWD)” from QuantSeq-Lexogen (Catn°015.96), following the manufacturer's protocol. Library quality and size range was assessed using a Bioanalyzer (Agilent Technologies, California, USA) with the DNA 1000 kit (Agilent Technologies, California, USA) according to the manufacturer's recommendations. Libraries were diluted to a final concentration of 2 nM and subsequently sequenced on a Illumina HiSeq4000 platform. Single-end reads of 50 bp length were generated with a minimum of around 2 million reads per sample. Quality control of raw reads was performed with FastQC v0.11.7 [38]. Adapters were filtered and trimmed off with ea-utils fastq-mcf v1.05 [39]. Seed based alignment was performed with HISAT2 version 2.1.0 [40] against the pig reference genome Sscrofa11.1, using the parameters (hisat2 -f -x genome -U reads.fa -S output.sam –no-spliced-alignment). Reads mapping to multiple loci in the reference genome were discarded. Resulting SAM alignment files were handled with Samtools v1.5 [41]. Quantification of reads per gene was performed with HT-seq Count v0.10.0 [42]. Analysis of Variance (One way ANOVA) and post hoc analysis (Tukey’s Honest Significant Difference (HSD) test) of the alignment data was used to assess the significance of difference between group means of the fixed factors. All sequences generated in this study is deposited to the Sequence Read Archive (SRA), with accession numbers SRR25184123 to SRR25184181 (BioProject number: PRJNA992446).

Differential expression analysis

Count-based differential expression analysis was done with the Bioconductor package DESeq2 [43] in R. Samples representing different tissue types and conditions were normalized for differences in their sequencing depth. The size factor for each sample was calculated and the count data were divided by the size factor for normalization. The count data were further transformed using the function ‘vst’ for variance stabilization, before visualization and clustering. Reported p-values were adjusted for multiple testing with the Benjamini–Hochberg procedure, which controls false discovery rate (FDR) [44]. The cut-off values for log2FoldChange (LFC), p-values and p-adjusted values for a sample size ranging between 4–12 is |0.5 – 2.0|, 10e-6 and 0.1 respectively, as recommended by Love et al. [43] and Schurch et al. [30]. The statistical model used in the analysis included the main effects of euthanasia method and storage condition, as well as the term indicating the interaction effect (euthanasia:storage). However, the model was re-run without the interaction term when it was found to be statistically non- significant.

Availability of data and materials

All transcriptome data generated in this study are deposited in the NCBI repository and are available in SRA under the bioproject PRJNA992446 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA992446).

Abbreviations

RL:

RNAlater™

LN2 :

Liquid nitrogen

RIN:

RNA integrity number

CO2 :

Carbon dioxide

SD:

Standard deviation

IQR:

Interquartile range

snRNAs:

Small nuclear RNAs

LFC:

Log2 fold change

snRNP:

Small nuclear ribonucleoproteins

cDNA:

Complementary DNA

mtCO1 :

Mitochondrial cytochrome C oxidase subunit 1

CCK :

Cholecystokinin

QSOX1 :

Quiescin sulfhydryl oxidase 1

GO:

Gene ontology

References

  1. Jin L, Tang Q, Hu S, Chen Z, Zhou X, Zeng B, et al. A pig BodyMap transcriptome reveals diverse tissue physiologies and evolutionary dynamics of transcription. Nat Commun. 2021;12:3715. https://doi.org/10.1038/s41467-021-23560-8.

  2. Schroyen M, Tuggle CK. Current transcriptomics in pig immunity research. Mamm Genome. 2015;26:1.

    Article  CAS  PubMed  Google Scholar 

  3. Yang Y, Yan J, Fan X, Chen J, Wang Z, Liu X, et al. The genome variation and developmental transcriptome maps reveal genetic differentiation of skeletal muscle in pigs. PLoS Genet. 2021;17:e1009910.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Wu T, Zhang X, Tian M, Tao Q, Zhang L, Ding Y, et al. Transcriptome analysis reveals candidate genes involved in splay leg syndrome in piglets. J Appl Genet. 2018;59:475–83.

    Article  CAS  PubMed  Google Scholar 

  5. Kobayashi E, Hishikawa S, Teratani T, Lefor AT. The pig as a model for translational research: overview of porcine animal models at Jichi Medical University. Transplant Res. 2012;1:1–9.

    Article  CAS  Google Scholar 

  6. Dawson HD. A Comparative Assessment of the Pig, Mouse and Human Genomes Structural and Functional Analysis of Genes Involved in Immunity and Inlammation. In: McAnulty PA, Dayan A, Hastings KH, Ganderup N-C, editors. The Minipig in Biomedical Research. Boca Raton, FL: CRC Press; 2011. p. 321–41.

    Google Scholar 

  7. Dawson HD, Loveland JE, Pascal G, Gilbert JGR, Uenishi H, Mann KM, et al. Structural and functional annotation of the porcine immunome. BMC Genomics. 2013;14:332. https://doi.org/10.1186/1471-2164-14-332.

  8. Pehböck D, Dietrich H, Klima G, Paal P, Lindner KH, Wenzel V. Anesthesia in swine: Optimizing a laboratory model to optimize translational research. Anaesthesist. 2015;64:65–70.

    Article  PubMed  Google Scholar 

  9. Lowe JWE. Humanising and dehumanising pigs in genomic and transplantation research. Hist Philos Life Sci. 2022;44:1–27.

    Article  Google Scholar 

  10. Swindle MM, Makin A, Herron AJ, Clubb FJ, Frazier KS. Swine as models in biomedical research and toxicology testing. Vet Pathol. 2012;49:344–56.

    Article  CAS  PubMed  Google Scholar 

  11. Close B, Croft B, Lane B, Knoll B. Recommendations for euthanasia of experimental animals: Part 1. Lab Anim. 1996;30:293–316.

    Article  CAS  PubMed  Google Scholar 

  12. Close B, Banister K, Baumans V, Bernoth EM, Bromage N, Bunyan J, et al. Recommendations for euthanasia of experimental animals: Part 2. DGXT of the European Commission. Lab Anim. 1997;31(1):1–32. https://doi.org/10.1258/002367797780600297.

  13. AVMA. AVMA guidelines for the euthanasia of animals: 2020 edition. Members of the Panel on Euthanasia AVMA Staff Consultants. Schaumburg, IL: AVMA; 2020. https://www.avma.org/sites/default/files/2020-02/Guidelines-on-Euthanasia-2020.pdf. Accessed 11 May 2022.

  14. Bray SE, Paulin FEM, Fong SC, Baker L, Carey FA, Levison DA, et al. Gene expression in colorectal neoplasia: modifications induced by tissue ischaemic time and tissue handling protocol. Histopathology. 2010;56:240–50.

    Article  PubMed  Google Scholar 

  15. Ozkan H, Kerman E. Comparative evaluation of RNAlater solution and snap frozen methods for gene expression studies in different tissues. Rev Rom Med Lab. 2020;28:287–97.

    Google Scholar 

  16. Passow CN, Kono TJY, Stahl BA, Jaggard JB, Keene AC, McGaugh SE. Nonrandom RNAseq gene expression associated with RNAlater and flash freezing storage methods. Mol Ecol Resour. 2019;19:456–64.

    Article  CAS  PubMed  Google Scholar 

  17. Kruse CPS, Basu P, Luesse DR, Wyatt SE. Transcriptome and proteome responses in RNAlater preserved tissue of Arabidopsis thaliana. 2017. https://doi.org/10.1371/journal.pone.0175943.

    Article  Google Scholar 

  18. Chakkingal Bhaskaran B, Meyermans R, Gorssen W, Maes GE, Janssens S, Buys N. A comparative study on the effect of euthanasia methods and sample storage conditions on RNA yield and quality in porcine tissues. Animals. 2023;13:698.

    Article  PubMed  PubMed Central  Google Scholar 

  19. Ko MJ, Mulia GE, van Rijn RM. Commonly used anesthesia/euthanasia methods for brain collection differentially impact MAPK activity in male and female C57BL/6 mice. Front Cell Neurosci. 2019;13:96.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  20. Zhang H, Good DJ. Comparison of hypothalamic mRNA levels in mice euthanized by CO2 inhalation and focused-beam microwave irradiation. Lab Animal. 2011;40(10):313–8. https://doi.org/10.1038/laban1011-313.

  21. Staib-Lasarzik I, Kriege O, Timaru-Kast R, Pieter D, Werner C, Engelhard K, et al. Anesthesia for euthanasia influences mRNA expression in healthy mice and after traumatic brain injury. J Neurotrauma. 2014;31:1664–71.

    Article  PubMed  PubMed Central  Google Scholar 

  22. Overmyer KA, Thonusin C, Qi NR, Burant CF, Evans CR. Impact of Anesthesia and Euthanasia on Metabolomics of Mammalian tissues: studies in a C57BL/6J mouse model. PLoS ONE. 2015;10:e0117232.

    Article  PubMed  PubMed Central  Google Scholar 

  23. Meyer RE, Whitley JT, Morrow WEM, Stikeleather LF, et al. Effect of physical and inhaled euthanasia methods on hormonal measures of stress in pigs. J Swine Health Prod. 2013;21(5):261–9. https://www.aasv.org/shap/issues/v21n5/v21n5p261.pdf.

  24. Mutter GL, Zahrieh D, Liu C, Neuberg D, Finkelstein D, Baker HE, et al. Comparison of frozen and RNALater solid tissue storage methods for use in RNA expression microarrays. BMC Genomics. 2004;5:88. https://doi.org/10.1186/1471-2164-5-88.

  25. Weber DG, Casjens S, Rozynek P, Lehnert M, Zilch-Schöneweis S, Bryk O, et al. Assessment of mRNA and microRNA stabilization in peripheral human blood for multicenter studies and biobanks. Biomark Insights. 2010;2010:95–102.

    Google Scholar 

  26. Hatzis C, Sun H, Yao H, Hubbard RE, Meric-Bernstam F, Babiera GV, et al. Effects of Tissue Handling on RNA Integrity and Microarray Measurements From Resected Breast Cancers. JNCI Journal of the National Cancer Institute. 2011;103:1871–83.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  27. Lohman BK, Weber JN, Bolnick DI. Evaluation of TagSeq, a reliable low-cost alternative for RNAseq. Mol Ecol Resour. 2016;16:1315–21.

    Article  CAS  PubMed  Google Scholar 

  28. Ma F, Fuqua BK, Hasin Y, Yukhtman C, Vulpe CD, Lusis AJ, et al. A comparison between whole transcript and 3’ RNA sequencing methods using Kapa and Lexogen library preparation methods. BMC Genomics. 2019;20(1):9. https://doi.org/10.1186/s12864-018-5393-3.

  29. Moll P, Ante M, Seitz A, Reda T. QuantSeq 3′ mRNA sequencing for RNA quantification. Nat Methods. 2014;11:i–iii. https://doi.org/10.1038/nmeth.f.376.

  30. Schurch NJ, Schofield P, Gierliński M, Cole C, Sherstnev A, Singh V, et al. How many biological replicates are needed in an RNA-seq experiment and which differential expression tool should you use? RNA. 2016;22:839–51.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. Jawdekar G, Henry R. Transcriptional regulation of human small nuclear RNA genes Biochimica et Biophysica Acta (BBA). Gene Regulatory Mechanisms. 2008;1779:295–305.

    CAS  Google Scholar 

  32. Guiro J, Murphy S. Regulation of expression of human RNA polymerase II-transcribed snRNA genes. Open Biol. 2017;7(6):170073. https://doi.org/10.1098/rsob.170073.

  33. Matera AG, Terns RM, Terns MP. Non-coding RNAs: lessons from the small nuclear and small nucleolar RNAs. Nature Reviews Molecular Cell Biology. 2007;8:209–20.

    Article  CAS  PubMed  Google Scholar 

  34. Pollard TD, Earnshaw WC, Lippincott-Schwartz J, Johnson GT. Chapter 11 - Eukaryotic RNA Processing. In: Cell Biology (Third Edition). Elsevier; 2017. p. 189–207. https://doi.org/10.1016/B978-0-323-34126-4.00016-5.

  35. Grainger RJ, Beggs JD. Pre-mRNA Splicing. In: Maloy S, Hughes K, editors. Brenner’s Encyclopedia of Genetics. 2nd ed. San Diego: Academic Press; 2013. p. 442–5.

    Chapter  Google Scholar 

  36. Cunningham F, Allen JE, Allen J, Alvarez-Jarreta J, Amode MR, Armean IM, et al. Ensembl 2022. Nucleic Acids Res. 2022;50:D988–95.

    Article  CAS  PubMed  Google Scholar 

  37. O’Reilly D, Dienstbier M, Cowley SA, Vazquez P, Drozdz M, Taylor S, et al. Differentially expressed, variant U1 snRNAs regulate gene expression in human cells. Genome Res. 2013;23:281.

    Article  PubMed  PubMed Central  Google Scholar 

  38. Andrew S. FastQC: A Quality Control tool for High Throughput Sequence Data. 2010. https://www.bioinformatics.babraham.ac.uk/projects/fastqc/. Accessed 25 Mar 2022.

  39. Aronesty E. ea-utils: Command-line tools for processing biological sequencing data. 2011. https://github.com/ExpressionAnalysis/ea-utils. Accessed 25 Mar 2022.

  40. Kim D, Paggi JM, Park C, Bennett C, Salzberg SL. Graph-based genome alignment and genotyping with HISAT2 and HISAT-genotype. Nat Biotechnol. 2019;37(8):907–15. https://doi.org/10.1038/s41587-019-0201-4.

  41. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The sequence alignment/map format and SAMtools. BIOINFORMATICS APPLICATIONS NOTE. 2009;25:2078–9.

    Google Scholar 

  42. Anders S, Pyl PT, Huber W. HTSeq—a Python framework to work with high-throughput sequencing data. Bioinformatics. 2015;31:166–9.

    Article  CAS  PubMed  Google Scholar 

  43. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Benjamini Y, Yekutieli D. The control of the false discovery rate in multiple testing under dependency. Ann Stat. 2001;29:1165–88.

    Article  Google Scholar 

Download references

Acknowledgements

The authors wish to thank KU Leuven TRANSfarm facility for their help in procuring the animals and also during sampling activities. Technical help provided by Jess Bouhuijzen Wenger during RNA extraction and other lab procedures is hereby duly acknowledged. The authors also would like to acknowledge the Genomics Core facility of KU Leuven for their technical help in sequencing and preliminary analysis of the data.

Funding

This research was funded by KU Leuven C2 project, grant number C24/18/036″. Roel Meyermans was also partially funded by an SB PhD fellowship (1S37119N) and Wim Gorssen was funded by a FR PhD fellowship (1104320N) of the Research Foundation Flanders (FWO). The funding bodies did not play a role in the design of the study, and collection, analysis, and interpretation of data or in writing the manuscript.

Author information

Authors and Affiliations

  1. Department of Biosystems, Centre for Animal Breeding and Genetics, KU Leuven, Kasteelpark Arenberg 30, Box 2472, Leuven, 3001, Belgium

    B. Chakkingal Bhaskaran, R. Meyermans, W. Gorssen, S. Janssens & N. Buys

  2. Centre for Human Genetics, Genomics Core, UZ-KU Leuven, Leuven, Belgium

    G. E. Maes

  3. Department of Biosystems, Laboratory of Livestock Physiology, KU Leuven, Kasteelpark Arenberg 30, Box 2472, Leuven, 3001, Belgium

    J. Buyse

Authors
  1. B. Chakkingal Bhaskaran
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. R. Meyermans
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. W. Gorssen
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. G. E. Maes
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. J. Buyse
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. S. Janssens
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. N. Buys
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

Conceptualization: N.B., S.J., G.M., J.B., B.C.B., R.M. and W.G.; Formal analysis: B.C.B.; Funding acquisition: N.B., SJ; methodology: N.B., S.J., G.M., J.B., B.C.B., R.M. and W.G.; project administration: N.B and S.J.; supervision: N.B and S.J.; original draft: B.C.B.; writing-review and editing, N.B.,S.J.,B.C.B.,G.M., J.B., R.M. and W.G. All authors have read and agreed to the published version of the manuscript.

Corresponding authors

Correspondence to B. Chakkingal Bhaskaran or N. Buys.

Ethics declarations

Ethics approval and consent to participate

The experiment was conducted in accordance with the guidelines and approval of the Ethical Committee Animal Experimentation (ECD), KU Leuven [Ref No: P198/2018]. Further, the study is reported in compliance with the ARRIVE guidelines.

Consent for publication

Not applicable.

Competing interests

The authors have declared that no conflict of interest exists.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1:

Fig. S1. Comparison of mean ‘Unique Alignment Percentage’ between different tissue types. This figure illustrates the comparison of the mean ‘Unique Alignment Percentage’ across various tissue types - Pituitary, Hypothalamus, Lungs and Liver. The results of the post hoc analysis confirm significant differences between these tissue types, with the average ‘Unique Alignment Percentage’ in the hypothalamus being notably different from other tissue types. In the figure, Ranges (bars) and means (dots) indicate the spread and the average values of the ‘Unique Alignment Percentage’ respectively. Groups carrying different superscript letters (a, ab, b and c) and colours differ significantly (p < 0.05).

Additional file 2:

Fig. S2. Comparison of mean ‘Unique Alignment Percentage’ between different storage conditions. This figure illustrates the comparison of the mean ‘Unique Alignment Percentage’ across different storage conditions - RNAlater (RL) and LN2 (Liquid Nitrogen). The result of the post hoc analysis confirms significant differences between the two storage conditions, with the average ‘Unique Alignment Percentage’ for samples snap frozen in LN2 being significantly lower. In the figure, Ranges (bars) and means (dots) indicate the spread and the average values of the ‘Unique Alignment Percentage’ respectively. Groups carrying different superscript letters (a and b) and colours differ significantly (p < 0.05).

Additional file 3:

Table S1. Results of analysis of variance (One Way ANOVA) for ‘Unique Alignment’ of reads to reference genome, according to different sources of variation.

Additional file 4:

Table S2. Detailed summary of samples used in the study, with RNA quality parameters as well as complete information of sequence alignment.

Additional file 5:

Table S3. List of differentially expressed genes ((LFC>= |0.5 - 2.0|, padj <= 0.1), detected on comparison of storage methods (RL vs LN2) in hypothalamus tissues.

Additional file 6:

Table S4. List of differentially expressed genes ((LFC>= |0.5|, padj <= 0.1), detected on comparison of storage methods (RL vs LN2) in liver tissues.

Additional file 7:

Table S5. Overview of samples selected for Quantseq Analysis - Samples were grouped by combining the factors euthanasia method and storage Condition.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Chakkingal Bhaskaran, B., Meyermans, R., Gorssen, W. et al. The forgotten variable? Does the euthanasia method and sample storage condition influence an organisms transcriptome – a gene expression analysis on multiple tissues in pigs. BMC Genomics 24, 769 (2023). https://doi.org/10.1186/s12864-023-09794-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-023-09794-4

Keywords

BMC Genomics

ISSN: 1471-2164

Contact us