http://orcid.org/0000-0002-1458-4540
Figures
Abstract
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.
Citation: Xu C, Dolby GA, Drake KK, Esque TC, Kusumi K (2020) Immune and sex-biased gene expression in the threatened Mojave desert tortoise, Gopherus agassizii. PLoS ONE 15(8): e0238202. https://doi.org/10.1371/journal.pone.0238202
Editor: Y-h. Taguchi, Chuo University, JAPAN
Received: April 27, 2020; Accepted: August 11, 2020; Published: August 26, 2020
This is an open access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.
Data Availability: Sequence data is available via NCBI SRA accession numbers SRS3600670 to SRS3600682, BioSample accession numbers SAMN09727116 to SAMN09727140, and BioProject accession number PRJNA483175. All other results are available as appendices or bundled in Harvard Dataverse: https://doi.org/10.7910/DVN/XOKGGW.
Funding: This research was supported in part by the National Science Foundation EID grant #1216054 (K.K.D. and T.C.E.), U.S. Bureau of Land Management California Desert, California and Las Vegas, Nevada Districts agreement #L11PG00370 (K.K.D. and T.C.E), U.S. Geological Survey, Western Ecological Research Center [GX16ZC00BQAP2, GX16ZC00BQAP4] (G.A.D and K.K.), and Coyote Springs Investment LLC (K.K.D. and T.C.E). Any use of trade, product or firm names in this publication is for descriptive purposes only and does not imply endorsement by the U.S. government.
Competing interests: I have read the journal's policy and the authors of this manuscript have the following competing interests: G.A.D. is a member of the Board of Directors for the Desert Tortoise Council, a nonprofit conservation organization that protects desert tortoises and their habitats and funding support was provided by Coyote Springs Investment LLC. The DTC and Coyote Springs Investment LLC had no involvement in the design, implementation, or interpretation of this study. These do not alter our adherence to PLOS ONE policies on sharing data and materials.
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–7]), these resources are newer and less well-developed than those of model systems, such as human, chicken, and African clawed frog [8–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–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–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.
Materials and methods
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, 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 RNeasy® 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 Microtainer® 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.
Quality control and differential expression
We assessed reads for quality using FastQC v0.11.7 and MultiQC v1.5 [45] and performed read trimming using BBDuk v38.00 [46] with a Q28 Illumina quality score to remove low quality sequences; reads of at least 36 bp were retained. Using STAR v2.5.3a [47], the trimmed paired reads were mapped to the Gopherus agassizii 1.0 genome [2] and the gopAga1.0 annotation was converted to GTF format via gffreadv0.10.5 (https://daler.github.io/gffutils/api.html). Uniquely mapped reads were used to obtain gene-wise counts from exon sequences using featureCountsv1.6.1 [48]. Paired reads were counted as fragments using the “-p” flag. Multi-mapping and multi-overlapping reads were excluded by default [48].
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 log2 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 accessibility
Sequence data is available via NCBI SRA accession numbers SRS3600670 to SRS3600682, BioSample accession numbers SAMN09727116 to SAMN09727140, and BioProject accession number PRJNA483175 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA483175). All other results are available as appendices or bundled in Harvard Dataverse: https://doi.org/10.7910/DVN/XOKGGW.
Results
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.
Shapes are colored according to experimental groups (NI-No Infection, MI- Medium Infection, SI-Severe Infection). Circles represent female tortoises, squares represent male tortoises. Variables with differentially expressed genes include immune experimental group, sex, and venipuncture site (see S2 Fig).
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).
Genes are clustered by the Ward method according to Manhattan distance. All genes shown (rows) are statistically significantly expressed (adjusted α < 0.05). The tree on the left shows four clusters that reflect qualitatively different expression patterns, depicted in graphical schematics on the right-hand panel. These patterns can be further examined in future work. Color scale presents the amount of expression. Expression values are mean-centered, regularized log counts and colors are represented as z-score values. (^ ETS domain−containing protein Elk−1).
REVIGO treemaps showing semantically clustered enriched GO terms (colored tiles) with box size proportional to normalized adjusted p-value and clusters are shown for (A) experimental immune groups (top), and (B) sex-biased among males and females (bottom). The white labels on panel B are cell-cell recognition, cholangiocyte proliferation, and endothelial cell apoptotic process.
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).
Genes are clustered by the Ward method according to Manhattan distance. All genes shown (rows) are statistically significantly expressed (adjusted α < 0.05). Color scale presents the amount of expression. Expression values are mean-centered, regularized log counts and colors are represented as z-score values.
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.
Dotted line points to the three genes differentially expressed in all comparisons.
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 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 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–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–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 male-derived 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–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 non-avian 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–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 temperature-based 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 Table. Sequencing and mapping statistics for the 25 samples included in this study.
*Denotes samples removed from analysis due to low sequencing depth.
https://doi.org/10.1371/journal.pone.0238202.s001
(DOCX)
S2 Table. Enriched Gene Ontology (GO) terms for biological processes that are uniquely differentially expressed based on venipuncture site.
https://doi.org/10.1371/journal.pone.0238202.s002
(DOCX)
S3 Table. Enriched Gene Ontology (GO) terms for biological processes that are uniquely differentially expressed among experimental groups (No Infection, NI; Medium Infection, MI; Severe Infection, SI).
https://doi.org/10.1371/journal.pone.0238202.s003
(DOCX)
S4 Table. Enriched KEGG pathways for genes uniquely differentially expressed based on sex.
https://doi.org/10.1371/journal.pone.0238202.s004
(DOCX)
S1 Fig. Principal Components Analysis (PCA) of gene expression data based on different variables.
Color-coded by additional variables relevant to the experimental design: (A) low coverage samples, (B) venipuncture site, (C) wild vs. captive, (D) collection year. Only high coverage samples are shown for venipuncture site, wild vs. captive, and collection year. Venipuncture site showed a pattern and was analyzed through the DESeq2 pipeline (see S2 and S3 Figs, and S1 Appendix).
https://doi.org/10.1371/journal.pone.0238202.s005
(TIF)
S2 Fig. Heatmap of 38 genes that are uniquely differentially expressed by venipuncture site.
Subcarapacial blood draws may have aspirated lymph fluid due to proximal lymphatic vessels, however this is the preferred phlebotomy technique by management agencies. Expression values are mean-centered, regularized log counts and colors are represented as z-score values. Jugular, light blue; subcarapacial, tan.
https://doi.org/10.1371/journal.pone.0238202.s006
(TIF)
S3 Fig. REVIGO treemap for genes differentially expressed based on venipuncture site.
https://doi.org/10.1371/journal.pone.0238202.s007
(TIF)
S1 Appendix. Differentially Expressed Genes (DEGs) by venipuncture site (all, including non-unique genes).
https://doi.org/10.1371/journal.pone.0238202.s008
(XLSX)
S2 Appendix. Differentially Expressed Genes (DEGs) by experimental immune group (all, including non-unique genes).
https://doi.org/10.1371/journal.pone.0238202.s009
(XLSX)
S3 Appendix. Differentially Expressed Genes (DEGs) by sex (all, including non-unique genes).
https://doi.org/10.1371/journal.pone.0238202.s010
(XLSX)
S4 Appendix. Enriched Gene Ontology (GO) terms for biological processes that are uniquely differentially expressed based on sex.
https://doi.org/10.1371/journal.pone.0238202.s011
(XLSX)
Acknowledgments
We thank many people for their support with this research including Christina Aiello, Josephine Braun, Patrick Emblidge, Rachel Foster, Peter Hudson, Sydney Kelly, Kenneth Nussear, and Margarete Walden.
References
- 1. Zimmerman LM, Vogel LA, Bowden RM. Understanding the vertebrate immune system: insights from the reptilian perspective. Journal of Experimental Biology. 2010;213: 661–671. pmid:20154181
- 2. Tollis M, DeNardo DF, Cornelius JA, Dolby GA, Edwards T, Henen BT, et al. The Agassiz's desert tortoise genome provides a resource for the conservation of a threatened species. PLoS ONE. 2017;12: e0177708. pmid:28562605
- 3. Quesada V, Freitas-Rodríguez S, Miller J, Pérez-Silva JG, Jiang Z-F, Tapia W, et al. Giant tortoise genomes provide insights into longevity and age-related disease. Nat Ecol Evol. Springer US; 2018;3: 1–12. pmid:30510174
- 4. Georges A, Li Q, Lian J, O'Meally D, Deakin J, Wang Z, et al. High-coverage sequencing and annotated assembly of the genome of the Australian dragon lizard Pogona vitticeps. Gigascience. 2015;4: 45. pmid:26421146
- 5. Liu Y, Zhou Q, Wang Y, Luo L, Yang J, Yang L, et al. Gekko japonicus genome reveals evolution of adhesive toe pads and tail regeneration. Nature Communications. Nature Publishing Group; 2015;6: 1–11. pmid:26598231
- 6. Miller JK. Escaping senescence: demographic data from the three-toed box turtle (Terrapene carolina triunguis). Exp Gerontol. 2001;36: 829–832. pmid:11295516
- 7. Lind A, Lai YYY, Mostovoy Y, Holloway AK, Iannucci A, Mak AC, et al. A high-resolution, chromosome-assigned Komodo dragon genome reveals adaptations in the cardiovascular, muscular, and chemosensory systems of monitor lizards. 2019;: 1–48.
- 8. Hellsten U, Harland RM, Gilchrist MJ, Hendrix D, Jurka J, Kapitonov V, et al. The genome of the western clawed Frog Xenopus tropicalis. Science. American Association for the Advancement of Science; 2010;328: 633–636. pmid:20431018
- 9. Collins FS, Lander ES, Rogers J, Waterston RH, Conso IHGS. Finishing the euchromatic sequence of the human genome. Nature. Nature Publishing Group; 2004;431: 931–945. pmid:15496913
- 10. Zhang G, Li B, Li C, Gilbert MTP, Jarvis ED, Wang J. Comparative genomic data of the Avian Phylogenomics Project. Gigascience. 2014;3: 1063–8. pmid:25671091
- 11. Zapata AG, Varas A, Torroba M. Seasonal variations in the immune system of lower vertebrates. Immunol Today. 1992;13: 142–147. pmid:1580995
- 12. Sandmeier FC, Horn KR, Tracy CR. Temperature-independent, seasonal fluctuations in immune-function in a reptile, the Mojave desert tortoise (Gopherus agassizii). Canadian Journal of Zoology. 2016;94: 583–590.
- 13. Lotter H, Altfeld M. Sex differences in immunity. Seminars in Immunopathology; 2019;: 1–3. pmid:30242450
- 14. Bull JJ, Vogt RC. Temperature-dependent sex determination in turtles. Science. American Association for the Advancement of Science; 1979;206: 1186–1188.
- 15. Ellegren H, Hultin-Rosenberg L, Brunström B, Dencker L, Kultima K, Scholz B. Faced with inequality: chicken do not have a general dosage compensation of sex-linked genes. BMC Biol. 2007;5: 689–12. pmid:17883843
- 16. Grath S, Parsch J. Sex-biased gene expression. Annu Rev Genet. 2016;50: 29–44. pmid:27574843
- 17. Oertelt-Prigione S. The influence of sex and gender on the immune response. Autoimmunity Reviews. Elsevier B.V; 2012;11: A479–A485. pmid:22155201
- 18. Klein SL, Flanagan KL. Sex differences in immune responses. Nature Publishing Group. Nature Publishing Group; 2016;16: 626–638. pmid:27546235
- 19. Klein SL. The effects of hormones on sex differences in infection: from genes to behavior. Neuroscience adn Biobehavioral Reviews. 2000;: 627–368.
- 20. Jaillon S, Berthenet K, Garlanda C. Sexual dimorphism in innate immunity. Clinical Reviews in Allergy & Immunology; 2019;: 1–14. pmid:28963611
- 21. Ngo ST, Steyn FJ, McCombe PA. Gender differences in autoimmune disease. Frontiers in Neuroendocrinology. Elsevier Inc; 2014;35: 347–369. pmid:24793874
- 22. Meester I, Manilla-Muñoz E, León-Cachón RBR, Paniagua-Frausto GA, Carrión-Alvarez D, Ruiz-Rodríguez CO, et al. SeXY chromosomes and the immune system: reflections after a comparative study. Biology of Sex Differences; 2020;11: 1–13. pmid:31900228
- 23. Jacobson ER, Brown MB, Wendland LD, Brown DR, Klein PA, Christopher MM, et al. Mycoplasmosis and upper respiratory tract disease of tortoises: A review and update. The Veterinary Journal; 2014;201: 257–264. pmid:24951264
- 24. Smith R. Endangered and threatened wildlife and plants; determination of threatened status for the Mojave population of the desert tortoise. Federal Registrar. 1990;55: 12178–12191.
- 25. Saab JB, Losa D, Chanson M, Ruez R. Connexins in respiratory and gastrointestinal mucosal immunity. FEBS Letters. Federation of European Biochemical Societies; 2014;588: 1288–1296. pmid:24631537
- 26. Brown MB, Schumacher IM, Klein PA, Harris K, Correll T, Jacobson ER. Mycoplasma agassizii causes upper respiratory tract disease in the desert tortoise. Infect Immun. American Society for Microbiology (ASM); 1994;62: 4580–4586. pmid:7927724
- 27.
USFWS. Desert tortoise (Mojave population) recovery plan. Portland, Oregon, USA: U.S. Fish and Wildlife Service; 1994.
- 28.
USFWS. Revised recovery plan for the Mojave population of the desert tortoise (Gopherus agassizii). Sacramento, California, USA: U.S. Fish and Wildlife Service; 2011.
- 29. Aiello CM, Nussear KE, Esque TC, Emblidge PG, Sah P, Bansal S, et al. Host contact and shedding patterns clarify variation in pathogen exposure and transmission in threatened tortoise Gopherus agassizii: implications for disease modelling and management. Montgomery I, editor. J Anim Ecology. 2016;85: 829–842. pmid:26935482
- 30. Aiello CM, Esque TC, Nussear KE, Emblidge PG, Hudson PJ. The slow dynamics of mycoplasma infections in a tortoise host reveal heterogeneity pertinent to pathogen transmission and monitoring. Epidemiol Infect. Cambridge University Press; 2018;147: 1–10. pmid:30251621
- 31. Drake KK, Aiello CM, Bowen L, Lewison RL, Esque TC, Nussear KE, et al. Complex immune responses and molecular reactions to pathogens and disease in a desert reptile (Gopherus agassizii). Ecol Evol. 2019;9: 2516–2534. pmid:30891197
- 32. Braun J, Schrenzel M, Witte C, Gokool L, Burchell J, Rideout BA. Molecular methods to detect Mycoplasma spp. and testudinid herpesvirus 2 in desert tortoises (Gopherus agassizii) and implications for disease management. Journal of Wildlife Diseases. 2014;50: 757–766. pmid:25121400
- 33. Sandmeier FC, Weitzman CL, Maloney KN, Tracy CR, Nieto N, Teglas MB, et al. Comparison of current methods for the detection of chronic mycoplasmal URTD in wild populations of the Mojave desert tortoise (Gopherus agassizii). Journal of Wildlife Diseases. 2017;53: 91–101. pmid:27788056
- 34. Weitzman CL, Sandmeier FC, Tracy CR. Prevalence and diversity of the upper respiratory pathogen Mycoplasma agassizii in Mojave desert tortoises (Gopherus agassizii). Herpetologica. 2017;73: 113–9.
- 35. Allison LJ, McLuckie AM. Population trends in Mojave desert tortoises (Gopherus agassizii). Herpetological Conservation and Biology. 2018;13: 433–452.
- 36.
U.S. Fish and Wildlife Service. Health assessment procedures for the Mojave desert tortoise (Gopherus agassizii): A handbook pertinent to translocation. 2016 May pp. 1–75.
- 37. Hernandez-Divers S, Hernandez-Divers SJ, Wyneken J. Angiographic, anatomic and clinical technique descriptions of a subcarapacial venipuncture site for chelonians. Journal of Herpetological Medicine and Surgery. 2002;12: 32–37.
- 38. Drake KK, Bowen L, Lewison RL, Esque TC, Nussear KE, Braun J, et al. Coupling gene-based and classic veterinary diagnostics improves interpretation of health and immune function in the Agassiz’s desert tortoise (Gopherus agassizii). Conservation Physiology. 2017;5: 3429–17. pmid:28835840
- 39. Drake KK, Esque TC, Nussear KE, Defalco LA, Scoles-Sciulla SJ, Modlin AT, et al. Desert tortoise use of burned habitat in the Eastern Mojave desert. Jour Wild Mgmt. 2015;79: 618–629.
- 40. Désert C, Merlot E, Zerjal T, Bed'Hom B, Härtle S, Le Cam A, et al. Transcriptomes of whole blood and PBMC in chickens. Comparative Biochemistry and Physiology—Part D: Genomics and Proteomics. Elsevier Inc; 2016;20: 1–9. pmid:27442111
- 41. Waits DS, Simpson DY, Sparkman AM, Bronikowski AM, Schwartz TS. The utility of reptile blood transcriptomes in molecular ecology. Molecular Ecology Resources. 2019;20: 308–317. pmid:31660689
- 42. Stacy NI, Alleman AR, Sayler KA. Diagnostic hematology of reptiles. clinics in laboratory medicine. Elsevier Ltd; 2011;31: 87–108. pmid:21295724
- 43. Jacobson ER, Schumacher J, Green M. Field and clinical techniques for sampling and handling blood for hematologic and selected biochemical determinations in the desert tortoise, Xerobates agassizii. Copeia. American Society of Ichthyologists and Herpetologists (ASIH); 1992;1992: 237–241.
- 44. Wendland LD, Zacher LA, Klein PA, Brown DR, Demcovitz D, Littell R, et al. Improved enzyme-linked immunosorbent assay to reveal Mycoplasma agassizii exposure: a valuable tool in the management of environmentally sensitive tortoise populations. Clinical and Vaccine Immunology. 2007;14: 1190–1195. pmid:17626160
- 45. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32: 3047–3048. pmid:27312411
- 46. Bushnell B. BBMap: A fast, accurate, splice-aware aligner. United States; 2014. p. 2. pmid:22988817
- 47. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2012;29: 15–21. pmid:23104886
- 48. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general purpose program for assigning sequence reads to genomic features. Bioinformatics. 2014;30: 923–930. pmid:24227677
- 49. Anders S, McCarthy DJ, Chen Y, Okoniewski M, Smyth GK, Huber W, et al. Count-based differential expression analysis of RNA sequencing data using R and Bioconductor. Nat Protoc. 2013;8: 1765–1786. pmid:23975260
- 50. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biology. 2014;15: 31–21. pmid:25516281
- 51. Anders S, Huber W. Differential expression analysis for sequence count data. Genome Biology. 2nd ed. BioMed Central; 2010;11: R106–12. pmid:20979621
- 52. Reimand J, Arak T, Adler P, Kolberg L, Reisberg S, Peterson H, et al. g:Profiler-a web server for functional interpretation of gene lists (2016 update). Nucleic Acids Research. Oxford University Press; 2016;44: 83–89. pmid:27098042
- 53. Supek F, Bošnjak M, Škunca N, Šmuc T. REVIGO summarizes and visualizes long lists of gene ontology terms. Gibas C, editor. PLoS ONE. 2011;6: e21800–9. pmid:21789182
- 54.
Rios FM, Zimmerman LM. Immunology of reptiles. Chichester, UK: eLS, John Wiley & Sons, Ltd; 2015. https://doi.org/10.1002/9780470015902.a0026260
- 55. Zimmerman LM. Reptilia: Humoral immunity in reptiles. In: E C, editor. Advances in comparative immunology. 2018. pp. 751–772.
- 56. Sandmeier FC, Tracy CR, duPré S, Hunter K. A trade-off between natural and acquired antibody production in a reptile: implications for long-term resistance to disease. Biology Open. The Company of Biologists Ltd; 2012;1: 1078–1082. pmid:23213387
- 57. Grey HM. Phylogeny of the immune response. J Immunol. 1963;91: 819. pmid:14106306
- 58. Shaffer HB, Minx P, Warren DE, Shedlock AM, Thomson RC, Valenzuela N, et al. The western painted turtle genome, a model for the evolution of extreme physiological adaptations in a slowly evolving lineage. Genome Biology. BioMed Central Ltd; 2013;14: R28. pmid:23537068
- 59. Drake KK, Bowen L, Nussear KE, Esque TC, Berger AJ, Custer NA, et al. Negative impacts of invasive plants on conservation of sensitive desert wildlife. Ecosphere. 2016;7: e01531–20.
- 60. Bowen L, Miles AK, Drake KK, Waters SC, Esque TC, Nussear KE. Integrating gene transcription-based biomarkers to understand desert tortoise and ecosystem health. ecohealth. Springer US; 2015;: 1–12. pmid:25561383
- 61. Dolby GA, Morales M, Webster T, DeNardo D, Wilson M, Kusumi K. Discovery of a new TLR gene and gene expansion event through improved desert tortoise genome assembly with chromosome-scale scaffolds. Genome Biology and Evolution. 0AD;: 1–52.
- 62. Zhou X, Guo Q, Dai H. Identification of differentially expressed immune-relevant genes in Chinese soft-shelled turtle (Trionyx sinensis) infected with Aeromonas hydrophila. Veterinary Immunology and Immunopathology. 2008;125: 82–91. pmid:18572252
- 63. Xu J, Zhao J, Li Y, Zou Y, Lu B, Chen Y, et al. Evaluation of differentially expressed immune-related genes in intestine of Pelodiscus sinensis after intragastric challenge with lipopolysaccharide based on transcriptome analysis. Fish & Shellfish Immunology. 2016;56: 417–426.
- 64. Brownlie R, Allan B. Avian toll-like receptors. Cell Tissue Res. Springer-Verlag; 2010;343: 121–130. pmid:20809414
- 65. Sandmeier FC, Weitzman CL, Tracy CR. An ecoimmunological approach to disease in tortoises reveals the importance of lymphocytes. Ecosphere. 3rd ed. 2018;9: e02427–12.
- 66. Zimmermann NFA, Ritz CM, Hellwig FH. Further support for the phylogenetic relationships within Euphorbia L. (Euphorbiaceae) from nrITS and trnL–trnF IGS sequence data. Plant Syst Evol. 2010;286: 39–58.
- 67. Brait M, Ling S, Nagpal JK, Chang X, Park HL, Lee J, et al. Cysteine dioxygenase 1 is a tumor suppressor gene silenced by promoter methylation in multiple human cancers. Wong C-M, editor. PLoS ONE. 2012;7: e44951–19. pmid:23028699
- 68. Arimoto K-I, Takahashi H, Hishiki T, Konishi H, Fujita T, Shimotohno K. Negative regulation of the RIG-I signaling by the ubiquitin ligase RNF125. Proceedings of the National Academy of Sciences. National Academy of Sciences; 2007;104: 7500–7505. pmid:17460044
- 69. Shoji-Kawata S, Zhong Q, Kameoka M, Iwabu Y, Sapsutthipas S, Luftig RB, et al. The RING finger ubiquitin ligase RNF125/TRAC-1 down-modulates HIV-1 replication in primary human peripheral blood mononuclear cells. Virology. 2007;368: 191–204. pmid:17643463
- 70. Chu P, Pardo J, Zhao H, Li CC, Pali E, Shen MM, et al. Systematic identification of regulatory proteins critical for T-cell activation. J Biol. BioMed Central; 2003;2: 21–16. pmid:12974981
- 71. Giannini AL, Gao Y, Bijlmakers M-J. T-cell regulator RNF125/TRAC-1 belongs to a novel family of ubiquitin ligases with zinc fingers and a ubiquitin-binding domain. Biochemical Journal. 2008;410: 101–111. pmid:17990982
- 72. Zhao H, Li CC, Pardo J, Chu PC, Liao CX, Huang J, et al. A novel E3 ubiquitin ligase TRAC-1 positively regulates T cell activation. The Journal of Immunology. 2005;174: 5288–5297. pmid:15843525
- 73. Wu G, Fang Y-Z, Yang S, Lupton JR, Turner ND. Glutathione metabolism and its implications for health. J Nutr. 2004;134: 489–492. pmid:14988435
- 74. Arimoto K-I, Miyauchi S, Stoner SA, Fan J-B, Zhang D-E. Negative regulation of type I IFN signaling. Wallet M, Justement L, Bishop G, Rane M, Iragavarapu-Charyulu V, editors. J Leukoc Biol. John Wiley & Sons, Ltd; 2018;103: 1099–1116. pmid:29357192
- 75. De Vito P, Incerpi S, Pedersen JZ, Luly P, Davis FB, Davis PJ. Thyroid hormones as modulators of immune activities at the cellular level. Thyroid. Mary Ann Liebert, Inc., publishers; 2011;21: 879–890. pmid:21745103
- 76. Manda K, Glasow A, Paape D, Hildebrandt G. Effects of ionizing radiation on the immune system with special emphasis on the interaction of dendritic and T cells. 2012;2: 1–9.
- 77. Clayton JA. Applying the new SABV (sex as a biological variable) policy to research and clinical care. Physiology & Behavior. Elsevier; 2018;187: 2–5. pmid:28823546
- 78. Klein SL, Schiebinger L, Stefanick ML, Cahill L, Danska J, de Vries GJ, et al. Opinion: Sex inclusion in basic research drives discovery: Fig 1. Proceedings of the National Academy of Sciences. 2015;112: 5257–5258. pmid:25902532
- 79. Zakiniaeiz Y, Cosgrove KP, Potenza MN, Mazure CM. Balance of the sexes: Addressing sex differences in preclinical research. Yale J Biol Med. 2016;89: 255–259. pmid:27354851
- 80. Tschirren B, Fitze PS, Richner H. Sexual dimorphism in susceptibility to parasites and cell-mediated immunity in great tit nestlings (vol 72, pg 839, 2003). J Anim Ecology. John Wiley & Sons, Ltd (10.1111); 2004;73: 814–814.
- 81. Mondal S, Rai U. Sexual dimorphism in phagocytic activity of wall lizard's splenic macrophage and its control by sex steroids. General and Comparative Endocrinology. 1999;116: 291–298. pmid:10562459
- 82. Berghöfer B, Frommer T, Haley G, Fink L, Bein G, Hackstein H. TLR7 ligands induce higher IFN-α production in females. The Journal of Immunology. 2006;177: 2088–2096. pmid:16887967
- 83. Abdullah M, Chai P-S, Chong M-Y, Tohit ERM, Ramasamy R, Pei CP, et al. Gender effect on in vitro lymphocyte subset levels of healthy individuals. Cellular Immunology. Elsevier Inc; 2012;272: 214–219. pmid:22078320
- 84. Marriott I, Bost KL, Huet-Hudson YM. Sexual dimorphism in expression of receptors for bacterial lipopolysaccharides in murine macrophages: A possible mechanism for gender-based differences in endotoxic shock susceptibility. Journal of Reproductive Immunology. 2006;71: 12–27. pmid:16574244
- 85. Aomatsu M, Kato T, Kasahara E, Kitagawa S. Gender difference in tumor necrosis factor-α production in human neutrophils stimulated by lipopolysaccharide and interferon-γ. Biochemical and Biophysical Research Communications. Elsevier Inc; 2013;441: 220–225. pmid:24140406
- 86. Moxley G, Stern AG, Carlson P, Estrada E, Han JF, Benson LL. Premenopausal sexual dimorphism in lipopolysaccharide-stimulated production and secretion of tumor necrosis factor. J Rheumatol. 2004;31: 686–694. pmid:15088292
- 87. Rubtsova K, Marrack P, Rubtsov AV. Sexual dimorphism in autoimmunity. J Clin Invest. 2015;125: 2187–2193. pmid:25915581
- 88. Adrianto I, Wang S, Wiley GB, Lessard CJ, Kelly JA, Adler AJ, et al. Association of two independent functional risk haplotypes in TNIP1with systemic lupus erythematosus. Arthritis & Rheumatism. 2012;64: 3695–3705. pmid:22833143
- 89. Gateva V, Sandling JK, Hom G, Taylor KE, Chung SA, Sun X, et al. A large-scale replication study identifies TNIP1, PRDM1, JAZF1, UHRF1BP1 and IL10 as risk loci for systemic lupus erythematosus. Nat Genet. Nature Publishing Group; 2009;41: 1228–1233. pmid:19838195
- 90. Allanore Y, Saad M, Dieudé P, Avouac J, Distler JHW, Amouyel P, et al. Genome-Wide Scan Identifies TNIP1, PSORS1C1, and RHOB as Novel Risk Loci for Systemic Sclerosis. McCarthy MI, editor. PLoS Genet. 2011;7: e1002091–13. pmid:21750679
- 91. Cox RM, Cox CL, McGlothlin JW, Card DC, Andrew AL, Castoe TA. Hormonally mediated increases in sex-biased gene expression accompany the breakdown of between-sex genetic correlations in a sexually dimorphic lizard. Am Nat. 2017;189.
- 92. Rettew JA, Huet-Hudson YM, Marriott I. testosterone reduces macrophage expression in the mouse of toll-like receptor 4, a trigger for inflammation and innate immunity. Biology of Reproduction. 2008;78: 432–437. pmid:18003947
- 93. D'Agostino P, Milano S, Barbera C, Di Bella G, La Rosa M, Ferlazzo V, et al. Sex hormones modulate inflammatory mediators produced by macrophages. Ann NY Acad Sci. 1999;876: 426–429. pmid:10415638
- 94. Fargallo JA, Martínez-Padilla J, Toledano-Díaz A, Santiago-Moreno J, Dávila JA. Sex and testosterone effects on growth, immunity and melanin coloration of nestling Eurasian kestrels. J Anim Ecology. 2nd ed. 2007;76: 201–209. pmid:17184369
- 95. Rostal DC, Lance VA, Grumbles JS, Alberts AC. Seasonal reproductive cycle of the desert tortoise (Gopherus agassizii) in the eastern Mojave Desert. Herpetological Monographs. Herpetologists' League; 1994;8: 72–82.
- 96. Lance VA, Rostal DC. The annual reproductive cycle of the male and female desert tortoise: Physiology and endocrinology. Chelonian Conservation and Biology. 2002;4: 1–12.
- 97. Rostal DC, Lance VA, Grumbles JS, Alberts AC. Seasonal reproductive cycle of the desert tortoise (Gopherus agassizii) in the eastern Mojave Desert. Herpetological Monographs. 1994;8: 72–82.
- 98. Drake KK, Nussear KE, Esque TC, Barber AM, Vittum KM, Medica PA, et al. Does translocation influence physiological stress in the desert tortoise? Acevedo-Whitehouse K, Kristensen TN, editors. Anim Conserv. 2012;15: 560–570.
- 99. O'Connor MP, Zimmerman LC, Ruby DE, Bulova SJ, Spotila JR. Home range size and movements by desert tortoises, Gopherus agassizii, in the eastern Mojave Desert. herpetological monographs. Herpetologists' League; 1994;8: 60–71.
- 100. Duda JJ, Krzysik AJ, Freilich JE. Effects of drought on desert tortoise movement and activity. Journal of Wildlife Management. 1999;63: 1181–1192.
- 101. Freilich JE, Burnham KP, Collins CM, Garry CA. Factors affecting population assessments of desert tortoises. Conservation Biology. 3rd ed. John Wiley & Sons, Ltd (10.1111); 2000;14: 1479–1489.
- 102. Harless ML, Walde AD, Delaney DK, Pater LL, Hayes WK. Home Range, Spatial Overlap, and Burrow Use of the Desert Tortoise in the West Mojave Desert. Copeia. The American Society of Ichthyologists and Herpetologists; 2009;2009: 378–389.
- 103. Franks BR, Avery HW, Spotila JR. Home range and movement of desert tortoises Gopherus agassizii in the Mojave Desert of California, USA. Endang Species Res. 2011;13: 191–201.
- 104. Lajbner ZXK, Pnini R, Camus MF, Miller J, Dowling DK. Experimental evidence that thermal selection shapes mitochondrial genome evolution. Scientific Reports. Springer US; 2018;8: 1–12. pmid:29311619
- 105. Turner FB, Hayden P, Burge BL, Roberson JB. Egg production by the desert tortoise (Gopherus agassizii) in California. Herpetologica. Allen Press; 1986;42: 93–104.
- 106. Goessling JM, Guyer C, Mendonca MT. More than Fever: Thermoregulatory responses to immunological stimulation and consequences of thermoregulatory strategy on innate immunity in gopher tortoises (Gopherus polyphemus). Physiol Biochem Zool. 2017;90: 484–493. pmid:28437174
- 107. Kluger MJ, Rothenburg BA. Fever and reduced iron: their interaction as a host defense response to bacterial infection. Science. American Association for the Advancement of Science; 1979;203: 374–376. pmid:760197