Transmissible cancer influences immune gene expression in an endangered marsupial, the Tasmanian devil (Sarcophilus harrisii).

Understanding the effects of wildlife diseases on populations requires insight into local environmental conditions, host defence mechanisms, host life-history trade-offs, pathogen population dynamics, and their interactions. The survival of Tasmanian devils (Sarcophilus harrisii) is challenged by a novel, fitness limiting pathogen, Tasmanian devil facial tumour disease (DFTD), a clonally transmissible, contagious cancer. In order to understand the devils' capacity to respond to DFTD, it is crucial to gain information on factors influencing the devils' immune system. By using RT-qPCR, we investigated how DFTD infection in association with intrinsic (sex and age) and environmental (season) factors influences the expression of 10 immune genes in Tasmanian devil blood. Our study showed that the expression of immune genes (both innate and adaptive) differed across seasons, a pattern that was altered when infected with DFTD. The expression of immunogbulins IgE and IgM:IgG showed downregulation in colder months in DFTD infected animals. We also observed strong positive association between the expression of an innate immune gene, CD16, and DFTD infection. Our results demonstrate that sampling across seasons, age groups and environmental conditions are beneficial when deciphering the complex ecoevolutionary interactions of not only conventional host-parasite systems, but also of host and diseases with high mortality rates, such as transmissible cancers.

INTRODUCTION

Wildlife diseases, notably (emerging) infectious diseases, are increasingly recognized as drivers of wildlife population dynamics with the potential to negatively impact the long‐term survival of species and biodiversity (Cunningham et al., 2017; Daszak et al., 2000; Fisher et al., 2012; Robinson et al., 2010). Among these infectious diseases are transmissible cancers, where the pathogen is a clonal malignant cell line, a cell type derived from the host that has crossed the threshold of contagiousness. Although transmissible clonal cell lines represent a rare type of pathogen, their ecological and conservation consequences can be significant (Metzger & Goff, 2016).

An emblematic example of transmissible cancer being a selective force, is the Tasmanian devil, whose population has been reduced by >85% in the last 20 years due to a contagious clonal cell line, named devil facial tumor disease (DFTD) (Cunningham et al., 2021; Hamede, Owen, et al., 2020; Lazenby et al., 2018). DFTD was first detected in 1996 (Hawkins et al., 2006) with a second transmissible cancer, DFT2, recently discovered (Pye, Pemberton, et al., 2016). The emergence of two independently evolved cancer lineages in Tasmanian devils suggests that devils are prone to transmissible cancers and that these diseases may occur more often than previously reported. Although devils have a functional immune system with a complete immune gene repertoire (Kreiss et al., 2008; Morris, Cheng, et al., 2015; Patchett et al., 2015; van der Kraan et al., 2013), DFTD has, until recently, been lethal to most devils (Hamede, Madsen, et al., 2020).

Over the last decade devils have begun to show adaptations to the disease (Epstein et al., 2016), devils surviving longer with visible tumours (Wells et al., 2017) and tumour regressions occurring in some populations (Pye, Hamede, et al., 2016). Importantly, some of the devils with tumour regressions have showed modifications in the activation of tumour suppressor genes (Margres et al., 2020). Although adaptative defences against DFTD have been emerging in Tasmanian devils, we currently lack detailed information on specific host immune responses, which, if known, could assist in conservation efforts, for example, development of a targeted vaccine (Owen & Siddle, 2019), or help in the identification of an early biomarker of DFTD (and DFT2) that would allow the detection of infection prior to the appearance of visible symptoms. Information on how extrinsic and intrinsic factors (e.g., season, age and sex) affect the hosts’ immune activity would assist with understanding specific host reactions to DFTD.

In order to understand how DFTD infection affects Tasmanian devils, it is not only crucial to compare DFTD infected and noninfected individuals, but also to account for other environmental factors and life‐history traits that can influence immune gene expression, and potentially devil immune responses. For example, changes in blood biochemistry and immune cell counts (Peck et al., 2016), decreased T cell repetoir (Cheng et al., 2019) and decrease in body condition (Ruiz‐Aravena et al., 2018) have been reported in DFTD infected animals compared to healthy devils. Furthermore, previous research detected higher leucocyte counts in the blood of healthy female devils (Stannard et al., 2016), and also that female and male devils showed different ability to cope with DFTD (Ruiz‐Aravena et al., 2018). These findings suggest possible underlying sex specific differences in devil immunity, similar to those reported across numerous other species (Klein & Flanagan, 2016). Additionally, Fraik et al. (2019) used a transcriptome wide approach (from ear biopsies of 20 non‐DFTD affected animals, across three geographic locations) to identify environmental and sex driven patterns of gene expression in Tasmanian devils. While they did not detect transcriptome wide variation in gene expression levels across geographic locations, their gene set enrichment analyses identified differentially expressed genes implicated in local adaptation to abiotic environment (i.e., adaptation to coastal environments including salinity tolerance). However their study did not investigate how gene expression profiles change in relation to DFTD infection status across seasons and age groups. Immune gene expression profiling of devil whole blood samples has previously demonstrated that devils are subjected to immunosenescence (Cheng et al., 2017). Immunosenescence has frequently been observed in both domestic (Day, 2010) and wild animals (Peters et al., 2019) and has been attributed to a reduction in the force of natural selection beyond peak reproductive age or a trade‐off with early life reproductive investment (Peters et al., 2019). Seasonal changes in white blood cell numbers (Peck et al., 2015), haemoglobin and total blood protein levels (Stannard et al., 2016) have been reported in devils, but seasonal variations in immune gene expressions have not yet been determined, despite being observed in other vertebrates (Brown et al., 2016; Dopico et al., 2015; Maher & Higgins, 2016). The underlying stimuli of seasonal variations in immune profiles have been proposed to be either photoperiod or resource trade‐offs with reproduction (Saino et al., 2000), as well as parasite infection or resource availability (Nelson, 2004).

Our study, therefore, aimed to investigate how the expression of immune genes in the blood of Tasmanian devils vary as function of season, age, sex and how these variations are modified by DFTD infections. As RNA from stabilized whole blood has been shown to be suitable for comprehensive immune gene expression profiling (van der Sijde et al., 2020), we measured the expression of immune genes from whole blood samples with qRT‐PCR. Although using peripheral blood and qRT‐PCR admittedly have their limitations, we chose these techniques because peripheral blood has been considered as a viable method for noninvasive, sequential sampling for measuring gene expression in endangered animals (Byrne et al., 2019; Huang et al., 2016; Whitney et al., 2003).

