A novel multi-variate immunological approach, reveals immune variation associated with environmental conditions, and co-infection in the koala (Phascolarctos cinereus)

External signs of disease are frequently used as indicators of disease susceptibility. However, immune profiling can be a more effective indicator to understand how host responses to infection may be shaped by host, pathogen and environmental factors. To better inform wildlife health assessment and research directions, we investigated the utility of a novel multivariate immunophenotyping approach examining innate and adaptive immune responses in differing climatic, pathogen co-infection and demographic contexts across two koala (Phascolarctos cinereus) populations in New South Wales: the Liverpool Plains (LP), and Southern Highlands to South-west Sydney (SHSWS). Relative to the comparatively healthy SHSWS, the LP had greater and more variable innate immune gene expression (IL-1β, IL-6), and KoRV transcription. During extreme heat and drought, koalas from the LP displayed upregulation of a stress pathway gene and reduced adaptive immune genes expression, haematocrit and plasma protein, suggesting the possibility of environmental impacts through multiple pathways. In those koalas, KoRV transcription status, Chlamydia pecorum infection loads, and visible urogenital inflammation were not associated with immune variation, suggesting that immune markers were more sensitive indicators of real-time impacts than observed disease outcomes.


Description of study populations
The LP study site comprised a 15 km radius from − 31.159855,150.118456, within the Liverpool Plains of northern NSW, where koalas exist across a fragmented landscape consisting of agricultural and other privately owned land 62 .This population faces periodic climatic pressures such as heat waves and drought 63 , as well as a high prevalence of chlamydial disease 64 associated with infertility and population decline 65 .The other population: SHSWS, is located across the Southern Highlands within the Wingecarribee Local Government Area and Southwest Sydney within the Wollondilly Local Government Area.In contrast, this population contains a variety of habitat types ranging from dense continuous forest through to forest-urban landscape and a lower prevalence of chlamydial infection and clinical disease and greater fecundity have been recorded (unpublished data).

Sample collection
Samples were collected over 2019-2021 as part of larger ongoing population health studies.Within the LP, 108 samples were collected from 90 koalas across three sampling time points: time point A (July of 2019, n = 38), time point B (January 2020, n = 26) and time point C (September 2020, n = 44) (Table 1).This included 33 repeated sample sets of 15 koalas.Within the SHSWS, 41 koalas were sampled across two time points: time point D (September 2019, n = 9) and time point E (February, March and April 2021, n = 32), with no repeated sampling of individuals (Table 1).Mean monthly maximum and minimum temperatures for the month prior to sampling, and the average rainfall (mm) in the three months prior to sampling, were obtained from landholder data for the LP site as well as the Bureau of Meteorology (BOM) 66 (Table 1).Wild koalas were captured from trees using the 'noose and flag' technique 67 and sedated with alfaxalone (Alfaxan, Jurox; 1.8-2mg/kg) for sample collection.The following were recorded: Body condition score (BCS), ranging from 1 (emaciated) to 5 (excellent); sex; body mass (kg); and age class (A to I), estimated by tooth wear as in Gordon, (1991) 68 .Among females, the presence of back or pouch young was recorded.At each sampling, swabs were collected from each conjunctiva and from the urogenital sinus or penile urethra.However, only the urogenital swabs were considered in our study due to the low prevalence of ocular chlamydial infection (6.6% prevalence across time points A, B and C and 4.3% prevalence across time points D and E).The severity of chlamydial clinical signs at the urogenital site was graded by a Wet Bottom Score (WBS) based on Griffith, (2010) 69 ; each koala was graded across a scale of 0-10, whereby we considered a grade of 0-0.5 to represent an absence of clinical disease and grades 1-10 to represent clinical disease.WBS grading was performed by the same researchers across sampling field trips, with reference to a photographic standard and descriptive criteria.Approximately 1ml of blood was collected into an EDTA tube for measuring packed cell volume (PCV), total plasma protein (TPP) and haematological analysis.To measure PCV, a microhematocrit capillary tube was filled to approximately 80% of its length with EDTA blood, sealed with plasticine and centrifuged (12,000 rpm) for 3 min.The PCV was determined by placing the capillary tube on a microhematocrit line reader.TPP was determined using a refractometer.To determine the absolute counts of neutrophils and lymphocytes (× 10 9 /L), blood smears were made from the EDTA blood and 100µl of EDTA blood separated into a 1.5ml eppendorf, fixed in Streck Cell Preservative (Streck Cell Preservative, Omaha, United States) at a 1:10 ratio, and stored at 3°C for subsequent differential analysis in the lab.Total white cell counts were determined by Sysmex Hematology Analyzer (XN-1000 series) and differentials were determined from the blood smears by light microscopy.Within 8 h of collection, additional EDTA blood was centrifuged (4000 × g), and the white cell fraction (buffy coat) aspirated, stored in liquid nitrogen and then subsequently fixed in RNAlater™ solution (Invitrogen™) to minimise risk of freeze/thaw artefacts and maximise consistency with SHSWS samples.Within the SHSWS, EDTA Buffy coats were only collected from those animals sampled in 2019 and EDTA whole blood from animals collected in 2021.Of the EDTA whole blood, 200 µL was fixed in RNAlater™ solution for extraction of RNA (Invitrogen™) and stored at − 20°C until processing.

