Immune and sex-biased gene expression in the threatened Mojave desert tortoise, Gopherus agassizii

The immune system of ectotherms, particularly non-avian reptiles, remains poorly characterized regarding the genes involved in immune function, and their function in wild populations. We used RNA-Seq to explore the systemic response of Mojave desert tortoise (Gopherus agassizii) gene expression to three levels of Mycoplasma infection to better understand the host response to this bacterial pathogen. We found over an order of magnitude more genes differentially expressed between male and female tortoises (1,037 genes) than differentially expressed among immune groups (40 genes). There were 8 genes differentially expressed among both variables that can be considered sex-biased immune genes in this tortoise. Among experimental immune groups we find enriched GO biological processes for cysteine catabolism, regulation of type 1 interferon production, and regulation of cytokine production involved in immune response. Sex-biased transcription involves iron ion transport, iron ion homeostasis, and regulation of interferon-beta production to be enriched. More detailed work is needed to assess the seasonal response of the candidate genes found here. How seasonal fluctuation of testosterone and corticosterone modulate the immunosuppression of males and their susceptibility to Mycoplasma infection also warrants further investigation, as well as the importance of iron in the immune function and sex-biased differences of this species. Finally, future transcriptional studies should avoid drawing blood from tortoises via subcarapacial venipuncture as the variable aspiration of lymphatic fluid will confound the differential expression of genes.