The limitations of our approach include: (i) mRNA expression and protein expression not being the same, (ii) blood‐borne cells frequently not being functionally active until they reach the sites of their functional effect, and (iii) mRNA levels not being clear indicators whether they are produced by more or less cells, or by up/downregulation of genes in a single cell. Nevertheless, Chaussabel et al. (2010) described “blood being the pipeline of the immune system” and showed the potential to use blood transcriptomics to assess the immune system. Immune‐related gene expression profiles have been used regularly to identify infection related transcriptional signatures (Ramilo et al., 2007), to monitor changes of the immune system (van der Sijde et al., 2020), and to detect cancer (Aarøe et al., 2010; Sakai et al., 2019; Showe et al., 2009). Human studies have shown positive correlations between mRNA and protein levels (R2: 0.39–0.79) (Buccitelli & Selbach, 2020; Edfors et al., 2016; Schwanhäusser et al., 2011 but see Reimegård et al., 2021), as well as the possibility to trace back variation in gene expression patterns to the relative proportions of specific blood cell subsets in peripheral blood (Whitney et al., 2003). While most of the human studies use transcriptomics and/or microarrays (but see Showe et al., 2009 who used qRT‐PCR to detect pancreatic ductal adenocarcinoma), recent studies have shown the applicability of qRT‐PCR to measure immune system associated gene expression and immune capacity in Tasmanian devil blood (Cheng et al., 2017; Ujvari et al., 2016). Cheng et al. (2019) demonstrated that, similar to humans (Chomczynski et al., 2016), blood RNA levels positively correlate with both white blood cell and lymphocyte counts in Tasmanian devils, as well as immune gene expression profiles being an appropriate proxy to measure changes in the immune system of these endangered animals. The technique qRT‐PCR also allows the focused, sensitive and cost effective monitoring of gene expression changes from large sample sizes (Liew et al., 2006). In addition, previous studies have shown similar proportion of expressed protein coding genes in blood and various tissue types (60%–80%) (Huang et al., 2016; Liew et al., 2006), showing the applicability of this technique in endangered species where access to tissue is a challenge.

As our intent was to investigate whether DFTD infection modifies the expression of immune system associated genes across seasons, age groups and sex of Tasmanian devils, we decided to use a targeted approach and measured the expression profiles of genes associated with the innate and adaptive immune system (VanGuilder et al., 2008) across 87 Tasmanian devils, sampled over a period of 2 years. To gain a snapshot of overall immune status of Tasmanian devils, we consequently used RNA extracted from whole blood samples and the expression levels of 10 immune genes CD4, CD8α (henceforth referred to as CD8), IgG, IgM, IgA, IgE, NKG2D (alternative name KLRK1), MHC‐II, CD11b (alternative name ITGAM, henceforth referred to as CD11) and CD16 (alternative name FcγRIII[a&b]) were measured (detailed information in Table 1).

TABLE 1

Brief description of genes used in the study, including their functions, cell types they are associated with and a selection of techniques used to detect expression in blood

Gene name (abbreviation)	Cell types the gene is expressed on in blood	Technique used to detect expression in blood	Protein function
Cluster of differentiation 4 (CD4)	All T helper cells1, small % of NK cells, macrophages2, neutrophils3	Also see gene related page on gene cards or for all genes in table4, 5, scRNAseq6, combination scRNAseq, flow cytometry and RT‐qPCR7	Communicates with antigen presenting cells, in T cells usually MHC‐ll, activates a range of immune pathways and responses ultimately leading to lymphokine production, adhesion, motility and activation of T helper cells4. Originally identified as a marker for CD4+ T cells1
Cluster of differentiation 8α (CD8)	Cytotoxic T cells6	scRNAseq6	Mediates cell–cell interactions with immune cells, communicates with antigen presenting cells, marker for CD8+ T cells1,6
Cluster of differentiation 11b or integrin subunit Alpha M (CD11)	Dendritic cells, monocytes, granulocytes, macrophages8	Flow cytometry and RT‐qPCR9, NGS‐multiple RNA data sets8	Mediates leucocyte adhesion and migration, phagocytosis, cell‐mediated killing, chemotaxis and cellular activation10,11
Cluster of differentiation 16 or Fc fragment of IgG Receptor IIIa (CD16)	Natural killer cells6, Neutrophil subset12, Monocyte subset6	scRNAseq6 using mAb13, Immuflorescence and PCR, neutrophil isolation and RT‐qPCR12	Marker for natural killer cells, activates the antigen dependent cytotoxicity cascade (ADCC), stimulates phagocytosis, recognises unknown tumours1,14,15
Immunoglobulin G (IgG)	Circulating antibodies, as well as B cells16,17	scRNAseq16,17	Antitoxin, stimulates phagocytosis, activates compliment pathway and ADCC1
Immunoglobulin M (IgM)	Circulating antibodies, as well as B cells16,17	scRNAseq16,17	Responds to infectious organisms, activates complement and ADCC1,18
Immunoglobulin A (IgA) (serum)	Circulating antibodies, as well as B cells16,17	scRNAseq16,17	Neutralizes pathogens (bacteria and virus) and exotoxins, weak activator of complement, activates ADCC via neutrophils1,19,20 enhances phagocytosis, both pro and anti‐inflammatory (depends on binding receptor)21, release of cytokines, immune cell recruitment and induction of necrosis22,23
Immunoglobulin E (IgE)	Circulating antibodies, as well as B cells16,17	scRNAseq16,17	Activation of allergies, parasite resistance16, activates ADCC and antigen dependent cytotoxicity phagatosis1,24
Major histocompatibility complex class ll (MHC‐ll)	Antigen presenting cells25, both innate and adaptive, including B cells, monocytes, macrophages, and dendritic cells 26,27	Combination scRNAseq, flow cytometry and RT‐qPCR7	Presents peptides to CD4+ T cells1,26
Natural killer group 2D (NKG2D)	Natural killer cells, subset of CD8+ T cells28	Using cytotoxic assays, flow cytometry and mAb28,29	Binds to ligands upregulated in response to cellular stress (including malignant or infected cells)30,31 and stimulates cytotoxic pathways32,33

1(Roitt et al., 2001), 2(Lodge et al., 2017), 3(Biswas et al., 2003) 4(Weizmann Institute of Science, 2021), 5(Harvard, 2021), 6(Lei et al., 2021), 7(Katzenelenbogen et al., 2020), 8(Dang et al., 2020), 9(Swirski et al., 2009), 10(Solovjov et al., 2005), 11(Schmid et al., 2018), 12(Di Fulvio & Gomez‐Cambronero, 2005), 13(Ravetch & Perussia, 1989), 14(Saleh et al., 1995), 15(Yeap et al., 2016), 16(Croote et al., 2018), 17(Schraven et al., 2021), 18(Scott, Wolchok, & Old, 2012), 19(Staff et al., 2012), 20(Macpherson et al., 2008), 21(Davis et al., 2020), 22(Breedveld & van Egmond, 2019), 23(Heineke & van Egmond, 2017), 24 (Karagiannis et al., 2007), 25(Goldinger et al., 2015), 26 (Rock et al., 2016), 27(Bovin et al., 2004), 28(Rosen et al., 2004), 29(Salih et al., 2003), 30(Raulet et al., 2013), 31(López‐Larrea, López‐Soto, & González, 2010), 32(López‐Soto et al., 2015), 33(Moretta et al., 2001).