Extraction of RNA and analysis of immune, co-infection and stress markers using NanoString technology
The RNA from buffy coat and whole blood samples was extracted using the RiboPure™ RNA Purification Kit, Blood (Invitrogen, #AM1928).Briefly, each sample was centrifuged, the RNAlater supernatant aspirated, and the RNA extracted from the cell fraction and stored at − 80°C until further processing.A spectrophotometer (Nanodrop 1000) was used to quantify nucleic acid in each sample prior to analysis by a custom NanoString PlexSet at the Ramaciotti Centre for Genomics, University of New South Wales, Sydney, Australia.The following numbers of samples were not screened by NanoString due to inadequate RNA quantity: time point A, 4/38; Table 1.Sample sizes of koalas in each time point analysed by PCA and average maximum, minimum temperatures (°C) of the month prior to sampling and the average rainfall (mm) of the preceding 3 months prior to each time point across time points A, B and C in the LP and D and E in the SHSWS.; and four housekeeping reference genes representing high to medium expression levels (GAPDH, ACTB, Stx2, Nckap1) 71,72 .To account for any bias or differences that might arise between buffy coat and whole blood samples, the expression ranges of all 35 targets including housekeeping genes were validated between the buffy coats and whole blood samples from time points D and E of the SHSWS population, with no differing expression identified between sample types.

Bioinformatic analysis of NanoString raw data
The raw Fastq RCC files generated from NanoString were analysed in NanoString NSolver 4.0 software (NanoString Technologies).The raw copy number counts of mRNA of targets were normalised against the housekeeping genes, without a threshold set.Normalisation allowed for comparison of the relative expression of targets between samples.Briefly, normalisation was performed in NSolver using the following formula (Housekeeping average geometric mean of the row / Housekeeping average geometric mean of individual well) to generate a unique normalisation factor for each individual well.This normalisation factor was then multiplied by the raw counts of each target for each animal.Following normalisation, targets with > 50% of mRNA counts above detectable limits, across all samples and all 5 time points, were included in further analysis.Among these, samples with mRNA counts below the Limit of Quantification (LOQ) for the assay were still included in the analysis as the error likely within this range was considered negligible in comparison to the amplitude of variation overall.Some rejected targets were subsequently evaluated for the full sample set using RT-qPCR if previously published assays were available (see below).

Real-time qPCR of immune targets: CD4, CD8β, IFN-y, IL-10 and IL-17A
Due to limited detection by NanoString, the relative expression of CD4, CD8β, IFN-γ, IL-10 and IL-17A, were determined using previously developed and validated qPCR assays 71,[73][74][75] .The target GAPDH 71 , was included as a housekeeping gene.Briefly, for removal of genomic DNA, the extracted RNA was treated using the DNase I, RNase-free kit (1U/µL) (Thermo Scientific™, #EN0521).cDNA synthesis was then performed on the DNase treated RNA using the RevertAid First Strand cDNA synthesis kit (Thermo Scientific™, #K1622) with a Revertaid negative included for each sample to confirm an absence of genomic DNA.PCR efficiencies for each assay were between 90-100% and the inter-assay variance was less than 5% across all targets.Reactions were made to a final volume of 20 µL consisting of 10 µL of SYBR Green Supermix (SsoAdvanced™ Universal SYBR® Green Supermix-BIORAD, 300 nM of each primer, 6.8 µL of dH 2 O and 2 µL of DNA with samples run in duplicates.PCR cycling conditions were a two-step PCR of 40 cycles consisting of an initial denaturation (3 min at 95°C), denaturation (15 s at 95°C) and combined annealing and extension step (30 s at 59°C) for 40 cycles with a melt curve between 55°C-95°C with 5 s stops.The relative expression of each target was then determined using the ΔCt method (Cq average gene of interest ratio − Cq average housekeeping gene ratio) 76,77 .For more intuitive interpretation, the ΔCt was then inverted to 1/ΔCt to make a high value indicate high expression and a low value indicate low expression.Samples with an absence of quantification for the qPCR were excluded from further analysis.

Quantification of urogenital shedding of C. pecorum and PhaHV-1
DNA was extracted from urogenital swabs using the MagMAX™ CORE Nucleic Acid Extraction Kit (Ther-moFisher Scientific).To determine the chlamydial shedding at the time of sampling, an established probe-based multiplex qPCR assay was used, targeting the beta-actin mRNA as a housekeeping gene, chlamydial genus 23S rRNA and C. pecorum ompB gene 78 .C. pneumoniae was not tested for detection in the swabs given its low prevalence on a population level 79  www.nature.com/scientificreports/

Principal components analysis (PCA) and immunological changes between populations
To visualise initial distributions and patterns of expression of the immune and stress genes and co-infection loads across time points A, B, C, D and E, box and whisker plots were created using R 81 .The variables that yielded quantifiable values above the detection limits in more than 50% of the samples (Nanostring, 9/21 immune variables; qPCR, 4/5 variables), were selected for inclusion in the PCA (CD3, CD79b, IL-1β, IL-6, IL-8, PhciDAB, PhciDBB, 1/ΔCD4, 1/ΔCD8β, 1/ΔIFN-γ, 1/ΔIL-10), along with absolute counts of neutrophils and lymphocytes (× 10 9 /L), to ascertain whether the variation in cytokine expression, especially low cytokine gene expression relative to housekeeping genes, might be attributed to elevated neutrophil counts.Before conducting the Principal Components Analysis (PCA), a Spearman's correlation analysis was performed in R 81 , to investigate potential correlations among the variables with all the continuous variables to be input in the PCA.A total of 108 observations across the LP and SHSWS with no missing data across the 13 immune variables were used in the analysis, which included LP time point A (n = 28), time point B (n = 15), time point C (n = 29) and SHSWS time point D (n = 7), time point E (n = 29) (Table 1).Proportion of females were: time point A (17/28), time point B (6/15), time point C (14/29), time point D (5/7) and time point E (12/29).No females were identified to be carrying joeys across time points A, B and C but young were identified in 2/5 females in time point D and 4/12 females in time point E. Within the LP, for individuals sampled at multiple time points, only data from one time point was used in the PCA.Selection of time points for inclusion was dependent on creating the best balance of sample sizes across time points A, B and C. For example, time point B had had a smaller sample size compared to time point A and C and therefore the data within this time point was selected for over the other time points.To visualise the grouping of the immune variables between populations, across time points A, B, C, D and E and sex, biplots for each of these parameters were created from the Principal Components (PCs) explaining the most variation for both PCAs.The PCs that explained the most variation, were subsequently extracted and analysed by a one-way analysis of variance (ANOVA) with the emmeans package 82 , to identify whether PC values significantly differed across time points between the populations.

PCA and immunological changes across timepoints within the LP population
To determine if immune profiles differed among time points within a population, when isolated from the effect of significant between population differences, a second PCA investigating the same 13 immune variables within the LP population across time points A, B and C was performed.As above, prior to performing the PCA, a Spearman's correlation was performed to confirm that correlations were present among the continuous immune variables within the population.This PCA included a total of 72 observations from time point A (n = 28), time point B (n = 15) and time point C (n = 29).As above, the PCs that explained the most variation were extracted and analysed by a one-way analysis of variance (ANOVA) and the emmeans package 82 , to identify whether PC values significantly differed among time points.

Statistical analysis: repeated measures multivariate analysis of variance (MANOVA) and Paired t-test of recaptured koalas over time points B and C
The expression of the innate and adaptive responses across the longitudinally sampled events of koalas over time points A, B and C in LP was analysed by a repeated measures multivariate analysis of variance (MANOVA) and paired t-tests using R 81 .This analysis was performed as a comparison to the PCA of immune variables over time points A, B and C, to determine whether the relative expression of immune targets within the cross-sectional dataset paralleled the longitudinal dataset.Due to a large quantity of missing data across variables in the 15 koalas recaptured over time points A, B and C in LP, a subset of 10 koalas with a more complete dataset from sampling events in time points B and C were analysed in the paired t-test analysis.Prior to performing the repeated measures MANOVA and paired t-test analysis, the data were checked for normality using Histograms, Kernel Density and q-q plots and the Shapiro-Wilk test.Log and Tukey transformation methods were applied to nonnormally distributed variables with the transformation method producing the closest to normal distribution of data being used.If the data for a variable were not normally distributed in one recapture event but normal in the other, the data from both recapture events were transformed.Two repeated measures MANOVAs were run, each incorporating innate and adaptive variables together based off the grouping of innate and adaptive variables from the PCA over time points A, B, C, D and E, to identify whether a significant difference in these two arms of the immune response existed across recapture events.Where the MANOVA were significant, post-hoc paired t-tests were performed for each individual innate or adaptive variable across time-points, using Bonferroni adjustment of p values to account for Type I errors.

Fold difference and change in expression of immune variables between time points A, B, C, D and E and MANOVA
Between time points A, B, C, D and E, the relative fold change in expression of immune variables that contributed the most to the PCs identified to significantly differ between time points, was calculated.For immune variables analysed by NanoString, the fold change for variables were calculated by dividing the mean expression in one time point by the mean expression in another time point.For immune variables analysed by qPCR, the 2 −ΔΔCt method 76,77 was applied using the mean ΔCt values for each time point.Between recaptured koalas over time points B and C in the LP, the relative fold changes of expression of individual innate and adaptive variables identified to be significantly different between time points from the MANOVA were calculated.The fold change in immune variables analysed by NanoString was performed as described above.For immune variables analysed by qPCR, the 2 −ΔΔCt method was used 76,77  www.nature.com/scientificreports/

Associations of immunological changes with co-infection, disease and environmental factors
To identify whether transcription of KoRV-A env, KoRV-D env, KoRVenvCKS17 and KoRV pol differed significantly among all time points A, B, C, D and E across both populations, a one-way ANOVA was used if variances of expression were equal across time points, followed by Tukey's HSD post hoc tests within the emmeans package 82 ; or a Welsh and Brown-Forsythe test if variances unequal.Across the populations, we were unable to explore associations between immune profiles with KoRV transcription, due to the bimodal distribution of KoRV expression between populations (Supplementary Figure S4), or with loads of infection with C. pecorum and PhaHV-1, or inflammation at the UGT site, due to the low prevalence of infection with these pathogens in the SHSWS relative to the LP population.We refrained from investigating correlations to age and stress due to substantial population differences between these groups, clear associations could not be established with respect to these factors.Associations between co-infection parameters, and clinical data such as urogenital inflammation, age and stress with immunophenotype were performed within the LP alone to remove the effect of between population differences.Among time points A, B and C in the LP population, variation of KoRV-A env, KoRV-D env, KoR-VenvCKS17 and KoRV pol transcription was evaluated by ANOVA if variances of expression were equal across time points followed by Tukey's HSD post hoc tests within the emmeans package, or a Welsh and Brown-Forsythe tests if variances were not equal.Individual linear regression models were performed to determine the relationship between the PC values of those variables that explained the most variation in the PCA, and transcription of KoRV-A env, KoRV-D env, KoRVenvCSK17 and KoRV pol, and mucosal loads of C. pecorum (shedding) and PhaHV-1.To determine associations between the PC values and the presence and absence of urogenital C. pecorum, WBS and age, independent t-tests were used.Prior to running the t-tests, assumptions of normality and homogeneity of variance of the PCs were checked and data transformed if necessary.Presence/absence of urogenital C. pecorum shedding was classed as 'Shedding' and 'Not shedding' , respectively.WBS was grouped into two outcomes: 'No clinical disease' (WBS 0-0.5) and 'Clinical disease' (WBS 1-4) and included only those individuals identified to be shedding C. pecorum at the urogenital site at the time of sampling to ensure only infection-challenged individuals were included.As immune responses can vary across age groups, particularly in older individuals 10 , t-tests included age groups representing 'adult' (classes D-E) and 'old' koalas (classes F-I) with juvenile koalas (classes A-C) excluded due to inadequate sample size.PCs were first evaluated for normality and homogeneity of variance, and transformation performed if necessary.
One-way ANOVAs followed by Tukey's HSD post hoc tests using the emmeans package 82 , was performed to identify differences in PCV, TPP and transcription of FKBP5 across time points A, B and C if variances of expression were equal.Individual linear regression models were then performed to determine the relationship between the PC values of the PCs which explained the most variation in the PCA, and PCV, TPP, and transcription of FKBP5.Linear regression models were also run to determine whether there was any relationship between transcription of FKBP5, and PCV or TPP.

Detection and quantification of loads of co-infecting pathogens: KoRV, C. pecorum, PhaHV-1 and trypanosomes
By NanoString, KoRV-A env, KoRV-D env, KoRVenvCKS17 and KoRV pol were detected in all koalas across both populations; KoRV-B env gene transcripts were not detected in any sample from either population.Urogenital ).Urogenital swabs were not available for one animal in the LP and therefore it was excluded from these analyses.Among the trypanosome species only T. irwini was detected, in 15/41 koalas from the SHSWS; no trypanosome RNA was detected in samples from the LP and so trypanosomes were not considered for further analysis.

PCA and immunological changes between populations
Box and Whisker plot visualisation of the counts of neutrophils and lymphocytes (× 10 9 /L) across all time points A, B, C, D and E, indicated no differences across time points (Supplementary Fig. S1).Across all time points A, B, C, D and E, the Spearman's rank correlation matrix identified correlations between a substantial number of immune variables.The first two PCs explained a total of 46.9% of the variance with PC1 explaining 26.7% and PC2 explaining 20.2%.Here PC1 represented the innate response with variables IL-8, IL-6, IL-1β, PhciDAB, PhciDBB contributing the most to this component (Supplementary Table S1).PC2 mostly represented the adaptive response with variables CD3, 1/ΔIFN-γ, 1/ΔCD8β, CD79b, 1/ΔCD4, and neutrophils, contributing the most (Supplementary Table S1).Biplots suggested differences between the time points in LP and SHSWS populations (Fig. 1

PCA and immunological changes across timepoints within the LP population
Within the LP alone, across time points A, B and C, the first two PCs explained a total of 48.6% of the variance, with PC1 explaining 26.8% and PC2 explaining 21.8% (Supplementary Fig. S2).Here PC1 represented an increase in the adaptive response and a decrease in the innate response with PC2 not representing clear adaptive or innate pathways.In order of contribution, variables 1/ΔCD8β, IL-1β, 1/ΔCD4, IL-6, IL-8 and neutrophils contributed the most to PC1 with adaptive variables 1/ΔCD8β, 1/ΔCD4 having negative loadings and innate variables IL-1β, IL-6, IL-8 and neutrophils having positive loadings (Supplementary Table S2).The remaining variables did not contribute as strongly to this PC (Supplementary Table S2).There was a significant difference in adaptive-innate balance (PC1) among time points A, B and C (PC1: F (2, 69) = [5.433],p = 0.006) with a significantly lower expression of adaptive loci in time point B relative to time point A (t (69) = 3.291, p = 0.004).No significant differences between time point A to C (t (69) = 1.217, p = 0.447) or between time point B to C (t (69) = − 2.297, p = 0.063) were identified.
The fold changes in the relative expression of immune variables contributing the most to PC1 and PC2 identified to be significantly different between time points A, B, C, D and E, and adaptive immune variables in the MANOVA across recapture events in time points B and C, are presented in Supplementary Fig. S3 and Table S3, S4. www.nature.com/scientificreports/

Associations of immunological changes with co-infection, disease and environmental factors
Box and Whisker plots indicated marked similarity in distributions for transcription of KoRV (particularly KoRV pol), IL-1β and IL-6 within the LP and SHSWS populations and these differed markedly between these populations (Supplementary Fig. S4).Expression of IL-1β and IL-6 and transcription of KoRV, were greater and more variable in the LP samples compared to the SHSWS (Supplementary Fig. S4).Using the one-way Welch and Brown-Forsythe tests due to differences in variance, significant differences in transcription of KoRV-A env (F ( 4

Discussion
This is the most extensive study investigating both innate and adaptive immune responses together with climatic, pathogen and host data in the koala.This study demonstrates the utility of a multi-variate immunological analysis to evaluate population health in a wildlife species, across diverse geographical, climatic and co-infection contexts.As an association-based study, the aim was to determine how, and to what degree, immune profiles of koalas changed under these differing conditions.Our multivariate approach covering both the innate and adaptive arms of the immune system, greatly facilitated interpretation of observed changes in immunophenotype and their association with concurrent co-infection and climatic events in a way that would not have been possible using single variate immune markers and disease outcome as measures of host susceptibility.The LP population, a population experiencing drought and a greater retroviral expression and chlamydial disease, had a significantly greater and more variable expression of innate immune parameters compared to the SHSWS.Marked similarity in distributions of KoRV transcription and the innate immune markers IL-1β and IL-6 identifies these immune mediators as candidates for investigation in mechanistic studies of KoRV pathogenesis.Within the LP, significantly greater expression of the stress co-chaperone, FKBP5, and reduced PCV, TPP and transcription of adaptive immune markers coincided with a sustained heat event but limited evidence of correlation between these parameters within individuals, suggests that this variation may arise from multiple mechanisms.Urogenital shedding of C. pecorum, PhaHV-1 and urogenital infection status, inflammation and age were not associated with immunological changes, suggesting that immune variation was driven by processes other than, for example, chlamydial inflammation or old age.The outcomes of this study are a novel approach to investigate drivers of wildlife disease and more clearly identified associative relationships to direct mechanistic or replicated field studies to elucidate causative relationships.
The multi-variate immunological approach used in this study was novel in that it provided perspectives and interpretation of immunological variation across the different arms of the immune system that would not have been observed using single variables alone or using disease outcomes as indicators of susceptibility.Given the complexity of the immune response, variations or associations detected in individual immune variables may not be representative of overall pathways being affected, or determine which variations are biologically meaningful 6,18,20 .Past immunological studies of koalas, in relation to KoRV infection 12,15 , have shown variation in individual parameters but interpretation of relevance of these to the performance of the immune response as a whole has been challenging.The use of PCA facilitated a systems consensus approach that more clearly www.nature.com/scientificreports/associated changes with the innate or adaptive responses.The use of immunophenotype, including innate and adaptive immune function as an outcome, demonstrated immune changes in real time, which were not evident in external or visual manifestations of disease, due to its chronic and often permanent nature.In koala health surveys, the WBS grading system is utilised as an indicator of inflammation at the urogenital site 69 , however the approach is limited as it is only possible to observe outcomes in animals that have contacted the pathogen and, among those, some outcomes, such as reproductive disease, are not always clinically evident 5 .In addition, due to the nature of chlamydial disease, observed signs may represent either current inflammation or chronic fibrotic pathology of the bladder which may have occurred months to years prior.Use of immunophenotyping appeared to circumvent these issues.
The greater expression of innate immune genes and the similarity in distributions of KoRV, IL-1β and IL-6, highlight the degree of variation between koala populations in different contexts, though the mechanisms involved require further work to elucidate, due to multiple differences between the populations.Koalas within the LP have experienced and continue to endure more extreme climatic conditions such as heatwaves and drought 63 and also live across a moderately fragmented habitat 62 .At the time of the study, the SHSWS population experienced greater rainfall and more moderate ambient temperatures as shown in Table 1 66 , and koalas appear to live across more continuous habitat in the north of the population.The difference in KoRV expression between the populations might have been related to the north-south cline reported by Blyton et al. 34 and may or may not be influencing immunophenotypic variation.Alternatively, the differing immune status may have influenced KoRV expression, or the two may be linked by some yet unidentified common factor.Numerous studies of koalas have demonstrated associations between KoRV infection and variable expression of cytokines across Th1, Th2 and Th17-A pathways 12,15 .The similar distributions of innate cytokines IL-1β and IL-6 and KoRV pol transcription loads (Supplementary Fig. S4) are of interest considering a previous study in mice where increased secretion of IL-1β was induced by, and enhanced, retroviral infection, replication and persistence in vivo by recruiting susceptible cells to the sites of infection creating a positive feedback loop for retroviral infection 83 .Induction of other cytokines IL-10, IL-6 and IL-8 following infection with KoRV 84,85 and other retroviruses 86 , has also been described; we found associations with IL-6, some correlation with IL-8 but not IL-10.Given that these findings are based on associative relationships, further mechanistic studies should be performed investigating the directionality and causality of these associations.
Within the LP alone, no association between the decreased expression of the adaptive response to transcription of KoRV during the extreme weather event was identified, which suggests that within in this event, adverse environmental conditions were more important in driving immunological variation than KoRV was, though it appears likely this is manifesting through multiple pathways.The consistent differences between sampling events but weak correlation between expression of FKBP5, PCV and TPP and adaptive variables suggests complex interactions that may involve, for example, heat, hydration, nutrition or exposure to PSMs.The lower PCV (anaemia) and TPP values, and elevated expression of a stress receptor pathway gene FKBP5, occurring in time point B, are all consistent with expected impacts of extreme ambient temperatures.Anaemia has been identified in sheep 87 and humans experiencing prolonged heat exposure 88 .Erythropoiesis can be reduced under heat stress events due to an increased partial pressure of oxygen in the blood from increased respiratory rates of animals experiencing heat 87 .Although we did not identify any associations between elevated expression of FKBP5 and reduced PCV and TPP and the decreased adaptive response during this event, in humans, increased stress impacts immune function, specifically by increasing inflammation 89 and expression of innate cytokines: IL-1β and IL-6 [48][49][50] .Koalas in time point B were observed to be in poorer body condition, consistent with strategic decrease in food intake by endotherms in increased ambient temperatures to maintain their body temperature 90 and reduce diet induced thermogenesis 91 .This however may suppress immune investment, specifically adaptive responses 92 , as initiating and maintaining an immune response requires energy expenditure 93 and can elevate metabolic rates 93,94 .In a study of leukocyte profiles of bats across the Neotropics, bats at the edge of their distributions had increased variation and quantities of white blood cell counts, likely from chronic stress due to the greater variation in temperature and rainfall experienced at these range limits 17 .Greater ambient temperatures have also been suggested to increase production of secondary metabolites across a number of plant species 95 and PSMs have been shown to impede koala immune cell function in vitro at concentrations experienced in their plasma 96 .Increased CO 2 emissions as a result of climate change can also increase concentrations of PSMs such as phenol and tannins 97 .We did not measure PSMs, however we speculate that extended exposure of koalas to increasing PSMs as a result of environmental extremes which are predicted to occur more frequently within Australia 98 , may adversely impact innate and adaptive responses.These findings highlight the need for this study to be replicated across diverse populations, and coupled with mechanistic studies, to evaluate these alternative hypotheses and the thresholds at which these climatic extremes have significant impact on immune function.
It is unlikely that other measured variables such as population differences in the proportions of chlamydial infection and disease, sex, breeding season and age demographics, are associated with the differences in immune function observed as we found no associations in immune variation with these factors within populations.In comparison to the LP, the SHSWS population has lower prevalence of chlamydial infection and disease but, within the LP, immune profiles did not differ with respect to chlamydial shedding or disease.As a systemic immunophenotype was measured in this study, it is possible that it may not reflect local, mucosal, inflammatory responses in response to infection with Chlamydia at the UGT site, particularly as the spectrum of acute to chronic UGT disease, were grouped in this study together as 'Clinical disease' .We consider it reasonable that systemic challenges posed by varying environmental conditions are likely to exert greater effects on systemic immunophenotypes than would more localised mucosal challenges posed by chlamydiosis.Sex and breeding season are unlikely to have affected interpretation of the immunological variation as sampling time points B and E both occurred within the breeding season over November to February and still differed.Within each of the sampling time points, sexes were balanced.Although the upregulation of IL-6 expression observed within time www.nature.com/scientificreports/point B coincides with reports of increased expression of this cytokine in December 15 , the expression patterns of the remaining cytokines of the innate and adaptive responses did not correspond to previously reported changes across seasons.Despite that the SHSWS population had a greater proportion of younger animals, no significant differences in immune responses between adults to older koalas were identified within the LP, likely ruling out any age associated impacts.Overall, the associative relationships identified in our study will guide additional health studies in the koala to prove causative relationships between altered immune function with climatic conditions and co-infection with KoRV.Future studies should apply the multi-variate immunological approach used in our study, to a highly curated longitudinal dataset over multiple time points and populations in association with environmental and co-infection parameters to define causation and establish thresholds for detrimental impacts.In particular, to better define environmental and climatic conditions and their impacts, additional data such as diet composition, leaf moisture 43,99 and days since last rain 100 should be collected.Further development of assays for markers identified to be important in chlamydial infection, but which remained below the limit of detection in our study, in particular IL-17A, should be addressed and included for future studies.This may involve using mitogen stimulation assays to upregulate those markers below the limit of detection in baseline samples, or investigation of more abundant pathway targets such as receptor molecules.Including markers representing all arms of the immune system will facilitate more accurate interpretation of immunological variation 6,18,20 .

Conclusion
Our study provides a novel approach, in that we assessed variation across the various arms of the immune response, to better quantify immunological responses to external stressors, a useful tool in the context of climate change and co-infecting pathogens.Continued changes to Australia's climate are predicted to have adverse impacts on koala populations 43,101,102 .Koalas are listed by the IUCN as one of the top 10 species to be impacted by climate change 103 , but thresholds for detrimental impact on immune responses and disease susceptibility are yet to be identified for the species.Additional studies assessing immunological responses in koala populations experiencing climatic extremes are therefore needed.The distributions of upregulation of innate cytokines and increased KoRV transcription identified in our study is seminal work required for directing the focus and efforts for future mechanistic studies to identify the directionality of these relationships and how they may alter disease outcomes.The outcome would be an enhanced ability to predict habitat and population viability in the face of continued challenges placed by habitat degradation, climate change and co-infections.The approach used in our study is also transferable to other epidemiological studies of health in other wildlife species.

Total sample size of buffy coat/whole blood (n) Sample size (n) in PCA Number of females Month prior to sampling
time point B, 2/26; time point C, 4/44 and time point E, 1/32 samples.To represent innate and adaptive arms of the immune response, including both humoral and cell-mediated responses, as well as the concurrent stress and co-infection traits of the individuals, the NanoString panel included 35 targets comprising innate and adaptive immune markers (CD3, CD79b, CD4, CD8β, IL-10, IL-12A, IL-17A, IL-18, IL-1β, IL-22, IL-4, IL-6, IL-8, IFNy, MHCIUA, PhciDAB, PhciDBB, TLR2, TLR4, TLR7 and TNFα); co-infection targets (KoRV-A env, KoRV-B env, KoRV-D env, KoRVenvCKS17 and KoRV pol; and Trypanosome species T. irwini, T.gilletti, T.copemani); a stress pathway gene (FKBP5), that encodes FK506 binding protein 51, a co-chaperone protein and mediator of glucocorticoid receptor (GR) function and activity Vol:.(1234567890) Scientific Reports | (2024) 14:7260 | https://doi.org/10.1038/s41598-024-57792-7www.nature.com/scientificreports/ 80esent across both the 23S and ompB targets and shedding absent if amplification was not present across these targets.The average Ct value of the neat ompB target was subsequently used for further statistical analysis and if inhibitors were present, the 1:10 value was converted to a neat value based off a 3.3 cycle difference per tenfold dilution in a PCR with 100% efficiency.For the detection of PhaHV-1 in urogenital swabs, an established qPCR assay was used80.Briefly, samples were run in duplicate, and each PCR reaction was made up to final volume of 20 µL, consisting of 10 µL of SYBR Green Supermix (SsoAdvanced™ Universal SYBR® Green Supermix-BIORAD), 250 nM of each primer, 7 µL of dH 2 O and 2 µL of DNA.PCR cycling conditions were a two-step PCR of 40 cycles consisting of an initial denaturation (3 min at 98°C), denaturation (10 s at 98°C) and a combined annealing/extension (30 s at 56°C) with a melt curve of 81°C indicating a specific product.A positive control (synthetic plasmid control at 1 × 10 3 concentration) and No Template Control (NTC) were included in each PCR run with PCR efficiencies ranging between 90-100% and the inter-assay variation being less than 5% across both targets.The average Ct value for PhaHV-1 was used for further statistical analysis.Vol.:(0123456789) Scientific Reports | (2024) 14:7260 | https://doi.org/10.1038/s41598-024-57792-7