Introduction
The innate and adaptive immune systems of non-avian reptiles remain a challenge to characterize, particularly in how they function relative to avian and mammalian model systems [1]. Challenges to this characterization are fourfold. First, while the number of available reference genome assemblies for non-avian reptiles is increasing (e.g., [2][3][4][5][6][7]), these resources are newer and less well-developed than those of model systems, such as human, chicken, and African clawed frog [8][9][10]. As a consequence, it is still largely unknown how much of the "immune gene set" of human, chicken, and frog is conserved in non-avian reptiles. Second, there are fewer functional studies about the immunological mechanisms and pathways governing how non-avian reptiles respond to pathogens compared to mammal and bird systems [1]. Third, the closest-related functional models-human and chicken-are endotherms so their immune processes occur at a consistent internal temperature unlike non-avian reptiles whose ectothermy means they have a wider set of physiological and activity states. The biochemistry that mediates immune responses may range in efficiency under these different conditions, adding a layer of complexity to their immune function [11,12]. Finally, how the sex of an organism modulates immune functions is now well-documented in mammals [13], but how sex-biased differences manifest in the immune function of non-avian reptiles with temperature-dependent sex determination is less known [14]. In many species, males and females exhibit differences in anatomy, morphology, physiology, and behavior. To some extent, sexually dimorphic traits are the product of sex differences in gene expression, which allow phenotypic differences from a common autosomal genome. Genes differentially expressed by sex are known as sex-biased genes and can be female or male biased, depending on which sex exhibits higher levels of transcription [15,16]. Sex-biased gene expression is the result of differential gene regulation between males and females, but in taxa with specialized sex chromosomes, gene dosage also plays a role in unequal transcription. Sex affects immune function as well as disease prevalence and severity. In general, females mount stronger innate and adaptive immune responses than males [17,18]. Possible explanations for this are differential regulation of gene expression by sex hormones [19,20] or differences in behavior [17,21]. Many studies focus on sex-biased expression of immune genes on sex chromosomes, but the majority of immune genes are located on autosomal chromosomes and these warrant further study [22]. How temperature-based sex determination affects sex-biased gene expression is also understudied. Therefore, the role of environment on disease susceptibility and prevalence, as well as how these effects manifest differently by sex are topics of considerable interest that have broad implications for the management of wildlife.
One way to fill in these knowledge gaps is to assay the transcriptional response of how a non-avian reptile responds to infection where the host-pathogen relationship is reasonably well understood. Infection and disease of the Mojave desert tortoise (Gopherus agassizii) by the pathogenic bacteria Mycoplasma spp. is among the most extensively characterized in chelonians [23]. Furthermore, because desert tortoises are long-lived animals, the cumulative effects of exposure to stressors such as prolonged infection may have considerable long-term impact for these animals. Knowledge gained about transcriptomic response to infection in this system will shed light on the innate and adaptive immune response of non-avian reptiles to infection broadly, including sex-biased effects.
Learning about this host-pathogen relationship is a major conservation priority for the desert tortoise, which was listed as threatened under the US Endangered Species Act in 1990 [24]. In this host-pathogen relationship, the Mycoplasma agassizii bacteria disrupt the tissue of the ciliated mucosal epithelium leading to upper respiratory tract disease (URTD). Although this disease has been extensively studied in desert tortoises, it remains unclear if or how M. agassizii bind to epithelial cells or if they just reside in the mucosal layer similar to other respiratory microbes [25]. In desert tortoises, URTD can lead to severe damage to the tissues of the upper respiratory tract and occlusion of nasal passages by thick mucosal discharges. The URTD is also thought to be a direct cause of mortality in Mojave desert tortoises [23,[26][27][28].
Higher pathogen loads from M. agassizii generally correlate with increased clinical signs and adaptive antibody responses [29,30]; however, antibody production can be delayed in this species by months to years after M. agassizii infection [30,31]. By the time most tortoises are classified has having URTD and an antibody response is detected, individuals have likely been infected for a long period of time. For some time initially after infection by pathogen, tortoises may not show clinical signs or an antibody response but would test positive for M. agassizii by qPCR [32][33][34]. As many of the recovery unit populations remain below viability levels [35], it is important to the health and viability of the species to understand how these bacterial infections and URTD impact the health of tortoises.
To learn more about the transcriptional response to this host-pathogen relationship, we used RNA-Seq to analyze the blood-based gene expression patterns of male and female tortoises with severe M. agassizii infection, tortoises inoculated with M. agassizii, and uninfected tortoises. Based on previous studies, we expect cytokines (e.g. IFN-γ, interleukins, TNFα1) and inflammatory signaling pathway genes to be expressed at higher levels in tortoises with bacterial infection, as they are important mediators of a host's defense to pathogens. Since the tortoises in this study were infected for long periods of time, we also expect signs of chronic infection and humoral immune response such as immunoglobulins or lymphocyte-specific cluster of differentiation genes, which are necessary for targeting, activation, and survival.

Experimental design
We analyzed 25 tortoise individuals across three experimental groups with varying degrees of bacterial infection including: 1) individuals discovered with severe M. agassizii infection and were not inoculated from experimentation (Severe Infection group, SI); 2) individuals inoculated with M. agassizii showing medium infection (Medium Infection group, MI); and 3) individuals serving as a control without any known infection (No Infection group, NI; Table 1). In our study, both medium and severe infection groups were sampled from captive colonies with documented infection for multiple years. We were not able to incorporate subclinical animals with early or low levels of infection in this study. All handling and experiments using animals were approved by the U.S. Geological Survey-Western Ecological Research Center Animal Care and Use Committee (WERC 2012-03) and covered under state (Nevada Division of Wildlife Permit #S33762) and federal (U. S. Fish and Wildlife Service TE-030659) permits.
For the severe infection (SI) group, we chose captive adults (N = 9; 6F:3M) from the Desert Tortoise Conservation Center (DTCC) in Clark County, Nevada, USA (35.975256, -115.251048) that were classified with severe infection based on long-term health evaluations by experienced veterinarians. Each tortoise in this category had a confirmed long-term M. agassizii infection and multiple clinical signs of potential illnesses associated with long-term weight loss and reduced or under-conditioned body condition (Table 1). Due to their poor overall health, consistent with captive herd management protocols and based on veterinary guidance, most tortoises (8 of 9) were euthanized following sample collection and immediately necropsied to evaluate tissue conditions morphologically and histologically. Tortoises were euthanized after this study period by licensed veterinarians using a mixture of ketamine (5 mg/kg) and dexmedotomidine (0.1 mg/kg) injected intramuscularly as an anesthetic. Once animals were non-responsive, Euthasol (2 ml/kg) was injected intravenously into the subcarapacial cranial plexus [37].
For the medium infection group, we also used captive adult tortoises (N = 7; 1F:6M) from the DTCC that were experimentally exposed to M. agassizii as part of a previous study [29]. These tortoises tested positive for the presence of M. agassizii bacteria for four years and exhibited targeted immune responses (specific antibody production measured using enzyme-linked immunosorbent assay (ELISA) tests) to M. agassizii as well as intermittent clinical signs associated with inflammatory responses to this infection for two years prior to sampling [31]. For the control group we chose clinically normal, adult tortoises (N = 9; 5F:4M) from a wild population that has been monitored since 2006 in Hidden Valley, Clark County, Nevada, USA (36.528008, -114.975905). Tortoises in this control group were clinically normal based on visual examination by veterinarians and tortoise biologists, and each tortoise had been evaluated and assessed as clinically free of M. agassizii infection for 11 consecutive years [38,39] prior to collecting samples for this study. Wild tortoises were not euthanized.
All tortoises were assessed and sampled in peak summer (July-early August) between 0500-0800 hours to minimize circadian and seasonal influences on measured blood analytes. Due to logistical constraints tortoises were sampled during the same season but in different years; Medium Infection (MI) tortoises were sampled in 2017, Severe Infection (SI) tortoises in 2013, and No Infection (NI) tortoises in 2015.

Choice of blood and venipuncture site
In nonmammalian vertebrates, whole blood is appropriate for gene expression studies for two reasons. First, white blood cells, which include granulocytes, monocytes, and lymphocytes,   . Numerical body condition scores (BCSs) were used to assess overall muscle condition and fat stores with respect to skeletal features of the head and limbs [36]. BCS scores were categorized as 'under (1-3),' 'adequate (4-6)' or 'over (7-9)' condition. https://doi.org/10.1371/journal.pone.0238202.t001 allow for the molecular characterization of the host's systemic response to mycoplasma infection. Second, in reptiles erythrocytes are nucleated and transcriptionally active [40,41] and thrombocytes remain as intact cells instead of producing anuclear platelet cytoplasmic fragments [42]. Genes involved in insulin signaling, electron transport chain, stress, and oxidative response are shared between whole blood and liver [41], making whole blood a non-invasive sample to assess immunological and physiological response to infection. We extracted~2.5 mL whole blood from all tortoises via jugular venipuncture [43] using either a 1.91-cm, 25-gauge needle-IV infusion set and 3 mL syringe, or subcarapacial venipuncture [37] using a 3.81-cm, 23-gauge needle and 3 mL syringe. All MI and NI tortoises were sampled using subcarapacial venipuncture while SI tortoises were sampled using both subcarapacial and jugular venipuncture (N = 9, 4 subcarapacial:5 jugular). Syringes were coated in sodium heparin to prevent coagulation and blood was collected from severely infected tortoises prior to euthanasia. Two aliquots of whole blood were made per sample. The first aliquot (0.1-0.5 mL blood) was placed immediately into RNeasy1 Animal Protect collection tubes (Qiagen, Valencia, CA) for RNA sequencing and gene expression analysis. The second aliquot of~1.5 mL whole blood was placed in BD Microtainer1 tubes with lithium heparin in order to assay blood counts, hematology, blood chemistry, trace elements, and vitamin A concentration (described in [38]). Samples were stored on wet ice for no more than four hours. We separated plasma from the second aliquot using centrifugation (1318 x g for 10 minutes) and stored in an ultracold freezer (-70˚C) until further processing. Aliquots of plasma (0.01 mL) were screened for antibodies specific to M. agassizii using an enzyme-linked immunosorbent assay (ELISA; [44]). Sloughed epithelial cells were also collected using sterile oral swabs and screened for M. agassizii and M. testudineum using a quantitative Polymerase Chain Reaction (qPCR) assay as described previously [32].

RNA extraction and sequencing for RNA-Seq
Total RNA was isolated from aliquots of whole blood with minor modifications to the total RNA isolation protocol. Briefly, whole blood samples were thawed on ice, homogenized in RNA lysis buffer, and aliquoted before extracting with acid phenol chloroform twice. Ethanol (100%) was added to each sample and passed through a glass-fiber filter, which binds RNA before eluting with nuclease-free water (mirVana miRNA Isolation Kit with phenol, Ambion, #AM1560, Carlsbad, California). The RNase-Free DNase Set (Qiagen, #79254, Valencia, CA) and RNeasy MinElute Cleanup Kit (Qiagen, #74204, Valencia, CA) were used to treat total RNA for residual DNA and salt contamination. Extracted total RNAs were sent to the Yale Center for Genomic Analysis (YCGA; West Haven, CT) to generate cDNA poly-A-enriched Illumina libraries that were run on two lanes of the Illumina HiSeq 2500 in High Output mode using 75-bp paired-end reads. To avoid batch effects, libraries for all individuals were split across the two sequencing lanes and then read files from the two lanes were concatenated for each sample.
To explore sample variance, we used Principal Components Analysis (PCA) of regularized log (rlog) count data generated by DESeq2 v1.20.0 [49] (Bioconductor v3.7), which transforms values onto a log 2 scale and normalizes for differences in sequencing depth. We assessed PCA results for variables of interest (experimental immune group, sex), as well as potentially confounding experimental variables: venipuncture site (subcarapacial vs. jugular), captive vs. wild, and sample collection year. In addition to experimental immune group and sex, PCA showed a result for venipuncture site as well. For this reason, we tested for differentially expressed genes (DEGs) for venipuncture site in addition to the DE analyses for experimental immune group and sex to remove it as a potentially confounding variable.
We identified DEGs with DESeq2 [49-51] after excluding four low coverage samples (see results section). Counts were normalized for library size internally, fitted to a negative binomial distribution, and corrected for multiple testing using the Benjamini-Hochberg method (FDR < 0.05). To control for variance associated with experimental group and sex, we used the multifactorial (two-factor) approach that employs the Wald test when evaluating differential expression. The NI experimental group and males were set as the reference levels for the immune and sex-based analyses, respectively. We did not have adequate sampling to analyze the venipuncture site in the multifactor analysis, so we separately ran a one-factor analysis to compare expression of jugular venipuncture versus subcarapacial venipuncture. This one-factor analysis also used the Wald test. Because venipuncture site is a confounding variable, we removed DEGs associated with venipuncture site from the final list of DEGs obtained from experimental group and sex. All results reported are excluding these venipuncture site-associated genes and are divided into genes that are only differentially expressed among experimental immune groups (unique immune genes), only differentially expressed among sexes (unique sex-biased genes), and those genes differentially expressed among both experimental immune and sex groups (sex-biased immune genes). Analyses and results were run in the R statistical programming environment (http://www.R-project.com). All heatmaps were generated using rlog transformed, mean-centered counts and genes were clustered by Manhattan distance using the Ward method.

Enrichment analysis
The unique gene lists for experimental immune group and for sex were individually ranked by adjusted p-values (padj <0.05) and tested for functional enrichment of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways in g:Profiler [52]. Functional profiling of differentially expressed genes were queried against the mouse database, excluding in silico curated terms, as an ordered list ranked by their adjusted p-values. Significant GO categories (padj < 0.05) were identified using Fisher's one-tailed test corrected for multiple testing using the Benjamini-Hochberg method (FDR<0.05) and hierarchically filtered by best-per-parent (moderate) parameters. We visualized significant GO terms (p < 0.05) using Reduce and Visualize Gene Ontology (REVIGO, [53]), which uses simRel scores as a measure of semantic similarity. A user-provided threshold value of 0.7 was selected for clustering.

Data generation and processing
All samples were sequenced across both Illumina HiSeq lanes to avoid batch effects. We obtained 153 Gb of sequence for the 25 individuals, averaging 30 million reads per individual (N = 25; 3-53 million reads/individual). Average sequence coverage was 2.6-fold lower for four individuals (two SI, two NI) than other individuals, which was evident in PCA (S1 Fig). Removing the four low-coverage samples resulted in a total of 21 samples with 25-53 million paired reads/individual. We only present results from this high-coverage dataset.
After read trimming, 552 million total reads were retained (11-43 million reads/individual) that ranged in length from 36-75 bp. Using STAR yielded a 91.7% mapping rate to the Gopherus agassizii 1.0 reference genome, of which 87.1% of reads were uniquely mapped, 4.1% of reads were multiply mapped, and 0.4% of reads were mapped to too many loci (only uniquely mapped reads were retained; S1 Table). Gene expression levels were quantified as the summation of unique fragments mapped to each exon. Overall, individuals in the male and Medium Infection (MI) groups exhibited the highest gene expression variance (Fig 1). However, six of the seven MI individuals were also male, so it is not possible to distinguish whether high variance is a trait specific to males or to having medium Mycoplasma spp. infection.

Data exploration and differential expression
Venipuncture site. Principal Components Analysis (PCA) showed no obvious patterns for two technical variables that were tested: captive vs. wild animals and collection year, suggesting these technical artifacts in our data are weak or absent (S1 Fig). We did observe a pattern associated with jugular vs subcarapacial venipuncture, which we attribute to the presence of varied amounts of lymph fluid in blood samples collected via subcarapacial venipuncture (S1 Fig). Blood collected via subcarapacial venipuncture can be collected within 1-2 min and without extensive animal handling, making it the preferred collection technique by managing agencies [36]. However, the subcarapacial plexus is proximal to lymphatic vessels, and may result in inadvertent aspiration of lymph when collecting blood [37]; jugular venipuncture does not result in any lymph admixture.
There were 53 total genes differentially expressed between venipuncture sites (S1 Appendix). Of these, 38 were unique to venipuncture analysis and showed some involvement in lymph-associated processes, such as regulation by host of viral transcription, T cell activation, eosinophil migration, and response to bacterium (S2 and S3 Figs; S2 Table). To attempt to mitigate the effects of lymph on our immune and sex-based analyses, we post hoc removed the 15 venipuncture-site DEGs that were also differentially expressed in the immune and/or sex comparisons. This was done to prevent interpretation of the lymph signal as a blood-based immune or sex-biased signal (all jugular venipuncture individuals were female). We note that there have not been RNA-Seq differential expression studies performed on desert tortoises, making this a novel and important result to account for in the design of future studies.

Immune and sex-biased results
After removing venipuncture site-associated genes, the multifactor Wald test yielded a total of 40 genes that were uniquely differentially expressed among experimental immune groups (Fig  2, S2A and S2B Appendix). Of these 40 genes, 14 were unique to the MI-NI comparison, 21 DEGs were unique to the SI-NI comparison, and five DEGs were differentially expressed in both comparisons (ABHD8, CDO1, RNF125, cell surface hyaluronidase-like, gopAga1_ 00017729). The 5 genes differentially expressed in both comparisons were upregulated relative to control (NI) with large log2-fold changes of 14.1-16.4 (SI-NI comparison values). In other words, these five genes are transcribed to differing degrees depending on whether a tortoise is infected with Mycoplasma spp. or not. There were 61 enriched GO terms for Biological Processes ( Table 2, S3 Table), but many of these contained only one or a few DEGs. The analysis produced zero enriched KEGG pathways. Semantic clustering of enriched GO terms using REVIGO yielded eight major clusters including cysteine catabolism, response to thyroid hormone, regulation of type 1 interferon production, and regulation of fibroblast apoptotic process ( Fig 3A).  There were 1,037 genes that were uniquely differentially expressed between females and males (Fig 4, S3 Appendix). Of these, 368 were upregulated in females relative to males with the greatest log2-fold change occurring in PACSIN3 (4.6), SH3GL3 (4.3), gopAga1_00019967 (2.9), gopAga1_00019277 (2.1), ZSWIM4 (2.1). There were 669 genes downregulated in females relative to males and the DEGs with the greatest log2-fold change were RGCC (-5.6), SMC2 (-3.0), CDK1 (-2.9), NEK2 (-2.8), gopAga1_00016417 (-2.7). The g:Profiler results for sex-biased DEGs produced 116 enriched Biological Processes and 79 enriched KEGG pathways (Tables 3 and 4, S4 Table, S4 Appendix). REVIGO produced 14 semantic clusters including iron ion transport, iron ion homeostasis, response to UV-C, regulation of interferon-beta production, cellular response to organonitrogen compound, and regulation of lysosomal protein catabolism (Fig 3B). Eight genes were differentially expressed both in the experimental immune and sex-biased analyses, which we refer to here as sex-biased immune genes (Fig 5). These sex-biased immune genes were identified as TMEM135, ST6GALNAC4, IGHE, TRIM3, and TRIM68-like, gopAga1_00007704, gopAga1_00017285, and gopAga1_00017523.

Discussion
Infections from pathogenic bacteria such as Mycoplasma spp. impact the morbidity and mortality of both wild and captive tortoise populations, which we use to understand more about the general immune system of non-avian reptiles. To do this we analyzed the transcriptional response of Mojave desert tortoises to infection by M. agassizii. We identified 40 uniquely differentially expressed genes among experimental immune groups (Fig 2), of which 5 genes were differentially expressed both in medium infection (MI) and severe infection (SI) groups relative to control (NI) animals. Given the strong health differences among these groups, it was surprising to discover that transcription in desert tortoises was foremost influenced by sex with 1,037 genes uniquely differentially expressed between males and females (Fig 4; S3 Appendix). We identified eight genes that were significantly differentially expressed in both analyses, which we consider to be sex-biased immune genes (Fig 5). Finally, results showed that the method most accepted to draw blood from tortoises (subcarapacial venipuncture) is not ideal for gene expression analysis due to varying amounts of lymph aspirate that biases the gene expression signal based on the amount of aspirate in the sample. Our results provide a systemic view of the effects of mycoplasmosis and sex-biased transcriptional differences in desert tortoises that will aid future management and conservation practices and shed new light on the immune response of non-avian reptiles broadly.

Infection-based differential expression
Tortoises, like many ectotherms, rely on broad non-specific innate immune responses such as non-specific leukocytes, lysozymes, antimicrobial peptides, the complement pathway and phagocytic B cells as primary lines of defense against pathogens [1,54,55]. Adaptive immune reactions mediated by T and B cells are induced in tortoises; however, their cell-mediated and humoral responses may be slow (weeks to years; [29,31]) or fail to develop into novel antigens [56], and do not consistently demonstrate evidence of memory response [1,57]. Indeed, we found many general immune-related responses to be enriched among experimental immune groups, including GO:0002275 (myeloid cell activation involved in immune response), GO:0009410 (response to xenobiotic stimulus), and GO:0002252 (immune effector process). Previous studies on turtles revealed expression of immune-related genes [2,31,38,58-63] associated with both innate and adaptive immune functions. We expected differential expression of cytokines (e.g. IFN-γ, interleukins, TNFα1) to be higher in tortoises with M. agassizii infection because they play important roles in modulating host defense responses to   :0044237 cellular metabolic process  433  TRIM25, PEMT, DAZAP2, C1D, SMARCB1, TCF25, BRPF1, DDX18,  CHORDC1, RAC1, RMND5A, GTF2F1, RELB, RAB8A, DVL3, RUVBL2 immediate and long-term pathogen exposure. While related categories GO:0002718 (regulation of cytokine production involved in immune response, Table 2) and regulation of type I interferon production ( Fig 3A) were enriched, we did not find these specific genes to be differentially expressed (Table 2). Moreover, western painted turtles (Chrysemys picta bellii; [58]) demonstrated a unique repertoire of toll-like receptors (TLRs) including TLR15-like receptor known in the response of birds to bacterial pathogens [64]. In our data we found related pattern recognition signaling categories such as GO:0039529 and GO:0039536 (RIG-I signaling pathway) as well as GO:1900745 (regulation of p38 MAPK cascade) and GO:1901224 (regulation of NIK/ NF-κB signaling), which are involved in innate and adaptive immune activation and maintenance. On the other hand, these pathways are important mediators of inflammation and if prolonged, may exacerbate the effects of M. agassizii infection and URTD. Chinese soft-shelled turtles (Pelodiscus sinensis) infected with pathogenic bacteria also demonstrated differential expression of innate immune genes including IL-8, serum amyloid A (SAA), CD9, CD59, activating transcription factor 4 (ATF4) and cathepsin L genes [62], pointing to an initial non-specific innate immune response, followed by later moderate adaptive responses to combat remaining pathogens. Given that this study was designed to assess chronic rather than acute immunological responses, it was not surprising that these five genes differentially expressed among the immune groups. However, we did observe that CD19 was downregulated in both the desert tortoise MI and SI groups relative to uninfected controls, although not significantly in the MI group. CD19 is an antigen expressed by both subsets of B lymphocytes. B1 cells are innate-like effectors that produce natural antibodies and exhibit phagocytic activity in fish, amphibians, and reptiles, including turtles/tortoises [1,65,66]. Additionally, B2 cells are responsible for generating antigen-specific antibodies against foreign pathogens. In Mojave desert tortoises, the mean infection intensity of M. agassizii is negatively correlated with the mean number of lymphocytes [65], which could provide protection via phagocytosis during early infection or long-term humoral immunity. Given that MI and SI individuals have been chronically infected with M. agassizii and tested positively for acquired antibodies, the latter is more likely in this study. Unlike the MI group, CD19 expression was significantly reduced in SI individuals, suggesting that the suppression of B lymphocytes may increase susceptibility to infection and morbidity. Taken together, CD19 may be a key B lymphocyte antigen in the chelonian immune response to infection, and may be a good gene target for future studies.
Of the 40 genes that were uniquely expressed among immune groups, 28 genes were previously annotated and are associated with a number of immune and metabolic processes. Genes ABHD8, CDO1, RNF125, cell surface hyaluronidase-like, and gopAga1_00017729 were significantly upregulated in MI and SI animals with high log2-fold change values (Fig 2). Most of the differentially expressed immune genes are involved in protein production, folding, and secretory domains, with additional broad functions related to both host defenses and mitigation of host-induced inflammatory responses. For example, cysteine dioxygenase type 1 (CDO1) is a   :05165 Human papillomavirus infection  29  DVL3, IKBKG, HDAC5, PPP2R5C, HDAC2, PSEN1, PPP2CA, RPS6KB1, PPP2R5E, SCRIB,  IFNAR1, HDAC3, TCF7L2, UBE3A, LAMC1, RBL1, PIK3R3, FADD, ATP6V1H, ITGA2B   5.41E-3 KEGG:04144 Endocytosis  23  AP2A2, RAB8A, KIF5B, GIT1, CYTH1, RAB5A, VPS26A, RUFY1, TFRC, USP8, CHMP5,  SH3GL3, WAS, VPS35, ARAP1, SNX2, PRKCI,  protein-coding gene that responds to glucagon, xenobiotic stimuli, bacterial insults, and organic substances, leading to enzymatic cysteine catabolic processes and dioxygenase activity ( Table 2; [67]). These processes aid in providing essential biomolecules such as vitamins, cofactors, antioxidants, and many defense compounds that frontline innate immune cells (e.g. macrophages and dendritic cells) need to protect cells or neutralize targeted antigens (e.g., invading bacteria). Relative to the control group, RNF125 was highly upregulated in both MI and SI groups with respective log2-fold changes of 14.9 and 16.5. As a ubiquitin ligase, RNF125 has been shown to reduce inflammatory signaling by targeting RIG-I [68] as well as regulating viral transcription in peripheral blood mononuclear cells [69]. RNF125 is primarily expressed in lymphoid tissues and is a positive regulator of T lymphocyte activation [70][71][72], which may be critical for modulating inflammatory pathways and adaptive host defense against infection in desert tortoises. Pathways with enriched DEGs included those involved in immune host defenses and regulatory processes. For example, cysteine catabolism plays an important role in many proteins and immune defenses such as mRNA transcription, regulation of macrophage chemotaxis, heme oxidation, etc. (Fig 3A; [73]). Other pathways such as regulation of type I interferon production (e.g. myeloid cell activation, basophil activation, immune effector processes) include a family of cytokines that are critically important in controlling host innate and adaptive immune responses to viral and bacterial infections, and other inflammatory responses (Fig 3A; [74]). There were also differential processes associated with regulation of thyroid hormones [75] and ionizing radiation [76] conditions. These enriched processes associate with modulating immune activities such as chemotaxis, phagocytosis, generation of reactive oxygen species (ROS), and cytokine synthesis at the cellular level based on hypo-and hyperthyroid conditions. Alternatively, they can interfere with the interactions of targeted cells such as dendritic cells and lymphocytes between innate and adaptive cell-mediated immunity, respectively [75].

Sex-biased differential expression
Sex is increasingly recognized as an important effector of the immune response [77][78][79]. Indeed, our results show striking differences in transcription between male and female tortoises evidenced by 1,037 sex-biased DEGs that yielded 116 enriched Biological Processes and 79 enriched KEGG pathways (Figs 3B and 4; Tables 3 and 4). Previous literature has shown that exposure to pathogens elicits sex-biased differential immunological responses. Females generally initiate stronger innate and adaptive immune responses, which can promote faster clearance of pathogens. This has been observed both in birds where females show increased T lymphocyte proliferation after parasite exposure compared to males [80] as well as in lizards where female-derived macrophages demonstrated greater phagocytic activity than malederived macrophages [81]. In humans, stimulation of TLR7 in plasmacytoid dendritic cells induces significantly greater expression of IFNα in females [82], and females show enhanced antibody responses with higher B cell numbers [83]. Consistent with the literature, our findings show that females exhibit upregulation of genes associated with cytokine production and host defense compared to males, including POLR3D, POLR3C, MAPKAPK2, MR1, ARID5A, SARS, SETD6, and PCBP2 (Table 3, S1 Appendix).
While females generally elicit stronger immune responses, this pattern is sometimes reversed. In this study, male tortoises displayed increased expression for biological processes involved in myeloid cell differentiation and activation, leukocyte mediated immunity, and positive regulation of tumor necrosis factor production (Figs 3B and 4; Tables 3 and 4). In mice, macrophages from males produced greater proinflammatory cytokines during the acute phase than their female counterparts [84]. Similarly, whole blood and neutrophils from human males produced greater levels of tumor necrosis factor than females [85,86]. Dysregulation of the immune response can lead to chronic infection or sepsis, which males are more prone to develop than females [87,88]. In contrast, females are more susceptible to inflammatory and autoimmune diseases as a result of stronger mounted immune responses [87]. Indeed, the TNIP1 gene is associated clinically with female-biased autoimmune diseases such as systemic lupus erythematosus and systemic sclerosis [88][89][90], which are diseases characteristic of an over-active immune system and here we also find TNIP1 to be significantly upregulated in female tortoises. While immune cell levels are not significantly different between male and female desert tortoises with no history of infection [12], our findings indicate that there is a sex-biased immune response following pathogen exposure. Whether immune cells vary between sex following infection is unknown but should be tested in future studies. Additionally, because immune cell function and levels change with seasonality and temperature in tortoises [12,106], it is possible that sex-biased gene expression is also affected by these factors. Overall, the sex-biased patterns in the literature are complex and sex-biased outcomes of Mycoplasma spp. infection and URTD warrants further study. Sex-biased immune responses depend on genetics and hormones. However, turtles (including the desert tortoise) lack sex chromosomes and instead have temperature-based sex determination [14]. How this fact affects immune gene expression, if at all, is unknown. Hormones, on the other hand, fluctuate seasonally and are responsive to changes in temperature in nonavian reptiles. Additionally, many sex hormone receptors directly regulate gene transcription by translocating to the nucleus and binding hormone response elements when activated. This suggests that hormones likely play a key role in sex-biased gene expression and immunity.
Indeed, female anole lizards suppress female-biased gene expression and exhibit higher levels of male biased genes when treated with testosterone [91]. Testosterone, which is higher in male than female vertebrates regardless of the mode of sex determination, also has immunosuppressive effects on immune cell activity [92]. Testosterone reduces the production of pro-inflammatory cytokines in mammals [93] and decreases cell-mediated immune reactions in birds, with greater immune suppression occurring in males than females [94]. In this study, the TMF1 gene, which is associated with the GO category for testosterone secretion, is significantly upregulated in male tortoises as expected. Whether testosterone causes immunosuppression in tortoises has not been investigated, but it is worth noting that testosterone levels in the desert tortoise vary seasonally [95]. This means that if there is an immunosuppressive effect of testosterone in tortoises then that immunosuppression may also vary seasonally. Testosterone levels begin to rise in April through July and peak in late August and September when male-male aggression, mating activity, and metabolism are greatest [96,97].
Additionally, testosterone is highly correlated with production of corticosterone [96], which also has immunosuppressive effects and has been shown to be higher on average in male tortoises relative to females [96,98]. Corticosteroids are generated as a stress response and we identified several genes upregulated in males that enriched for cellular response to stress (GO:0033554). These results together raise the question of whether immunosuppression in males during mating season predisposes them to Mycoplasma spp. infection due to a dampened immune response. Moreover, male tortoises also contact both sexes more frequently relative to females [30] and males travel greater distances, suggesting they may be more likely to spread M. agassizii (or other infections) among populations. Further investigation would be instructive for future wildlife management practices.
As expected, differentially expressed genes between male and female tortoises were enriched for metabolism, particularly those related to iron. Iron is an essential cofactor for many metabolic processes including cellular respiration, oxygen transport, and DNA synthesis. In this study, genes associated with iron homeostasis and iron ion import were expressed at higher levels in male tortoises (e.g., TFRC, SLC25A37). Transferrin receptor (TFRC) is required for the cellular uptake of iron through endocytosis and SLC25A37 is involved in transporting cytosolic iron into mitochondria via Mitoferrin-1. Sex differences in iron uptake may be associated with elevated energy demands in males because male tortoises have larger home range sizes and travel larger distances relative to females [99][100][101][102][103], and may also expend more acute energy requirements by engaging in combative activity over mates. Perhaps this is also why males exhibited higher gene expression variability than females (Fig 1). Additionally, female desert tortoises lay their eggs in April-mid July [104,105], which requires a large investment of energy and resources. Allocation of energetic resources during this period may deplete iron levels and could further contribute to transcriptional differences related to iron homeostasis between males and females. Cellular sequestration of iron is a recognized immunological response. Tortoises challenged with lipopolysaccharide demonstrated a reduction in plasma iron concentration [106]. Host import of iron prevents the pathogen from acquiring it, thereby limiting the rate of pathogen proliferation [107]. Given the diverse roles of iron in metabolism and immunity, future studies will be important to determine the functional effects of iron homeostasis. Fortunately, iron can be assayed easily and cheaply for wild and captive animals, lending itself as a good topic for future studies.

Conclusions and implications for future studies
We carried out RNA-Seq and differential expression analysis to identify the systemic host response of the Mojave desert tortoise to three levels of Mycoplasma agassizii infection. We identified 40 uniquely differentially expressed genes associated with infection. Among these were genes involved in protein production, secretory domains, host defenses, and mitigation of host-induced inflammatory responses. The genes ABHD8, CDO1, RNF125, cell surface hyaluronidase-like, and gopAga1_00017729 were upregulated in medium and severely infected animals relative to control, indicating these may be biomarkers of infection. A stronger result from these data is that sex plays a dominant role in determining gene expression, even among healthy and severely sick animals (1,037 sex-biased vs. 40 immune-based DEGs, respectively). We identified eight genes that were differentially expressed both by sex and infection status (i.e. sex-biased immune genes). Further assessment of gene expression before and during disease progression would be instructive. Because tortoises have temperaturebased sex determination, effects of sex hormones are of primary interest, including the potential immunosuppressive effect of testosterone and corticosterone during select seasonal periods of their activity. For future studies, subcarapacial venipuncture may result in aspirated lymph fluid, which will confound gene expression analysis, so jugular venipuncture would be a preferred method of blood collection for gene expression studies.
Supporting information S1