These genes are known to be highly expressed in various immune cells in blood (for detection techniques, see Table 1), including adaptive immune cells: B cells (IgG, IgM, IgA, IgE: Croote et al., 2018; Roitt et al., 2001; Schraven et al., 2021), T cells (CD4: Science, 2021, CD8: Lei et al., 2021) and innate immune cells: natural killer cells (NKG2D: Rosen et al., 2004; Salih et al., 2003, CD16: Houchins et al., 1991; Lei et al., 2021; Moretta et al., 2001; Ravetch & Perussia, 1989), macrophages (CD11: Dang et al., 2020; Swirski et al., 2009) and neutrophils (CD16: Di Fulvio & Gomez‐Cambronero, 2005; Ravetch & Perussia, 1989). The chosen genes are also known to be involved in recognition and response to abnormal cells through a variety of immune pathways, and some of the genes (NKG2D reviewed in Frazao et al., 2019, CD16: Cheng et al., 2021; Frazao et al., 2019), cells (CD4 & CD8 T cells, natural killer [NK] cells: Davis et al., 2017; Rittig et al., 2011; Tay et al., 2021) and antibodies (IgG, IgM, IgA, IgE: Almagro et al., 2018; Karagiannis et al., 2017) are also utilized or stimulated in cancer immunotherapy and vaccines. The relative binding (and subsequent pathway stimulation) of two pairs of the selected immune genes (CD4 and CD8, IgM and IgG), illustrated as a ratio, have also been used as an indicator of immune health and a predictor of patient prognosis in viral and autoimmune diseases as well as cancer (Li et al., 2019; Serrano‐Villar et al., 2014) and as predictors in cancer treatment outcomes (Ravindranath et al., 1998; Staff et al., 2012). While immune capacity is commonly determined by looking at immune cell counts (Beechler et al., 2012; Kreiss et al., 2008), measuring the expression levels of immune genes has also been used (as an indicative rather than absolute measure) in wildlife studies (Maher et al., 2019; Ujvari et al., 2016).

By analysing the expression profiles of the selected genes, we endeavoured to answer the following questions: Will DFTD infection influence devil immune gene expression profiles? Are there sex specific differences in the devils’ immune system activity? If so, do they explain the increased cancer tolerance and higher DFTD regression rates observed in females (Pye, Hamede, et al., 2016; Ruiz‐Aravena et al., 2018)? Do other variables (age, season, and sex) affect immune gene expression in Tasmanian devils? By elucidating the association between the aforementioned external and internal factors, we aimed to increase our knowledge on hosts’ susceptibility and response to DFTD infection, thus gaining insight into the immune system of the species.

MATERIALS AND METHODS

Fieldwork and sample collection

Fieldwork was undertaken at West Pencil Pine (WPP hereafter), a 252 km area west of Cradle Mountain in north‐western Tasmania (41º 31ʹ S, 145º 46ʹ E). The study site is located in the north‐west semi alpine region of Tasmania, on a large area used for timber production. Devils were caught, aged and scored for DFTD using protocols in Hamede et al. (2015). The first time devils were captured they were individually marked with a microchip and aged using a combination of molar eruption, molar wear and the distance from the dentine‐enamal junction to the gum. Juvenile devils have no visible distance between the dentine enamel junction, the distance becomes visible in the second year and increases as the devils age. As the sites are regularly monitored and the animals are being recaptured, their age can be easily determined based on previous captures (Hamede et al., 2013). Sampling occurred seasonally at three‐month intervals (February, May, August and November) between November 2016 and 2018. These sampling times represent different life‐history stages of Tasmanian devils: February (summer): just before the start of breeding season (as devils are not interested in traps during the breeding season); May (late autumn): period of pregnancy and/or females carrying very small pouch young; August (winter): period of extensive lactation with large pouch young and/or joeys in the den; November (spring): young devils becoming independent and start dispersing. A total of 87 different individuals (39 male, 48 female) were captured for the study, with 39 of these individuals caught either twice (21), three (14) or more (4) times. This low number of repeated samples would lead to over‐parameterisation when accounting for individual in mixed effect models. Thus, to avoid pseudoreplication, only one sample from each individual was selected for the analyses. Animals that only had a single capture were automatically included, and from animals with multiple captures one sample was selected randomly while aiming to achieve as uniformed and balanced distribution of data across the different age, sex and DFTD categories as possible (See Table S1 for additional details).

For a detailed breakdown of samples collected and used in the analyses, see Table S1. Blood was collected using a needle prick to a peripheral ear vein, ~0.5 ml of blood was dripped directly into RNAprotect blood tubes (Qiagen) and immediately stored in liquid nitrogen. Samples were transported to Deakin University in dry shippers via air and then were stored ≤–80°C until analysis.

Obtaining tissue where immune cells perform their immunological tasks (such as lymphatic tissues or the peripheral sites of infection) would have been highly informative. However, as our study used wild animals, we could not conduct invasive sample collection due to animal ethics considerations, and the necessity to resample the animals throughout their life. As detailed above, blood‐based gene expression profiles have previously been used for comprehensive immune gene expression profiling (see for example, Chaussabel et al., 2010; van der Sijde et al., 2020), as well as identifying potential biomarkers (Ramilo et al., 2007). Since invasive sample collection is not possible in our endangered species, we opted to use whole blood samples.

Due to the remote location of the field sites and extreme field conditions (i.e., snow, rain etc.), blood samples were collected directly into the RNAprotect blood tubes that not only preserve the integrity of RNA but also lyse the cells, therefore cell counts could not be determined. Cheng et al. (2017) showed that blood mRNA level positively correlates with white blood cell and lymphocyte counts in our study species, and further studies have also demonstrated a positive correlation between cell counts and gene expression of the genes (CD4, CD8, IgG, CD16, MHC‐II) used in our study (Fujii et al., 2013; Palmer et al., 2006). It is important to note that in our study, without cell number counts, a direct correlation could not be drawn between immune cell numbers and gene expression as immune cells may increase expression per cell rather than the cell count.

Sample collection was carried out according to guidelines and regulations approved by University of Tasmania's Animal Ethics Committee (approval permit A0013326; A0016789) and Deakin University's Animal Ethics Committee (approval number AEX03‐2017). Sampling permit was issued by the Department of Primary Industries, Parks, Water and the Environment.

Genes

Immune responses can be classified as innate/adaptive or Type 1 T helper (Th1)/Type 2 T helper (Th2) type of reactions (Clark & Kupper, 2005). The classical characterisation of innate versus adaptive immune responses is based on distinguishing between fast first line defence against infections (mediated by eosinophils, monocytes, macrophages, natural killer cells, toll‐like receptors) and slow immune responses mediated by the T and B lymphocytes (Clark & Kupper, 2005). In contrast, the Th1/Th2 classification is based on immune responses mediated by cytokines generated by T helper cells, resulting in proinflammatory and anti‐inflammatory reactions (Romagnani, 1997). In our study, we used the classification of innate versus adaptive responses, and classified our genes accordingly, for the following reasons. (i) A previous study focused on Th1/Th2 immune responses in captive Tasmanian devils, and measured the expression of genes in the Th1/Th2 pathway (Cheng et al., 2017). Under captive conditions, the authors were able to collect up to 5 ml of blood, a sample volume that is impossible to obtain under our field conditions (our permit allowed the collection of max. of 0.5 ml of blood). While an attempt was made to amplify the same genes as by Cheng et al. (2017), the expression of these genes was below the detection limit in wild devils. (ii) As described in the introduction, we aimed to target genes involved in general immune responses and cancer surveillance, and genes were chosen as proxy markers of specific immune cells, to use them as indicators of cellular responses, since we were not able to obtain cell count data. Therefore, genes were classified as adaptive or innate based on the immune cells where they are most commonly expressed. For example, T cells (CD4, CD8) (Lei et al., 2021) and B cells (IgG, IgM, IgA, IgE) (Croote et al., 2018) are adaptive immune cells, while NK cells (CD16, NKG2D) (Lei et al., 2021) and macrophages (CD11) are generally classified as innate immune cells (Roitt et al., 2001) (more details in Table 1). This classification does not mean the selected genes are only involved in either adaptive or innate immune responses, most of the chosen genes stimulate multiple immune pathways and can be involved in recruiting additional immune cells. In addition, the ratios of immune gene expression CD4:CD8 and IgM:IgG have been used as relative indicators of binding capacity, the activation of different immune pathways and general immune health (Ravindranath et al., 1998; Tancini et al., 1990).

Sample processing

RNA was extracted with the RNeasy Protect Animal Blood Kit (Qiagen) following the manufacturer's protocol, with some adjustments to improve RNA yield and quality. In steps 1, 3 and 6, 10,000 g centrifugation was used rather than 5,000 g. Step 17 was repeated twice, and 22 μl of buffer REB was directly pipetted onto the spin column, in the repeat, the elute was pipetted back onto the spin column rather than adding additional buffer REB to increase final concentration of RNA in the elute. DNase treatment is included in the extraction kit. Quality and concentration of RNA were determined using the Agilent Bioanalyser 2100 (Agilent Technologies) following the manufacturer's protocols. RNA concentration ranged between 10.5–1,347 ng/µl, with an average of 360.7 ng/µl, and only samples with a RIN of ≥7 were included in downstream analyses. RNA was reverse transcribed to cDNA using Bio‐Rad iScript Advanced cDNA Synthesis Kit for RT‐qPCR that relies on recognising the PolyA tail of mRNA (BioRad). Reverse transcription was conducted for 30 min at 42°C, inactivation was completed for 5 min at 85°C. Primers to amplify the specific genes were designed to span exon boundaries in conserved regions using the Primer3Plus website (Untergasser et al., 2007) (http://www.bioinformatics.nl/cgi‐bin/primer3plus/primer3plus.cgi) from Tasmanian devil genes available through Ensembl and Genbank. (Table 2 and Table S2). Primer quality was checked using the OligoAnalyser3.1 website (Integrated DNA Technologies, 2018) (http://sg.idtdna.com/calc/analyzer ). Reference gene selection was based on BestKeeper2 (Pfaffl et al., 2004) which tests candidate genes for stability within test samples. Altogether 11 samples were chosen across age, sex and DFTD status when conducting the reference gene selection. Similar to previous studies (Cheng et al., 2017; Ujvari et al., 2016) RPS29 had the highest correlation coefficient (r): 0.981 and hence was used in subsequent analyses (RPL4 r: 0.951, OAZ r: 0.694).

TABLE 2

Primers used in the RT‐qPCR reactions

Primer name	Primer sequence	Exon number	Melting temp (oC)	RT‐qPCR efficiency (%)	RT‐qPCR R 2	Amplified fragment length (bp)
qRPS29_F a	ATGGGTCATCAGCAGCTCTAC	1	59.3	105.1	0.996	107
qRPS29_R a	AGGCCGTATTTGCGGATTAG	2	61.3
qOAZ_F a	TGAAATTCCAGGTGGTGCTC	3	61.0	100.3	0.997	90
qOAZ_R a	GACGTGATCAACATGAAGCT	4	59.3
qRPL4_F	AGGCTTGTGTACGTCCCTTG	1	57.2	104.1	0.994	235
qRPL4_R	ACACGAGGAATTCGAGCAAC	2	55.4
qIgG_F a	CAGGTGATCAGCACTCTCTCTG	3	60.3	98.1	0.995	162
qIgG_R a	GGATGTGGGGCAAGACATA	4	59.6
qIgM_F a	TTTGATATCTGGGGCAAAGG	1	59.9	100.0	0.997	119
qIgM_R a	ACAGCAAAGGAGGCATCTTC	2	59.4
qIgA_F	ATCTTCCTGCAAGCCAGTG	1	59.0	99.9	0.992	90
qIgA_R	GTAGTTTCGCAAGGACAATCG	2	59.8
qIgE_F	GAAGACAGTGCCCAAAAGTG	2	58.3	107.7	0.993	123
qIgE_R	CGCTGACATAGAGGTCAAAGG	3	59.9
qCD8_F	CTGGGAAATGCAAGTCCATC	3	60.5	98.9	0.996	107
qCD8_R	GAAGCAAGACAACACAGACACC	4	59.8
qCD4_F	AGAGAACCGAAAGCAGGAAG	4	58.7	103.4	0.991	142
qCD4_R	GACCATTCCATTCCACCTTG	5	60.2
qMHC‐ll_F	GCCCGAGGTGACTGTGTATC	2	60.5	99.0	0.999	61
qMHC‐ll_R	AGACAAGCAGGTTGTGGTGTC	3	60.2
qCD11b_F	ACTGCACGCACTTTCCAAG	7	63.9	97.0	0.997	134
qCD11b_R	GAATTGCTCAAAACCCCTCAG	8	63.1
qCD16_F	CATCACAGCCGACAATATCAC	1	59.0	95.2	0.998	123
qCD16_R	CAGTTTCAAGGTTTGGAGCAG	2	59.9
qNKG2D_F	GACGTGGGAAGATGGTTCAC	7	60.4	94.6	0.994	85
qNKG2D_R	TGGAGCCATAGATTGCACAG	8	59.8

Primers that were developed in a previous study (Ujvari et al., 2016).

Alterations in mRNA levels between samples and across plates were quantitated with standard curves determined using a RT‐qPCR machine (BioRad CFX Connect Real‐Time System). We followed the experimental design of previous studies that generated a composite or “golden sample” to establish standard curves (Foley et al., 2020; Morris, Mathew, et al., 2015; Taylor et al., 2019; Ujvari et al., 2016). The “golden sample” (Foley et al., 2020) was generated by pooling 2 µl of cDNA from each sample, and diluted to a 50 ng/µl stock concentration. Standard curves were generated using serial 1:5 dilutions of the composite cDNA sample, run in triplicate. The 1:5 dilution covered the range from 50 to 0.04 ng concentration, a 625‐fold dilution span. No Cq values of individual values fell outside of the standard curves. Following the establishment of the standard curves, a lower concentration of the composite sample (5 ng/µl) was included in the follow‐up analyses. All standard curves had an R2 > 0.99, contained at least five dilutions from the dilution series with a linear dynamic range of at least three orders of magnitude and had RT‐qPCR efficiencies between 92.4% and 108%, as per standard procedure (Bustin et al., 2009) (Table 2). All RT‐qPCR reactions were performed in 15 µl total volume, containing 7.5 µl of Biorad SsoAdvanced Universal SYBR Green Super Mix (Biorad), 0.5 µM forward and reverse primers and 1 µl of cDNA. RT‐qPCR reactions were run at 95°C for 30 s denaturing and 45 cycles of 95⁰C for 10 s and 60°C for 30 s (annealing temperature). Fluorescence signal was acquired at the annealing temperature. To evaluate the specific amplification, a final melting curve analysis (from 65°C up to 99°C) was added under continuous fluorescence measurements.

The “golden sample” (5 ng/µl) was included on each plate to measure the expression of both the reference and the target genes, thus allowing the standardization of measurements across and between plates (as described by Cheng et al., 2017; Foley et al., 2020; Morris, Mathew, et al., 2015; Taylor et al., 2019; Ujvari et al., 2016). Relative gene expression ratio (herein gene expression) was calculated from the RT‐qPCR efficiencies and crossing point deviation of the sample of interest against the “golden sample” in both the target and reference gene analyses (Pfaffl, 2001). The Pfaffl technique is similar to the delta‐delta Ct method, but it also accounts for any differences in primer efficiencies, and thus increases reproducibility (Pfaffl, 2001). Throughout the design of the RT‐qPCR experiments we adhered to the standards described in Bustin et al. (2009), see Supporting Information (Figure S1 and RT‐qPCR reports). Gel electrophoresis and melt curve analyses validated the lack of primer dimers and genomic DNA contamination in our reactions. Cloning and sequencing of the cDNA fragments confirmed that our primers amplified the targeted regions (see Supporting Information for cloning protocols and GenBank accession MZ516519‐26).

Statistical analysis

All statistical analyses were conducted in R v3.5.2 (R Core Team, 2018). To visualize the Pearson correlation between the expressions of the studied genes we used Corrplot v0.84 (Wei & Simko, 2017). For each gene we used linear models to evaluate the effects of the following explanatory variables on gene expression: sex, season (summer, autumn, winter, spring) the sample was collected, DFTD infection status (presence/absence) and age. Body mass, structural size and condition were not included in the analyses since these variables correlate strongly with age, sex and presence/absence of DFTD in our sample set (Figure S2). Age was calculated as the difference in months between the first of April the year each animal was born and the capture date. While recent studies on Tasmanian devils have shown an expansion of breeding season and shift in breeding age, the data set used demonstrated <10% of females displaying either of these changes, so 1st of April was used as the birth date (Guiler, 1970; Hamede et al., 2012, 2015; McCallum et al., 2009). Previous immunosenescence studies have highlighted the existence of both linear and nonlinear (threshold, quadratic) effects between efficiency of the immune system and age (Cichon et al., 2003; Lund et al., 2002). Therefore, to determine the relationship most appropriate for our data, we ran the model selection step (described below) three times for all genes using (1) a linear relationship with age, (2) a linear relationship with log‐transformed age (for threshold or non linear effects) and (3) a quadratic relationship with age. In all 10 genes the linear relationship with log‐transformed age gave a similar (i.e., within 2 AIC units) or better result (IgA) than the two alternative models (Table S3). Therefore, we used log‐transformed age in all analyses outlined below. All continuous explanatory life‐history variables were scaled to make parameter estimates comparable.

Calculating one full factorial model with all the explanatory variables listed above, as well as their two‐way interaction would result in overparametrisation. Therefore, we excluded all interactions that did not include DFTD and with the remaining main effects and interactions, determined which explanatory variables and interactions were associated with variation in each gene using the function “dredge” within the package MuMIn (Barton, 2020; Loukola et al., 2020; Nagarajan‐Radha & Devaraj, 2021) which evaluated models with all possible combinations of explanatory variables and logical two‐way interactions while respecting the marginality constraints (i.e., no interactions are tested without also including their main effects in the model). Models were ranked and given a weight according to akaike information criterion (AICc) adjusted for sample size. Next, models with substantial empirical support (i.e., models within 2 AICc units) were retained for further analyses and model weights adjusted. For each explanatory variable (both main effects and interactions) in these top models we calculated the importance value by summing the top‐models’ weights in which these terms occurred. These summated importance values can take values between 0 and 1 and only terms with a cumulative importance value greater than 0.3 were selected to form the final, “average” model for each gene.

To determine similarities in the relationships with explanatory variables across all genes we clustered the genes by the importance values of all model terms and next visualized these in a heatmap created using the gplots v 3.0.1.1 (Warnes et al., 2019) and plotrix v 3.7–8 (Lemon, 2006) packages. Using similar heatmaps, we also showed the effects of each variable on gene expression by visualizing the slopes. To display the effect of a specific explanatory variable, accounting for all other variables in the model, marginal‐effect means were calculated in each final gene model using the ggeffects package (Lüdecke, 2018), to test significant differences between categories, post hoc tests were calculated using the emmeans v 1.5.2–1 (Lenth, 2020) and multcomp 1.4–14 (Hothorn et al., 2008) packages with a sidak adjusted p‐value. Marginal means data with post hoc test ouputs in the form of letters were plotted using the ggplots2 v 3.3.2 (Wickham, 2016) and patchwork packages (Pedersen, 2020).

RESULTS

Correlation between the expression of immune genes

The degree of correlation between the expression of different genes (Figure 1) may indicate genes activated or deactivated in the same pathway or responding to similar intrinsic/extrinsic factors. Strong negative correlations were absent in this comparison (i.e., all r ≥ –.3). Particularly high correlations (i.e., r > .7) were all positive and observed between immune genes IgG and IgM, IgA and CD8, NKG2D and CD8, and CD4 and MHC‐II. Otherwise noticeable results, some genes and gene ratios displayed only very low to medium correlations with any of the other genes/gene ratios. These included IgE, CD16 and the two gene ratios IgM:IgG and CD4:CD8.

FIGURE 1

Correlation matrix displaying pearson correlation coeffients calculated between immune gene expression profiles in Tasmanian devils. Matrix is ordered by similarity. Pearson's correlation coefficients are printed in matrix, strength of correlation is represented visually by circle size and colour. Larger circles, higher correlation coefficients (stronger linear relationship); smaller circles, lower correlation coefficients (weaker linear relationship); red, positive correlation; blue, negative correlation

Predictors influencing gene expression profiles in Tasmanian devils

To objectively identify the variables influencing the expression of genes (and thus eliminate researcher bias), model selection analyses were conducted. For each gene this first step provided us with an importance value for each of the seven explanatory variables used (i.e. four main effects and three interactions with DFTD). This part of the analysis showed that DFTD infection was of importance to at least some degree (i.e., >0) in all 12 genes and gene ratios, followed by season (10) and age (9) (Figure S3). Next, for each of the 12 genes a final linear model was created using all exploratory variables for that gene that had an importance value ≥0.3 (Table S4). The parameter estimates (slopes) for these 12 final linear models as well as their significance levels showed that season had a significant effect on the expression of as many as 10 genes with a clear tendency for gene expression to be significantly lower in both autumn and spring compared to summer (Figure 2). Age was the next most determining variable in gene expression with two genes being significantly downregulated with age (IgM and IgM: IgG) and two being upregulated with age (IgA, CD16, and NKG2D). DFTD had a significant effect (p ≤ .05) on gene expression in only four out of the 12 models. In one of these, CD16, the effect was positive, while it was negative for NKG2D, IgG, and IgE. However, DFTD in interaction with the three other investigated factors had a number of additional significant effects. For the interaction of DFTD with sex significant effects were found on IgG and IgE, while for the interaction between DFTD and season significant effects were found on IgM:IgG and IgE (Figure 2).

FIGURE 2

Heatmap as a visual representation of the model estimate (slope) for explanatory variables influencing gene expression profiles in Tasmanian devils, across all genes. For categorical explanatory variables the intercepts are: Season, summer; DFTD, healthy; Sex, female. The specific colour shows the direction of the effect blue, negative slope; red, positive slope. The intensity of colour illustrates the strength of the effect. p‐values calculated from models are indicated as: *** <.001, ** <.01, * <.05, . <.1

The parameter estimates in the final models were used to hierarchically cluster the genes (see dendrograms at top of Figure 2). Interestingly this clustering, depicting the similarity in the collection of explanatory variables affecting gene expression, deviated considerably from what was found in the correlation analyses. Of the four gene pairs showing highly correlated gene expression (Figure 1) only IgA and CD8, and CD4 and MHC‐II could be found side by side in the dendrogram (Figure 2b). The other two pairs (IgG and IgM, NKG2D and CD8) may thus well be highly correlated but still differentially affected by the explanatory variables investigated here.

To visualise the significant effects of the explanatory variables on gene expression, we generated marginal mean plots (i.e., plots where the y‐axis depicts the residual variation in the expression of a gene after correcting for all factors in the final linear model but the explanatory variable of interest) (Figures 3, 4, 5). The marginal mean plots showed similar trends in the effect of season on gene expressions in 10 models, with gene expression being upregulated in summer and often in winter and downregulated in autumn and/or spring (Figure 3). The marginal mean plots also showed trends for increasing gene expression of CD16, NKG2D and IgA and decreasing expression of IgM and IgM:IgG ratios in older devils (Figure 4). DFTD infection had an opposing effect on the expression of NKG2D and CD16, the first showing a decrease and the latter showing an increase in association with DFTD infection (Figure 5). The expression of IgG was influenced by sex in combination with DFTD, and it was downregulated in infected females compared to males. Similarly, the expression of IgE was influenced by the interaction between sex and DFTD, with healthy females displaying higher gene expression compared with males or diseased females. The interaction between DFTD infection and season had an impact on the expression of IgE as well as on IgM:IgG ratio, IgE showed an overall downregulation with DFTD infection, the most significant difference being in winter, while the IgM:IgG ratio showed differing seasonal patterns between healthy and DFTD infected animals.

FIGURE 3

(a) Marginal means effects of DFTD infection, red dots are raw data. (b) Marginal means effects of interactions with DFTD infection and sex for each gene, red, healthy animals; blue, visibly DFTD affected animals. (c) Marginal means effects of interactions with DFTD infection and season for each gene. Black, healthy animals; orange, visibly DFTD affected animals. 95% confidence intervals are displayed. Only statistically significant results are graphed

FIGURE 4

Marginal means effects of season, 95% confidence intervals are displayed. Only statistically significant results are graphed. Marginal means: The effect of a specific explanatory variable, accounting for all other variables in the model, Raw data: light blue dots

FIGURE 5

Marginal means effects of age, black dots are raw data. 95% confidence intervals are shown. Only statistically significant results are graphed

DISCUSSION

Despite substantial research and conservation efforts for more than a decade, Tasmanian devils remain on the list of endangered species, due to the dramatic consequences of the DFTD epizootic (but see Cunningham et al., 2021). Although some animals have shown recovery, most devils still succumb to the disease (Wells et al., 2017). Thus, there is still an urgent need for research, to increase our understanding of the host's immune system and the various intrinsic and extrinsic factors influencing devil immune capacity. Here, therefore, we analysed the expression profiles of immune genes associated with innate and adaptive immune responses to gain detailed understanding of predictors that may influence the Tasmanian devils’ ability to respond to challenges, including DFTD. The major findings of this work were that, DFTD infection had a differential effect on the expression of the studied genes, with some genes showing sex or season dependant changes in expression. Season affected the expression of most of the immune genes, while age influenced the expression about half of the genes analysed.

Complex associations between gene expression levels

Strong correlation was observed between the expression levels of IgA and CD8 genes. As our approach targeted IgA expressed in serum (not mucosal IgA), from both serum IgA antibodies and B‐cells expressing IgA, the observed correlation was expected. Associations between B‐cells expressing IgA and CD8+ T cells have been reported in previous studies, particularly following vaccination or inflammation due to pathogenic challenges (Brown et al., 2001; Polak et al., 2014; Zhu et al., 2017). In contrast, the observed negative correlations between the expression of IgG and CD16, genes that interact to stimulate the antibody‐dependent cellular cytotoxicity (ADCC) pathway (Roitt et al., 2001), was unexpected. The negative correlation is most probably due to age differently affecting the expression of these two genes, as in our observations expression of IgG decreased with age and with DFTD infection, while CD16 increased with both age and DFTD infection (see further below).

Seasonal changes in gene expression

Similar to other studies investigating immune gene expression profiles across seasons in humans (via transcriptomics: Dopico et al., 2015) and voles (Jackson et al., 2011), we observed significant fluctuations in gene expression levels in Tasmanian devils. In general, we observed higher gene expression levels in summer and sometimes winter, and lower in spring/autumn. Three month seasonal fluctuations in the expression of cytokines have also been seen in another marsupial, the koala, but with the opposite peaks and troughs (Maher & Higgins, 2016). The difference between the two studies could be driven by targeting different genes, species occupying different geographic locations and niches (Tasmania vs. Sydney, ground dwelling vs arboreal) or simply differences between the species. Similar to our study, higher expression of immune genes (i.e., transcription factors and cytokines) have been observed during winter in voles (Jackson et al., 2011) a pattern that was assumed to help compensate for the immunosuppressive effect of harsh weather and limited food supplies (Nelson, 2004). While winter also presents the harshest environement for Tasmanian devils in Tasmania, due to DFTD driven drastic reduction in devil numbers, it is unlikely that food is a limiting factor during winter. The observed high immune gene expression during summer may be linked to trade‐offs with reproduction or breeding season (Martin et al., 2008). During late summer, male devils show peaks in androgen (Hesterman & Jones, 2009) and cortisol levels prior to the breeding season (Keeley et al., 2012), thus higher investment in immune function may compensate for the increase in hormone levels (Roved et al., 2017). In female devils progesterone levels rise in March (breeding season) (Hesterman et al., 2008), but no information is available on immune suppressive stress hormones across seasons. Since the emergence of DFTD, the reproductive age of devils has shifted to start between the age of 12 and 24 months (Lachish et al., 2009; Russell et al., 2019), but during this age span it is impossible to distinguish mature from immature animals (unless females have pouch young) even based on size or weight (Russell et al., 2019). As half of devils analysed in this study were within this age span, our data set may contain some sexually immature animals. Follow up studies investigating hormone levels and immune gene expression profiles across sex and age groups would therefore be necessary to confirm whether the observed gene expression patterns during summer were the result of trade‐offs between reproductive effort and immune function maintenance.

In addition to climate and reproductive status, a plethorea of other environmental and physiological factors could underly the seasonal changes in immune gene expression and in immune function, including photoperiod, melanine levels, life history trade‐offs, parasite or pathogens (various topics reviewed in Altizer et al., 2006; Dowell, 2001; Martin et al., 2008; Nelson & Demas, 1996). Furthermore, the observed gene expression patterns could also be the result of other unknown factors such as leucopoiesis triggered by different pathogens, and/or acute stress (e.g., infectious status, injury, stress‐induced activity) (Santini et al., 1995). A notable caveat of our study is that as the expression of some of the genes were correlated, we can therefore not conclude whether these observed seasonal patterns were drivern by general or individual responses to changes across seasons. To further investigate whether the seasonal variation in devil gene expression was initiated by DFTD and/or other pathogens, either single pathogens stimulating multiple parts of the immune system or several pathogens evoking a range of immune responses (Becker et al., 2020), additional data on internal and external pathogen load and prevalence would be required. As it is currently not known if the observed seasonal changes in immune gene expression in devils affect their ability to mount immune responses to local parasites or to DFTD infection, future studies on immune function and DFTD should aim to account for any underlying seasonal changes.

Age associated changes in immune gene expressions

Attenuation of the immune system has been associated with increased parasite susceptibility and subsequent mortality in wild animals (Froy et al., 2019), and reduced IgM:IgG ratios have been linked to poor oncogenic outcomes in humans and mice (Lokshin et al., 2006; Ravindranath et al., 1998). A decrease in the expression of IgM and the IgG:IgM ratio (most probably due to the decrease in IgM) with increasing age in the present study, may indicate some genes are subject to immunosenescence and hence mirrors previous findings (Cheng et al., 2017; Ujvari et al., 2016). However, our results suggest a more complex picture as the expression of IgA, NKG2D and CD16 increased with devil age. This is also supported by Cheng et al. (2017), who showed both increases and decreases of immune gene expression measured with RT‐qPCR from blood through devil puberty. IgA can stimulate the ADCC pathway (via neutrophils) by binding to cell receptors with higher affinity than IgG (Macpherson et al., 2008), therefore the increased expression may indicate enhanced production of IgA as compensation for decreases in the expression of other genes. Futhermore, previous wildlife studies proposed that the stable expression of innate immune genes (e.g., CD16) could counteract some of the age related decreases in the adaptive immune system (Cichon et al., 2003; Palacios et al., 2007; Peters et al., 2019). The age specific changes seen in our study could also represent observed changes in leucocyte compositions, gene expression profiles, and receptor density on specific cells (Palmer et al., 2006; Yung, 2000). Although a caveat, IgA and NKG2D display correlated gene expression, so it is possible the increase in both genes is indicative of a single biological response. As discussed above, our data contains animals within an age range of two years, which may limit our ability to observe nonlinear patterns of ageing and immunosenescence (e.g., threshold and quadratic effects), especially as the lifespan of devils has only recently halved, due to the impacts of DFTD (Hamede et al., 2015). For a clearer picture of age‐related changes and the extent of devil immunosenescence, multiple approaches measuring different aspects of cellular and protein functioning, correlation between blood and tissue gene expression, and cell counts would be necessary.

DFTD related changes in immune gene expression levels

Human studies have previously shown reduced expression of NKG2D in blood of patients with breast cancer measured with RT‐qPCR (Roshani et al., 2016) and gastric cancer measured with flow‐cytomtery (Saito et al., 2012), a pattern also observed by us in Tasmanian devils affected by DFTD. While future studies are needed to determine exact mechanisms in Tasmanian devils, human studies show NKG2D proteins being involved in cancer immunosurvailance, protecting hosts from tumour initiation by detecting ligands found on transformed cells (Guerra et al., 2008). However, tumours use downregulation of NKG2D receptors to avoid immune recognition (Liu et al., 2019; López‐Soto et al., 2015) leading to lower expression of NKG2D in the blood of cancer patients (Roshani et al., 2016; Saito et al., 2012). Future studies examining NKG2D ligands and measuring NKG2D expression using flow cytometry could verify whether Tasmanian devils are experiencing downregulation of NKG2D via DFTD initiated mechanisms.

DFTD infection also altered the expression pattern of some of the immunoglobulins (IgE, IgG and IgM:IgG) across seasons, that is, IgE was downregulated in autumn in healthy animals and during winter in DFTD affected devils; as well as IgM:IgG ratio was significantly lower in autumn in sick devils. These patterns may indicate a shift in B cell composition in DFTD infected animals, and thus, altered response to DFTD progression (Ravindranath et al., 1998) and susceptibility to other pathognes (Møller et al., 2003) during the colder months. DFTD infection also altered the gene expression of IgG between the sexes. Interestingly, in contrast to humans, where females often have higher B cell and immunoglobulin counts than males (Klein & Flanagan, 2016), there was no difference in the expression of IgG between healthy female and male devils, but DFTD infected females had lower IgG expression levels. Compared to DFTD infected males, female devils have been shown to maintain body condition with small tumours (Ruiz‐Aravena et al., 2018), and most tumour regressions have been observed in females (Pye, Hamede, et al., 2016). Consequently, the lower female IgG expression observed in the present study was surprising. As immunoglobulins, including IgG, have been been linked to initiating immune responses against malignant and transformed cells, and to stimulating and recruiting other immune cells to the site of infection (Díaz‐Zaragoza et al., 2015; Frazao et al., 2019; López‐Soto et al., 2015; Palma et al., 2018; Schwartz‐Albiez et al., 2009), a potential explanation (although admittedly speculative and not supported here by data from cell counts) is that females that are able to maintain IgG expression levels are less impacted by DFTD (e.g., slower tumour pregression, lower susceptibility to DFTD etc). Future studies investigating antibody levels throughout the progression of DFTD would be necessary to validate these observations.

A consistent finding of our study was the strong positive correlation between CD16 expression levels and DFTD status, independent of any other factors. CD16 is known to activate the ADCC pathway when binding to IgG molecules and is expressed on the surface of NK cells, neutrophils and a subset of monocytes in blood (Yeap et al., 2016). Therefore, the observed high expression of CD16 could be either an increase in sensitivity ‐ higher number of CD16 receptors per cell (Yung, 2000; Yung & Mo, 2003) ‐ or an increase in the number of cells expressing the CD16 receptor in DFTD affected animals. As NK cells target tumour cells via the ADCC pathway (Arnould et al., 2006; Roberti et al., 2012) and CD16 is known to be expressed on NK cells, our findings could potentially indicate an increased NK cell activity in DFTD infected devils. However, we did not see associations between the expression of CD16 and NKG2D, the other gene that is also expressed on NK cells. To support the involvement of NK cells in responses to DFTD, we expected to observe similar patterns in the expression levels of both CD16 and NKG2D. Interestingly, previous studies also did not find NK cell reaction, but identifyied the ADCC pathway being active in responses to DFTD (Brown et al., 2011). Since CD16 can activitae the ADCC pathways by being expressed on neutrophils and on a subset of monocytes, our findings, in concordance with Brown et al. (2011), indicate that other innate immune cells than NK cells might be generating the observed response to DFTD infection, as seen in human cancers (Uribe‐Querol & Rosales, 2015). In our study, it is not possible to distinguish which mechanism, or cell type is responsible for the increase in CD16 expression and it is also unclear if this is a positive or negative immune response. As Tasmanian devils have been shown to have higher neutrophil counts when infected with DFTD (Peck et al., 2016), these cells might be interesting to target in future studies. As the immunological outcome will be different if the increased CD16 expression results from a single cell expressing high level of CD16 or many cells expressing low levels of the protein (e.g., few assumingly preactivated cells vs. many inactive cells) further studies should investigate the change in the number of cells that express CD16 (Yung & Mo, 2003). Furthermore, longitudinal samples and targeted immunological approaches would be necessary to fully understand the significance of CD16 expression levels in DFTD epidemiology. As CD16+ monocytes have also been identified as potential early diagnostic marker in breast (Feng et al., 2011) and gastric cancers (Eljaszewicz et al., 2012), focusing on these pathways may allow the development of DFTD biomarkers.

CONCLUSION

Our study demonstrated the importance of including both intrinsic (e.g., sex and age) and environmental factors (e.g., seasons) when studying the response of endangered species, such as the Tasmanian devils, to pathogen challenges. The observed trend in CD16 expression in relation with DFTD infection highlights the importance of innate immune system in association with DFTD (that has so far been mostly overlooked), and hence could open an avenue for biomarker development. Our study shows that to improve the detection of disease related changes, larger sample sizes and sampling across seasons, age, sex, and environmental factors would be beneficial to identify differentially expressed genes between biological or experimental conditions. In addition future studies investigating associations between cell counts, antibodies, protein levels and gene expression profiles would assist in deciphering the underlying mechanisms and aid the interpretation of the patterns observed here. Although our study could not distinguish whether the observed associations between DFTD infection and immune gene expressions were (i) an immune response to the cancer, (ii) DFTD exploiting host resources, or (iii) the cancer modulating gene expression, the observed changes occur across multiple genes, suggesting the need for future research on the effect of DFTD on immune gene expression. Consequenly, future studies would benefit from using longitudinal samples from different populations to examine immune gene expression changes over time and devil survival in relation to DFTD infection. Our results highlight the need to incorporate sex and age specific differences, as well as seasonal variations in immune gene expression profiles when conservation management strategies on this iconic marsupial are considered and developed.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

AUTHOR CONTRIBUTIONS

Nynke Raven, assisted in research design, performed research, analysed data, wrote, reviewed and edited the manuscript, secured research funding. Marcel Klaassen, assisted in research design, assisted in statistical support and knowledge, reviewed and edited the manuscript. Thomas Madsen, assisted in research design, reviewed and edited the manuscript. Frédéric Thomas, assisted in research design, reviewed and edited the manuscript, secured research funding. Rodrigo Hamede, assisted in research design, collected samples and field data, reviewed and edited the manuscript, secured research funding. Beata Ujvari, assisted in research design, assisted in laboratory support and knowledge, assisted writing and structuring the manuscript, reviewed and edited the manuscript, secured research funding.

Supplementary Material

Click here for additional data file.