Transcriptomic Profiling of Fibropapillomatosis in Green Sea Turtles (Chelonia mydas) From South Texas

Sea turtle fibropapillomatosis (FP) is a tumor promoting disease that is one of several threats globally to endangered sea turtle populations. The prevalence of FP is highest in green sea turtle (Chelonia mydas) populations, and historically has shown considerable temporal growth. FP tumors can significantly affect the ability of turtles to forage for food and avoid predation and can grow to debilitating sizes. In the current study, based in South Texas, we have applied transcriptome sequencing to FP tumors and healthy control tissue to study the gene expression profiles of FP. By identifying differentially expressed turtle genes in FP, and matching these genes to their closest human ortholog we draw on the wealth of human based knowledge, specifically human cancer, to identify new insights into the biology of sea turtle FP. We show that several genes aberrantly expressed in FP tumors have known tumor promoting biology in humans, including CTHRC1 and NLRC5, and provide support that disruption of the Wnt signaling pathway is a feature of FP. Further, we profiled the expression of current targets of immune checkpoint inhibitors from human oncology in FP tumors and identified potential candidates for future studies.

Introduction

Sea turtle fibropapillomatosis (FP) is a tumor promoting disease that occurs globally in all seven species of sea turtles. The species known for having the highest FP prevalence is the green sea turtle (Chelonia mydas) (1). FP was first reported over 80 years ago in captured green turtles from the Florida Keys, USA (2). Despite sustained and/or increased detection in Florida (3, 4), Hawaii (5–7), and other areas globally (8, 9), FP was not detected in Texas, USA until 2010 (10). Since then the prevalence of this disease in Texas has increased dramatically. From 2010 to 2015 the yearly prevalence of FP in Texas among stranded turtles was a maximum of 5% each year (11). After 2015, the prevalence of FP in Texas grew substantially to over 35% in green turtles in 2018 (11).

The most visibly prominent feature of FP is the appearance of external cutaneous tumors on affected turtles (Figure 1) (12–15). Ocular, oral and internal tumors are also a disease component, affecting the visceral organs including the lungs, kidneys, and also bone (16). Juvenile sea turtles are recognized as the most affected age category, but FP has also been reported in adult nesting olive ridley turtles (Lepidochelys olivacea) in Costa Rica (17) indicating that this disease is not limited to a life stage. Tumors may spontaneously regress, or increase in size and/or number to the point of debilitation (18). The primary impact of FP is the decrease in ability of sea turtles to forage for food, swim, and avoid predation, thus affecting the overall survival of affected turtles (16). Clinical pathology profiles of FP affected turtles have also identified systemic effects of the disease including anemia, particularly among the most severely affected turtles (19–22). As green turtles are recognized by the International Union for the Conservation of Nature and Natural Resources as threatened with extinction (23) health problems like FP that affect remaining populations are of major concern to conservation efforts.

FIGURE 1

Figure 1. Two fibropapillomatosis affected juvenile green turtles (Chelonia mydas) undergoing rehabilitation at Sea Turtle Inc., South Padre Island, Texas, USA. Turtle 1 (A,B) and turtle 2 (C,D) both display prominent external tumors on their ventral side (A,C) and also ocular tumors (B,D).

The primary cause of FP is unknown. Several lines of evidence point to a herpesvirus [chelonid alphaherpesvirus 5 (ChHV5)] as a potential causal agent of FP (4, 24–27), which is reminiscent of human cancers with known viral origins including Epstein-Barr virus infection in Burkitt's lymphoma (28), human herpesvirus-8 infection in Kaposi sarcoma (29), and human papillomavirus infection in cervical cancer (30). There is a strong suggestion that it is ChHV5 together with other environmental and anthropogenic pressures that has contributed to the growth of FP prevalence [reviewed in (1)].

For sea turtles debilitated by FP that are admitted to rehabilitation facilities for treatment, the current standard of care is preoperative screening for internal tumors, which carry poor prognosis, followed by surgical excision of external tumors (31). A recent report from (15) shows that postoperative regrowth is seen in 50% of treated turtles and that the rehabilitation survival rate of FP affected turtles is low (25%). Given the significant advancements in human precision oncology and the potential applications in wildlife medicine (32–34) it seems reasonable to expect increasing the understanding of the biological features of FP may result in additional therapeutic approaches to the treatment of this disease. In turn this may increase survival rates, reduce tumor reoccurrence and minimize rehabilitation time. Moreover, insight into pathogenesis of tumor development may inform hypotheses related to disease etiology and expression.

In 2018, Duffy et al. (35) reported the first application of human precision oncology techniques to FP. The study undertook a molecular characterization of FP tumors in two juvenile green turtles from Florida, performing RNA-sequencing of seven tumor samples and three matched skin control samples. This was the first effort to characterize the genetic disruptions in FP and was a significant advancement in the molecular understanding of this disease. More recent work from this group has explored the shared and distinct genomic and transcriptomic profiles of internal and external FP tumors (36).

Here, we continue the efforts to understand the genetic perturbations of FP tumors through transcriptome sequencing of healthy and tumor tissue from stranded green turtles undergoing rehabilitation at Sea Turtle Inc., on South Padre Island, Texas, USA. We aimed to study FP in a geographically distinct population of green turtles, hypothesizing that both similarities and differences will emerge between the transcriptomic profiles of FP tumors in Texas and Florida. By drawing upon the wealth of biological knowledge established in humans, this study aimed to further characterize FP to understand the biological pathways disrupted in this disease and to identify potential therapeutic avenues for future investigation.

Results

Patient Characteristics

During 2017–2018, 35 tissue samples were collected from nine green turtles undergoing rehabilitation at Sea Turtle Inc., after stranding at South Padre Island, Texas. As summarized in Table 1, 26 FP tumors and three healthy tissue biopsies were collected from three turtles with FP. Healthy tissue biopsies were also obtained from six turtles with no visible evidence of FP.

TABLE 1

Table 1. Characteristics of juvenile green sea turtles, stranded at South Padre Island, Texas.

The most common cause of stranding among the sampled sea turtles was cold stunning (4/9 turtles). The remaining turtles stranded due to trauma from injury (2/9), fishing hook ingestion/entanglement (2/9), and lethargy (1/9). All turtles were juveniles with straight carapace lengths (nuchal notch to caudal extent of the carapace) and weight ranging from 27.4 to 55.7 cm (μ = 35.6 cm) and 3.1–20.6 kg (μ = 7.5 kg), respectively. For the three turtles that had FP tumors at collection, using the NOAA tumor burden scale of 1–3 (where 1 is least affected and 3 is most affected), two turtles were classified as category 1 and one turtle was category 2 (37). All FP affected turtles underwent multiple tumor excisions by CO₂ laser, with turtle FP03 requiring two surgeries due to tumor load. Turtles FP01 underwent five surgeries and FP02 three surgeries each requiring multiple surgeries due to tumor regrowth. All turtles were rehabilitated successfully at Sea Turtle Inc., and released.

General features observed from histopathology profiling of tumors collected and sequenced in this study confirmed samples as FP and identified common concurrent features including mild to moderate lymphocytic infiltration in all tumors, ulceration in tumor FP01–02 and secondary bacterial infection in tumor FP02–05.

Differential Gene Expression Analysis of Healthy vs. Tumor Tissue

Our analysis of gene expression compared eight healthy control skin samples to 25 tumor samples across eight turtles total. Two samples, both tumors, from the originally sequenced 35 were excluded during quality control. Of the 21,432 transcripts annotated in the CheMyd 1.0 reference genome, 17,698 (82.6%) had sufficiently large enough counts in this sample set to be retained for statistical analysis. Using DESeq2, the design of this analysis controlled for two factors of unwanted variation calculated with RUVseq as well as the individual turtle (to account for differences in the number of tumor samples obtained from each FP turtle). As expected, principal components analysis of global gene expression profiles for each sample clustered distinctly healthy control skin samples separately from FP tumors (Figure 2), separating the two sample types in the first principal component which explained 24.62% of sample variation.

FIGURE 2

Figure 2. Principal components analysis of gene expression in healthy (green) and tumor (orange) samples, incorporating 2 factors of unwanted variation calculated using RUVSeq. PC 1 separates healthy from tumor samples and explains 24.62% of the variation in gene expression across samples.

For differential gene expression analysis, a log fold change threshold of 1 and an alpha of 0.05 were chosen. This analysis identified 283 transcripts that were upregulated, and 90 transcripts that were downregulated (Figure 3). The 10 most significantly upregulated, and 10 most significantly downregulated genes are shown in Figure 4 and described in Table 2. All gene expression differences are summarized in Supplementary Table 1.

FIGURE 3

Figure 3. Volcano plot of differentially expressed genes in FP tumor samples compared to healthy control skin. The x-axis plots the log₂ fold change from healthy control skin to tumor samples, while the y-axis plots the –log₁₀ adjusted P-value resulting from the differential expression analysis.

FIGURE 4

Figure 4. The top 10 upregulated (upper panels) and downregulated (lower panels) differentially expressed genes in FP tumors. The x-axis in each plot panel separates the box plots of tissue samples into the healthy control skin group and tumor group. The y-axis plots the log₂-transformed normalized read count for each sample.

TABLE 2

Table 2. Top 10 upregulated and downregulated differentially expressed genes.

Functional Enrichment and Pathway-Based Analysis of Differentially Expressed Genes

To extend beyond individual differentially expressed genes and to examine the broader biological processes involved in FP tissue, turtle genes were mapped to their closest human gene ortholog. This then enabled analyses that draw upon well-characterized biological networks from human-based data to infer understandings in sea turtle FP.

Among the 17,698 sea turtle genes that underwent differential expression analysis between healthy tissue and tumors, 91.2% were successfully matched to a human ortholog through BLAST-based analysis. For the 373 significantly differentially expressed genes, 327 (87.7%) were matched to 302 human orthologs.

GSEA using GO terms and KEGG pathways was conducted using clusterProfiler. For this analysis we relaxed the significance threshold and considered genes with an adjusted p < 0.1 as input on the basis that biologically relevant expression changes may be reflected in this expanded list. This increased the number of genes to 459 of which 407 (88.7%) were matched to 374 human orthologs. For each duplicate ortholog the most significant differentially expressed gene value was used in downstream analyses.

From the GO based analysis, 465 nominally significant (p < 0.05) terms were identified as enriched within the 374 genes (Supplementary Table 2). The top ranked terms were primarily immune system related. The top 20 most significant GO biological process ontology terms are shown in Figure 5. When the genes overlapping between these GO terms are considered as an enrichment map (Figure 6), two functional modules emerge. One module specifically involves biological processes related to responses to bacteria and interferon-gamma and the second module is related to immune cells and immune processes.

FIGURE 5

Figure 5. A dot plot of the top 20 most significant GO biological process ontology terms enriched among differentially expressed genes in FP. The x-axis shows the ratio of differentially expressed genes in each term relative to the total number of genes in that term. Dot size shows the number of differentially expressed genes in that term and color is the gene set enrichment analysis test statistic with Benjamini–Hochberg adjustment for multiple testing.

FIGURE 6

Figure 6. An enrichment map of genes overlapping the top 20 most significant GO biological process ontology terms enriched among differentially expressed genes in FP. The size of each node represents the number of genes differentially expressed in that ontology term, with thickness of edges between nodes representing the number of genes overlapping. Node color represents the gene set enrichment analysis test statistic with Benjamini–Hochberg adjustment for multiple testing.

From the KEGG pathway based analysis, 32 nominally significant pathways, shown in Figure 7 (Supplementary Table 3), were identified as over-represented within the differentially expressed genes. The most significantly enriched pathway was hsa04060, cytokine-cytokine receptor interaction, adjusted p = 8.72 × 10⁻⁷. Again, considering an enrichment map (Figure 8) between these pathways it is clear that immunological processes are driving the differential expression of genes between FP tumors and normal tissue.

FIGURE 7

Figure 7. A dot plot of nominally significant KEGG pathways enriched among differentially expressed genes in FP. The x-axis shows the ratio of differentially expressed genes in each term relative to the total number of genes in that term. Dot size shows the number of differentially expressed genes in that term and color is the gene set enrichment analysis test statistic with Benjamini–Hochberg adjustment for multiple testing.

FIGURE 8

Figure 8. An enrichment map of genes overlapping nominally significant KEGG pathways enriched among differentially expressed genes in FP. The size of each node represents the number of genes differentially expressed in that KEGG pathway, with thickness of edges between nodes representing the number of genes overlapping. Node color represents the gene set enrichment analysis test statistic with Benjamini–Hochberg adjustment for multiple testing.

Targeted Analysis of Immune Checkpoint Dysregulation in FP Tumors

Inhibition of immune checkpoint molecules, including CTLA4 and PD1, are current therapeutic targets in human precision oncology. Given their high therapeutic potential, we specifically assessed the gene expression of immune checkpoint molecule orthologs in this data including both established and emerging targets, as reviewed in (38, 39) and summarized in Supplementary Table 4. Of the 13 molecules targeted from the literature, 12 were present as orthologs in the sea turtle transcriptome based on protein sequence homology, with PD-L1, B7–H3, and CTLA-4 matching multiple sea turtle proteins. TIGIT was not detected. Three immune checkpoint molecules, BTLA, PD-L2, and LAG3, were significantly upregulated in FP tumors (Figure 9, Table 3). The most significant being PD-L2 (LOC102933512), with a log₂fold increase of 3.27, p_adj = 0.0001.

FIGURE 9

Figure 9. Significant upregulation of three immune checkpoint genes in FP. The x-axis in each plot panel separates the box plots of tissue samples into the healthy control skin group and tumor group. The y-axis plots the log₂-transformed normalized read count for each sample.

TABLE 3

Table 3. Differential expression of immune checkpoint molecules in FP.

Similarities and Differences Between Gene Expression Profiles of FP From Sea Turtles in Texas and Florida

Previously in the first reported transcriptomic profiling of FP, Duffy et al. sequenced 3 healthy control samples and 7 tumor samples from 2 green turtles in Florida (35), assembled a de novo transcriptome and conducted a differential gene expression analysis. Here we compare our findings where possible with these earlier findings. Duffy et al. reported the top 10 up-regulated and top 10 down-regulated transcripts in their analysis of which 12 were annotated to a matching human ortholog. Supplementary Table 5 compares the differential expression results for these 12 genes between the Duffy et al. (35) turtles from Florida and the turtles in this current study. Of the 12 genes, only four (FNDC1, S1PR3, NES, NTN3) were statistically replicated in this study, however all showed the same direction of expression change.

Recognizing that the results from the two sample cohorts are derived from different sequencing and analytical pipelines we can compare the trend in overall log₂fold changes in gene expression where possible. Duffy et al. report their top 300 up-regulated and top 300 down-regulated transcripts. Where a human gene annotation was assigned to a transcript, we have correspondingly matched that gene to our data, allowing for multiple matches where the same gene symbol is assigned to multiple transcripts. Figure 10A shows the direct comparison of log₂fold changes in gene expression between 449 overlapping observations from Duffy et al. and this study. Two transcripts from the Duffy study, both corresponding to the TBX20 gene were excluded as outliers as the TBX20 gene had a log₂fold change of 21.8 in our study, in comparison to 2.9 and 3.6 in the Duffy study. The overlap between studies was highly correlated (R = 0.76, p < 2.2 × 10⁻¹⁶), with a decrease in correlation to R = 0.68 with inclusion of TBX20 data points. This high correlation was increased when considering only the subset of transcripts (N = 98) that were significantly differentially expressed in both studies (R = 0.86, p < 2.2 × 10⁻¹⁶) as seen in Figure 10B.

FIGURE 10

Figure 10. Comparison of differentially expressed genes reported in Duffy et al. to the current study. (A) log₂fold correlations between 449 overlapping transcripts between the two studies. (B) log₂fold correlations for the subset of overlapping transcripts that were differentially expressed in both studies.

Detection of ChHV5 Transcript Expression in Both Tumor and Normal Tissue

We aligned the RNA sequence data from this study to the ChHV5 reference genome to identify sample level evidence of the presence of the ChHV5 virus. Of the 104 defined gene features in the ChHV5 genome, 83 were detected in this sample set. Total alignment counts across all 83 transcripts ranged from 0 to 483,216, with <10 counts in the majority of samples (22/33). Across the eight healthy tissue samples RNA molecules aligning to ChHV5 virus transcripts were detected in five, including all three healthy tissue samples from the FP affected turtles. In tumor samples alignment to virus transcripts was detected in 20/25 samples with two tumors from the same turtle (FP01) exhibiting high levels of RNA molecules aligning to ChHV5 transcripts. Five tumors, all from turtle FP03 had no evidence of ChHV5 transcript expression. The number of tumors in each turtle and healthy tissue in which ChHV5 transcripts were detected are indicated in Table 1. We report with interest that the two tumors from turtle FP01 with the apparent highest virus levels exhibited a smooth, raised surface with altered pigmentation in comparison to surrounding normal tissue, in comparison to the more typical verrucous character of the other tumors in this study.

Discussion

The prevalence of FP in Texas is growing among juvenile green turtles inhabiting nearshore habitat. As a wildlife health problem threatening an endangered species there is considerable need for further biological insight into this disease. Here, we have conducted transcriptome sequencing to identify differential gene expression profiles between healthy skin samples and FP tumors. Our analysis has identified 373 significantly differentially expressed sea turtle transcripts between healthy and tumor tissue. Through protein sequence homology analysis, we identified, where possible, the best matching human ortholog for each transcript. This facilitated exploration of both known biology at the gene level, but also allowed GSEA to identify aberrantly affected biological processes and pathways.

At the gene level, among the most significantly up-regulated and down-regulated genes resulting from this analysis are several genes currently implicated as important in human oncology. Expression of LOC102930327 which was found to be orthologous with human MR1 (major histocompatibility complex, class I-related) is upregulated in sea turtle FP tumors (log₂FoldChange = 2.92, p_adj = 1.22 × 10⁻²⁵). MR1 has a known role in the innate immune system as a selector and/or restrictor of mucosal-associated invariant T (MAIT) cells, which have a role in microbial defense (40–42). Through antigen presentation by MR1, MAIT lymphocytes are activated by microbially-derived intermediates of riboflavin synthesis (43). MAIT cells are the dominant lymphocyte population in human skin and also have a role in promoting tissue repair and wound healing (44). As FP tumors have been observed to have bacterial infections previously (5) as well as in the current study, we posit that the identified MR1 upregulation in FP tumors is an antimicrobial response and while not directly related to FP tumorgenicity may contribute to a permissive tumor microenvironment through inflammation upregulation. It is possible that as FP tumors develop and tissue becomes compromised, facilitating bacterial infection, the infection in turn contributes to continued tumor growth.

Expression of CTHRC1 (collagen triple helix repeat containing 1) is also upregulated in FP tumors (log₂FoldChange = 6.60, p_adj = 6.06 × 10⁻²²). In humans, CTHRC1 encodes an extracellular matrix protein that was originally identified from an upregulation in models of arterial injury (45). Further research has led to the recognition of CTHRC1 as an oncogenic protein with a role in the promotion of tumor growth, invasion and metastasis in multiple malignancies. CTHRC1 expression is increased in hepatocellular carcinoma (46, 47), melanoma (48) and breast cancer (49) and in several studies increased mRNA and protein expression in tumors was associated with poor disease prognosis and survival (46, 49–52). Further, increased CTHRC1 has been linked to the upregulation of matrix metalloproteinases (MMPs) which have established roles in tumor progression and metastasis in humans (53). In hepatoma cells (54), and colorectal cancer cells (50) increased CTHRC1 activated expression of MMP-9, which is known to contribute to extracellular matrix breakdown and tumor promotion (53). Increased expression of MMP9 was associated with poorer outcomes in nasopharyngeal carcinoma patients (55) and together with increased expression of CTHCR1 and MMP7, higher expression of MMP9 was associated with a lower survival rate after surgery in non-small cell lung cancer patients (56). This known connection between CTHCR1 and MMP7/MMP9 is reflected in the current study, with significant upregulation of CTHRC1 in FP tumors accompanied by increased expression of MMP7 (log₂FoldChange = 4.26, p_adj = 8.87 × 10⁻⁹) and MMP9 (log₂FoldChange = 4.24, p_adj = 4.05 × 10⁻⁸; Supplementary Table 1). Together this suggests that CTHRC1 may be a key oncogenic driver in FP.

The most significantly downregulated gene expressed in FP tumors was TBX4 (log₂FoldChange = −7.35, p_adj = 1.00 × 10⁻¹⁵). Encoding T-box transcription factor 4, TBX4, belongs to a highly evolutionarily conserved family of transcription factors responsible for limb development with fundamental roles during embryogenesis (57, 58). Given that Tbx4 expression in chick embryos determines hind limb development (59, 60), and that the healthy control tissue samples used in the current project were obtained from punch biopsies of the hind flippers, it is possible that the higher expression of TBX4 in healthy tissue in this study is due to tissue source. This finding warrants further investigation through gene expression comparisons of hind flipper tissue with non-hind flipper healthy skin. Indeed a limitation of this study is that control tissue was only sourced from hind flipper biopsies as collection from other external healthy epidermal tissue was not pursued to minimize impact to study turtles. The sustained expression of a gene known as an embryonic transcription factor in adult tissue initially seems counterintuitive. However, TBX4 has been implicated as a transcription factor of importance in human adult lung fibroblasts (61). Furthermore, TBX4 has been shown to have decreased expression in human lung cancer associated fibroblasts in comparison to matched normal lung fibroblasts (61), providing preliminary support that TBX4 expression may be involved in FP biology.

Through alignment of RNA sequence data generated in this study to the ChHV5 genome we also were able to identify evidence of the presence of the ChHV5 virus in 20/25 tumor samples and 5/8 healthy tissue samples. The majority of samples did not have large numbers of reads aligning to the ChHV5 genome, which may be a product of lowly expressed ChHV5 transcripts, or transcriptionally inactive virus, outcompeted during transcriptome sequencing. A limitation of this study was that a complementary ChHV5 qPCR-based assay was not conducted on the samples studied to robustly identify the presence or absence of ChHV5 virus DNA. Therefore, the data generated in this study are insufficient to conclude that samples without ChHV5 transcripts are free of ChHV5 infection. Given the consistent detection of ChHV5 virus in FP tumors (25, 27, 62) we expect that further qPCR-based analysis of these samples would identify the presence of this virus in each tumor sample.

In comparison to previous work conducted by Duffy et al. (35) we see a strong correlation with the gene expression differences identified in the current study, despite geographical and methodological differences. Duffy et al. previously implicated the Wnt signaling pathway as a potential therapeutic avenue for FP and our current findings reiterate this potential. While we did not replicate the finding of Duffy et al. of an upregulation of WNT5A in FP tumors, among the top differentially expressed genes in this study several are involved in Wnt signaling pathways. This includes upregulated genes CTHRC1, NLRC5 and downregulated LGR5. CTHRC1 has been shown to selectively activate the WNT/planar cell polarity pathway (63) and that this may have a role in cervical cancer (64). NLRC5 has been identified as a regulator of the Wnt/β-catenin signaling pathway in clear cell renal carcinoma (65), and in the current study we observe increased expression of NLRC5 in FP (log₂FoldChange = 3.39, p_adj = 3.71 × 10⁻¹²). Upregulation of NLRC5 in FP may represent a physiological anti-tumor response in sea turtles. Studies have shown that NLRC5 is a transcriptional regulator of MHC class I genes (66) and through this, activates cytotoxic CD8⁺ T cells as part of an anti-tumor immune response (67). The increased tumor expression of NLRC5 is associated with better outcomes in human cancer (68) and is a target of interest in anti-tumor immunity in transmissible tumors of another endangered species; the Tasmanian devil (69). Furthermore, downregulation of LGR5 in FP (log₂FoldChange = −4.56, p_adj = 1.03 × 10⁻⁹) may also demonstrate evidence of an anti-tumor response through regulation of the Wnt/β-catenin signaling pathway, as increased expression of LGR5 is associated with poor prognosis in glioma (70), breast cancer (71) and oral squamous cell carcinoma (72). The consistent identification of the alteration of genes involved in Wnt signaling in FP indicates that further in-depth exploration of this pathway in FP is warranted.

Moving beyond individual genes, we next took a systems biology based approach drawing upon the human centered knowledge bases of Gene Ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. We annotated all differentially expressed genes where possible with their closest orthologous human protein match to enable a comprehensive analysis within the human centered knowledge bases. The GO based enrichment analysis identified 465 nominally significant ontology terms within the gene set. These were primarily immune system based, and among the top 20 terms when considering the genes overlapping between them two functional modules emerged. One module involved the most significantly enriched GO biological process (BP), GO:0009617, response to bacterium (p_adj = 1.77 × 10⁻¹³), while the other module was comprised of immune related processes. Enrichment of genes involved in biological responses to bacterial infection is intuitive given that FP tumors, as an externally exposed and compromised tissue have potential to become infected with environmental bacteria. This bacterial infection interplays with tumor-based processes meaning that the enrichment of immune related GO terms is likely driven by tumor biology but also contributed to by bacterial infection. This is also reflected in the KEGG pathway-based analysis where the GSEA identified that the most significantly enriched pathway was hsa04060, cytokine-cytokine receptor interaction (p_adj = 8.72 × 10⁻⁷). It is well-established in humans that cytokines have a fundamental role in the tumor microenvironment as well as inflammation and infection control (73). Here in FP, it is likely that the contribution of cytokine pathways is again an interplay between immune processes related to bacterial infection, and those related to tumorgenicity. This is supported by histopathological examination of the tumors in this study showing various degrees of lymphocytic infiltration. Given these findings, a fundamental question arises: why is this immune response insufficient to result in tumor regression in individuals that exhibit progressive disease?

In relation to this question, given that the differentially expressed genes at both the gene level and the systems level implicate significant involvement of the immune system in FP tumors, we explored immune checkpoint molecules in this data as a potential future therapeutic avenue for FP treatment. The inhibition of immune checkpoint molecules as a therapeutic approach for the treatment of cancer is a rapidly growing field (74). We targeted a non-exhaustive set of 13 established and emerging immune checkpoint molecules identified from the literature for assessment in the current dataset. Gene expression of three checkpoint molecules, BTLA, PD-L2, and LAG3 was significantly upregulated in FP tumors. Each of these molecules represents a potential drug target opportunity to explore in the treatment of FP. Trialing an existing human based inhibitor that is injectable directly into FP tumors may substantially improve the burden of FP management in sea turtle rehabilitation but would require significant investment and carefully designed clinical trials to identify a chemotherapeutic regimen. Given the recent increases in stranded sea turtles with FP requiring care in some regions, a therapy that facilitates rapid rehabilitation and return to the environment could be very beneficial.

Here we have described the transcriptomic profiling of FP tumors in green turtles, building upon prior work from (35) but also identifying novel insights into this disease. While in comparison to human tumors much remains unknown about the basic biology of sea turtle FP there is significant opportunity to further apply alternative approaches previously established in humans and model organisms to increase our understanding of this disease. Research approaches involving single cell sequencing, long read sequencing, metagenomics and proteomics all stand to substantially benefit the FP knowledge base. Investment in this area may even offer returns to human oncology (75). Wildlife health problems, such as FP, are increasing and will require a multidisciplinary approach to their mitigation.

Methods

Sample Collection

Sampling of green turtles was carried out under permit number TE181762-4 from the U.S. Fish and Wildlife Service and with the ethical approval of the University of Texas Rio Grande Valley's Institutional Animal Care and Use Committee (IACUC). The nine turtles involved in this study were undergoing rehabilitative care at Sea Turtle Inc., South Padre Island. At the time of collection, three turtles had external cutaneous FP tumors and for each turtle tumor burden was quantified and scaled according to the National Oceanic and Atmospheric Administration (NOAA) Ordinal Scale for Grading Fibropapillomatosis by Photographic Comparison as described in (37). From these three turtles, external FP tumors were sampled after removal during rehabilitative surgery using CO₂ laser resection. From all turtles a single healthy control skin sample was collected using 3–6 mm punch biopsies of the dorsal proximal hind flipper. Tissue samples were collected into RNA-later (Qiagen) and stored at −80°C until processing. All turtles sampled in this study were juvenile green turtles of unknown sex. Table 1 summarizes the individual characteristics of turtles in this study.

Histopathological Studies of Tissue Samples

For each tumor sample collected, ~50% of the sample was fixed in a 10% formalin solution. Tissues were processed into paraffin blocks, sectioned by 5 μm, mounted onto glass slides and stained with hematoxylin and eosin using routine methods. Diagnosis of fibropapilloma was confirmed based on previously described histopathological features (76).

RNA Extraction and RNA-Seq

Tumor and healthy tissue samples were homogenized using a rotor-stator homogenizer and total RNA was extracted using the Qiagen AllPrep DNA/RNA Mini kit, according to the manufacturer's instructions. RNA quality was assessed on a 2200 TapeStation System (Agilent Technologies) with sample RIN values ranging from 6.4 to 9.6 (μ = 8.4). mRNA sequencing libraries were generated from up to 250 ng of total RNA and were prepared using a KAPA RNA HyperPrep kit (Roche Diagnostics) with polyA selection and indexed using a KAPA Dual-Indexed adapter kit (Roche Diagnostics), according to manufacturer's instructions. Prepared libraries were evaluated via the 2200 TapeStation System and were pooled for paired-end 100 bp sequencing on a HiSeq 2500 system (Illumina).

Quality Control of RNA Sequence Data

Sequence data was demultiplexed using bcl2fastq. Raw sequencing reads were processed with Trim Galore (version 0.6.5, https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/), using Cutadapt (77) to remove adapter sequences and indexes from reads and to exclude low quality sequences using a Phred score of 30 and discarding reads with lengths shorter than 25 bp, unpaired reads were not retained. FastQC (version 0.11.9, https://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to assess the data quality of the processed reads. To avoid potential bias in analysis, two individual samples were excluded due to overall differences in total sequencing content, one tumor tissue (FP03-05, Beach Bum) had over 83 million more paired-end reads than the next highest sample and one healthy tissue (FP08-01, Robocop) had over 2.5 million less paired-end reads than the next lowest sample and had an identified technical collection issue post sequencing, whereby a less than optimal amount of tissue had been obtained in the healthy tissue biopsy collection. In the remaining 33 samples, the mean and median number of paired-end sequencing reads was 26,607,522 and 23,683,506 respectively, ranging from 4,447,485 to 70,675,924 reads.

Sequence Alignment and Quantification of Turtle Gene Expression

The CheMyd_1.0 reference assembly for C. mydas was obtained from NCBI [GenBank assembly accession number: GCA_000344595.1] (78). Reads were aligned to the reference assembly using HISAT2 (version 2.2.0) (79) with an overall alignment per sample ranging from 75.26 to 88.29% (μ = 84.19%). Transcript abundance was quantified in each sample using htseq-count (version 0.11.2) (80) at the gene level according to the defined genomic features of the CheMyd_1.0 assembly.

Differential Expression Analysis

Differential expression analysis was conducted with DESeq2 (version 1.28.1) (81). Genes with low counts were filtered from the analysis, retaining genes with more than five counts in two or more samples. This corresponded to 17,698 genes (83%). To calculate factors of unwanted variation in this sample using RUVSeq (version 1.22.0) (82), a first pass differential expression analysis was conducted to identify a set of in-silico empirical negative control genes. This corresponded to a model comparing 25 tumor samples to 8 healthy tissue samples with inclusion of a turtle covariate in the model to account for multiple tumor samples obtained from individual sea turtles. Using a test statistic p-value of >0.5 to identify genes not differentially expressed in the comparisons there were 5,730 genes identified as empirical control genes from the analysis. Using RUVSeq, two factors of unwanted variation were calculated from the data. These factors were then incorporated into the DESeq2 analysis design and an analysis was conducted to identify differentially expressed genes with a log₂ fold change threshold of 1 and a false discovery rate threshold of 0.05. Relative log expression plots of the data unadjusted and after normalization, including the two factors of unwanted variations from RUVSeq are shown in Supplementary Figures 1, 2, respectively. The differential expression analysis was visualized as a volcano plot using the “EnhancedVolcano” package (version 1.6.0) (83). Correlations of previously published differential gene expression changes with the results from this study were calculated using a Pearson correlation test in R (version 4.0.0).

Human Ortholog Identification

To enable pathway based and gene set enrichment analyses in this data sea turtle genes were matched to their closest matching human gene orthologs based on protein amino acid sequence homology. To accurately draw upon human biological understandings for inference in sea turtle fibropapilloma the defined set of 28,672 proteins for the CheMyd_1.0 reference was compared to the human landmark protein database using NCBI's protein BLAST (84–86) to identify the closest human orthologous match. Initially the blastp query was run with the following parameters: –max_target_seqs 1 –max_hsps 1 –evalue 1e−6. This identified a match for 27,961 sea turtle proteins (97.5%). For the remaining 711 sea turtle proteins, the blastp analysis was repeated reducing the threshold of the expect value (evalue) to 0.01. This matched a further 165 human proteins to turtle proteins. Finally, the evalue threshold was reduced to 1, matching a further 182 proteins. In total 28,308 protein sequences (98.7%) from the CheMyd_1.0 reference were successfully matched with a human landmark database protein.

Turtle and human protein accession numbers were converted to gene symbols using bioDBnet (87). A subset of 211 turtle protein isoforms of 101 individual turtle genes matched to multiple human genes. When there was concordance between the human and turtle gene symbols for one or more of the turtle protein isoforms, that gene symbol was used as the match between turtle and human. For the remaining genes, the gene symbol corresponding to the blastp result with the lowest evalue was used. For evalue ties the gene symbol of the result corresponding to the longest turtle protein isoform was used.

For the remaining 364 turtle proteins that could not be matched to a human protein, the sea turtle genome gene annotation was used. For 2,233 turtle genes, including predicted tRNA molecules and other non-coding genes, a corresponding protein sequence was not available and for these genes the sea turtle genome gene annotation was retained. Supplementary Table 6 provides a mapping from CheMyd_1.0 gff_file gene identifiers through to the human based gene matches used in this project.

Gene Set Enrichment Analysis

Gene set enrichment analysis (GSEA) was conducted using clusterProfiler (version 3.16.1) (88). Human gene symbols from the orthologous matching analysis of sea turtle genes were annotated with Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. An adjusted p-value of < 0.1 was used to identify input genes for this analysis. For each duplicate ortholog the most significant differentially expressed gene value was used in the GSEA analyses. The background gene set used was all orthologous human genes identified as matching a sea turtle gene. Test statistics were adjusted for multiple testing using the Benjamini–Hochberg method.

Detection of Chelonid Alphaherpesvirus 5 Sequences

The chelonid alphaherpesvirus 5 (ChHV5) reference assembly was obtained from NCBI [GenBank assembly accession number: GCA_008792165.1] (89). Reads were aligned to the reference assembly using HISAT2 (version 2.2.0) (79) and transcript abundance was quantified in each sample using htseq-count (version 0.11.2) (80) at the gene level according to the defined genomic features of the ChHV5 assembly. Differential gene expression was not conducted for ChHV5 with results only reported as a presence or absence of viral transcripts in each sample.

Data Availability Statement

The RNA-Seq data presented in this study, including raw fastq reads, are deposited in NCBI under BioProject ID: PRJNA672698 (https://www.ncbi.nlm.nih.gov/bioproject/PRJNA672698) and are publicly available.

Ethics Statement

The animal study was reviewed and approved by the U.S. Fish and Wildlife Service under permit number TE181762-4 and with the ethical approval of the University of Texas Rio Grande Valley's Institutional Animal Care and Use Committee (IACUC).

Author Contributions

NB, AL, NN, MD, JP, IM, MK, JG, JB, and JC contributed to conception and design of the study. NB, NN, Td, and MD coordinated and conducted sample collection. AL, ML, JC, and BS coordinated and conducted laboratory analysis of samples, with AL and ML generating sequencing data. NB, JP, and IM conducted data analysis and interpretation with NB performing statistical analyses. NB wrote the first draft of the manuscript. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was funded by the University of Texas Rio Grande Valley Transforming Our World Strategic Plan.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

This work was conducted in part in facilities constructed under the support of NIH grant C06 RR020547. The co-authors wish to acknowledge the ongoing and generous contributions from supporters of Sea Turtle Inc., a non-profit sea turtle conservation, rehabilitation, and education facility located on South Padre Island, Texas. Without public support for the Sea Turtle Inc., facilities, employees, and associates the findings presented here would ultimately not have been possible. A version of this article (90) was published on the BiorXiv preprint server.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2021.630988/full#supplementary-material

Supplementary Figure 1. Relative log expression plot of the gene expression data per sample prior to adjustment for two factors of unwanted variation calculated using RUVSeq.

Supplementary Figure 2. Relative log expression plot of the gene expression data per sample after adjustment for two factors of unwanted variation calculated using RUVSeq.

Supplementary Table 1. Results of differential gene expression analysis of healthy vs. tumor tissue.

Supplementary Table 2. Results of GO based analysis for functional enrichment of pathways in differentially expressed genes.

Supplementary Table 3. Results of KEGG based analysis for functional enrichment of pathways in differentially expressed genes.

Supplementary Table 4. Immune checkpoint molecules assessed.

Supplementary Table 5. Comparison to previously identified top differentially expressed genes in FP from Duffy et al. (2018) Communications Biology (main article reference 35).

Supplementary Table 6. Gene ID mapping from CheMyd_1.0 gene identifiers to human gene matches.

References

1. Jones K, Ariel E, Burgess G, Read M. A review of fibropapillomatosis in Green turtles (Chelonia mydas). Vet J. (2016) 212:48–57. doi: 10.1016/j.tvjl.2015.10.041

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Smith GM, Coates CW. Fibro-epithelial growths of the skin in large marine turtles Chelonia mydas (Linnaeus). Zoologica. (1938) 23:93–8.

3. Foley AM, Schroeder BA, Redlow AE, Fick-Child KJ, Teas WG. Fibropapillomatosis in stranded green turtles (Chelonia mydas) from the eastern United States (1980-98): trends and associations with environmental factors. J Wildl Dis. (2005) 41:29–41. doi: 10.7589/0090-3558-41.1.29

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Lackovich JK, Brown DR, Homer BL, Garber RL, Mader DR, Moretti RH, et al. Association of herpesvirus with fibropapillomatosis of the green turtle Chelonia mydas and the loggerhead turtle Caretta caretta in Florida. Dis Aquat Organ. (1999) 37:89–97. doi: 10.3354/dao037089

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Aguirre AA, Balazs GH, Zimmerman B, Spraker TR. Evaluation of Hawaiian green turtles (Chelonia mydas) for potential pathogens associated with fibropapillomas. J Wildl Dis. (1994) 30:8–15. doi: 10.7589/0090-3558-30.1.8

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Chaloupka M, Balazs G. Modelling the effect of fibropapilloma disease on the somatic growth dynamics of Hawaiian green sea turtles. Mar Biol. (2005) 147:1251–60. doi: 10.1007/s00227-005-0026-1

CrossRef Full Text | Google Scholar

7. Chaloupka M, Balazs GH, Work TM. Rise and fall over 26 years of a marine epizootic in Hawaiian green sea turtles. J Wildl Dis. (2009) 45:1138–42. doi: 10.7589/0090-3558-45.4.1138

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Adnyana W, Ladds PW, Blair D. Observations of fibropapillomatosis in green turtles (Chelonia mydas) in Indonesia. Aust Vet J. (1997) 75:736–42. doi: 10.1111/j.1751-0813.1997.tb12258.x

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Williams EH, Bunkley-Williams L, Peters EC, Pinto-Rodriguez B, Matos-Morales R, Mignucci-Giannoni AA, et al. An epizootic of cutaneous fibropapillomas in green turtles Chelonia mydas of the caribbean: part of a panzootic? J Aquat Anim Health. (1994) 6:70–8. doi: 10.1577/1548-8667(1994)006<0070:AEOCFI>2.3.CO;2

CrossRef Full Text | Google Scholar

10. Tristan T, Shaver DJ, Kimbro J, deMaar T, Metz T, George J, et al. Identification of fibropapillomatosis in green sea turtles (Chelonia mydas) on the Texas Coast. J Herpetol Med Surg. (2010) 20:109. doi: 10.5818/1529-9651-20.4.109

CrossRef Full Text | Google Scholar

11. Shaver DJ, Walker JS, Backof TF. Fibropapillomatosis prevalence and distribution in green turtles Chelonia mydas in Texas (USA). Dis Aquat Organ. (2019) 136:175–82. doi: 10.3354/dao03403

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Brooks DE, Ginn PE, Miller TR, Bramson L, Jacobson ER. Ocular fibropapillomas of green turtles (Chelonia mydas). Vet Pathol. (1994) 31:335–9. doi: 10.1177/030098589403100306

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Flint M, Limpus CJ, Patterson-Kane JC, Murray PJ, Mills PC. Corneal fibropapillomatosis in green sea turtles (Chelonia mydas) in Australia. J Comp Pathol. (2010) 142:341–6. doi: 10.1016/j.jcpa.2009.10.012

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Norton TM, Jacobson ER, Sundberg JP. Cutaneous fibropapillomas and renal myxofibroma in a green turtle, Chelonia mydas. J Wildl Dis. (1990) 26:265–70. doi: 10.7589/0090-3558-26.2.265

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Page-Karjian A, Perrault JR, Zirkelbach B, Pescatore J, Riley R, Stadler M, et al. Tumor re-growth, case outcome, and tumor scoring systems in rehabilitated green turtles with fibropapillomatosis. Dis Aquat Organ. (2020) 137:101–8. doi: 10.3354/dao03426

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Work TM, Balazs GH, Rameyer RA, Morris RA. Retrospective pathology survey of green turtles Chelonia mydas with fibropapillomatosis in the Hawaiian Islands, 1993–2003. Dis Aquat Organ. (2004) 62:163–76. doi: 10.3354/dao062163

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Aguirre AA, Spraker TR, Chaves A, Toit L, Eure W, Balazs GH. Pathology of fibropapillomatosis in olive ridley turtles lepidochelys olivacea nesting in Costa Rica. J Aquat Anim Health. (1999) 11:283–9. doi: 10.1577/1548-8667(1999)011<0283:POFIOR>2.0.CO;2

CrossRef Full Text | Google Scholar

18. Herbst L. Fibropapillomatosis of marine turtles. Annu Rev Fish Dis. (1994) 4:389–425. doi: 10.1016/0959-8030(94)90037-X

CrossRef Full Text | Google Scholar

19. Work TM, Balazs GH. Relating tumor score to hematology in green turtles with fibropapillomatosis in Hawaii. J Wildl Dis. (1999) 35:804–7. doi: 10.7589/0090-3558-35.4.804

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Aguirre AA, Balazs GH. Blood biochemistry values of green turtles, Chelonia mydas, with and without fibropapillomatosis. Comp Haematol Int. (2014) 10:132–7. doi: 10.1007/s005800070004

CrossRef Full Text | Google Scholar

21. Hirama S, Ehrhart LM, Rea LD, Kiltie RA. Relating fibropapilloma tumor severity to blood parameters in green turtles Chelonia mydas. Dis Aquat Organ. (2014) 111:61–8. doi: 10.3354/dao02765

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Page-Karjian A, Norton TM, Krimer P, Groner M, Nelson SE Jr, Gottdenker NL, et al. Factors influencing survivorship of rehabilitating green sea turtles (Chelonia mydas) with fibropapillomatosis. J Zoo Wildl Med. (2014) 45:507–19. doi: 10.1638/2013-0132R1.1

PubMed Abstract | CrossRef Full Text | Google Scholar

23. IUCN. The IUCN Red List of Threatened Species. Version 2020-2 (2020). Available online at: https://www.iucnredlist.org.

Google Scholar

24. Herbst LH, Jacobson ER, Moretti R, Brown T, Sundberg JP, Klein PA. Experimental transmission of green turtle fibropapillomatosis using cell-free tumor extracts. Dis Aquat Organ. (1995) 22:1–12. doi: 10.3354/dao022001

CrossRef Full Text | Google Scholar

25. Jones K, Burgess G, Budd AM, Huerlimann R, Mashkour N, Ariel E. Molecular evidence for horizontal transmission of chelonid alphaherpesvirus 5 at green turtle (Chelonia mydas) foraging grounds in Queensland, Australia. PLoS One. (2020) 15:e0227268. doi: 10.1371/journal.pone.0227268

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Quackenbush SL, Work TM, Balazs GH, Casey RN, Rovnak J, Chaves A, et al. Three closely related herpesviruses are associated with fibropapillomatosis in marine turtles. Virology. (1998) 246:392–9. doi: 10.1006/viro.1998.9207

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Quackenbush SL, Casey RN, Murcek RJ, Paul TA, Work TM, Limpus CJ, et al. Quantitative analysis of herpesvirus sequences from normal tissue and fibropapillomas of marine turtles with real-time PCR. Virology. (2001) 287:105–11. doi: 10.1006/viro.2001.1023

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Brady G, Macarthur GJ, Farrell PJ. Epstein-barr virus and burkitt lymphoma. Postgrad Med J. (2008) 84:372–7. doi: 10.1136/jcp.2007.047977

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Boshoff C, Weiss RA. Epidemiology and pathogenesis of Kaposi's sarcoma-associated herpesvirus. Philos Trans R Soc Lond B Biol Sci. (2001) 356:517–34. doi: 10.1098/rstb.2000.0778

CrossRef Full Text | Google Scholar

30. Bosch FX, Manos MM, Munoz N, Sherman M, Jansen AM, Peto J, et al. Prevalence of human papillomavirus in cervical cancer: a worldwide perspective. International biological study on cervical cancer (IBSCC) Study Group. J Natl Cancer Inst. (1995) 87:796–802. doi: 10.1093/jnci/87.11.796

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Raiti P. Carbon dioxide (CO2) laser treatment of cutaneous papillomas in a common snapping turtle, Chelydra serpentina. J Zoo Wildl Med. (2008) 39:252–6. doi: 10.1638/2007-0055R.1

PubMed Abstract | CrossRef Full Text | Google Scholar

32. Duffy DJ, Martindale MQ. Perspectives on the expansion of human precision oncology and genomic approaches to sea turtle fibropapillomatosis. Commun Biol. (2019) 2:54. doi: 10.1038/s42003-019-0301-1

PubMed Abstract | CrossRef Full Text | Google Scholar

33. McAloose D, Newton AL. Wildlife cancer: a conservation perspective. Nat Rev Cancer. (2009) 9:517–26. doi: 10.1038/nrc2665

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Whilde J, Martindale MQ, Duffy DJ. Precision wildlife medicine: applications of the human-centred precision medicine revolution to species conservation. Glob Chang Biol. (2017) 23:1792–805. doi: 10.1111/gcb.13548

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Duffy DJ, Schnitzler C, Karpinski L, Thomas R, Whilde J, Eastman C, et al. Sea turtle fibropapilloma tumors share genomic drivers and therapeutic vulnerabilities with human cancers. Commun Biol. (2018) 1:63. doi: 10.1038/s42003-018-0059-x

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Yetsko K, Farrell JA, Blackburn NB, Whitmore L, Stammnitz MR, Whilde J, et al. Molecular characterization of a marine turtle tumor epizootic, profiling external, internal and postsurgical regrowth tumors. Commun Biol. (2021) 4:152. doi: 10.1038/s42003-021-01656-7

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Stacy BA, Foley AM, Work TM, Lauritsen AM, Schroeder BA, Hargrove SK, et al. Report of the Technical Expert Workshop: Developing Recommendations for Field Response, Captive Management, and Rehabilitation of Sea Turtles with Fibropapillomatosis. US Department of Commerce, National Marine Fisheries Service (NOAA Technical Memorandum NMFS OPR-60) (2018). p. 1–56.

Google Scholar

38. De Sousa Linhares A, Leitner J, Grabmeier-Pfistershammer K, Steinberger P. not all immune checkpoints are created equal. Front Immunol. (2018) 9:1909. doi: 10.3389/fimmu.2018.01909

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Waldman AD, Fritz JM, Lenardo MJ. A guide to cancer immunotherapy: from T cell basic science to clinical practice. Nat Rev Immunol. (2020) 20:651–68. doi: 10.1038/s41577-020-0306-5

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Treiner E, Duban L, Bahram S, Radosavljevic M, Wanner V, Tilloy F, et al. Selection of evolutionarily conserved mucosal-associated invariant T cells by MR1. Nature. (2003) 422:164–9. doi: 10.1038/nature01433

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Ussher JE, Klenerman P, Willberg CB. Mucosal-associated invariant T-cells: new players in anti-bacterial immunity. Front Immunol. (2014) 5:450. doi: 10.3389/fimmu.2014.00450

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Ussher JE, van Wilgenburg B, Hannaway RF, Ruustal K, Phalora P, Kurioka A, et al. TLR signaling in human antigen-presenting cells regulates MR1-dependent activation of MAIT cells. Eur J Immunol. (2016) 46:1600–14. doi: 10.1002/eji.201545969

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Kjer-Nielsen L, Patel O, Corbett AJ, Le Nours J, Meehan B, Liu L, et al. MR1 presents microbial vitamin B metabolites to MAIT cells. Nature. (2012) 491:717–23. doi: 10.1038/nature11605

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Constantinides MG, Link VM, Tamoutounour S, Wong AC, Perez-Chaparro PJ, Han SJ, et al. MAIT cells are imprinted by the microbiota in early life and promote tissue repair. Science. (2019) 366:eaax6624. doi: 10.1126/science.aax6624

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Pyagay P, Heroult M, Wang Q, Lehnert W, Belden J, Liaw L, et al. Collagen triple helix repeat containing 1, a novel secreted protein in injured and diseased arteries, inhibits collagen expression and promotes cell migration. Circ Res. (2005) 96:261–8. doi: 10.1161/01.RES.0000154262.07264.12

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Chen YL, Wang TH, Hsu HC, Yuan RH, Jeng YM. Overexpression of CTHRC1 in hepatocellular carcinoma promotes tumor invasion and predicts poor prognosis. PLoS One. (2013) 8:e70324. doi: 10.1371/journal.pone.0070324

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Tameda M, Sugimoto K, Shiraki K, Yamamoto N, Okamoto R, Usui M, et al. Collagen triple helix repeat containing 1 is overexpressed in hepatocellular carcinoma and promotes cell proliferation and motility. Int J Oncol. (2014) 45:541–8. doi: 10.3892/ijo.2014.2445

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Eriksson J, Le Joncour V, Nummela P, Jahkola T, Virolainen S, Laakkonen P, et al. Gene expression analyses of primary melanomas reveal CTHRC1 as an important player in melanoma progression. Oncotarget. (2016) 7:15065–92. doi: 10.18632/oncotarget.7604

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Kim JH, Baek TH, Yim HS, Kim KH, Jeong SH, Kang HB, et al. Collagen triple helix repeat containing-1 (CTHRC1) expression in invasive ductal carcinoma of the breast: the impact on prognosis and correlation to clinicopathologic features. Pathol Oncol Res. (2013) 19:731–7. doi: 10.1007/s12253-013-9636-y

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Kim HC, Kim YS, Oh HW, Kim K, Oh SS, Kim JT, et al. Collagen triple helix repeat containing 1 (CTHRC1) acts via ERK-dependent induction of MMP9 to promote invasion of colorectal cancer cells. Oncotarget. (2014) 5:519–29. doi: 10.18632/oncotarget.1714

PubMed Abstract | CrossRef Full Text | Google Scholar

51. Ma MZ, Zhuang C, Yang XM, Zhang ZZ, Ma H, Zhang WM, et al. CTHRC1 Acts as a prognostic factor and promotes invasiveness of gastrointestinal stromal tumors by activating Wnt/PCP-Rho signaling. Neoplasia. (2014) 16:265–78.e13. doi: 10.1016/j.neo.2014.03.001

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Wang CN, Li ZT, Shao F, Yang XY, Feng XL, Shi SS, et al. High expression of collagen triple helix repeat containing 1 (CTHRC1) facilitates progression of oesophageal squamous cell carcinoma through MAPK/MEK/ERK/FRA-1 activation. J Exp Clin Cancer Res. (2017) 36:84. doi: 10.1186/s13046-017-0555-8

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Egeblad M, Werb Z. New functions for the matrix metalloproteinases in cancer progression. Nat Rev Cancer. (2002) 2:161–74. doi: 10.1038/nrc745

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Zhang R, Cao Y, Bai L, Zhu C, Li R, He H, et al. The collagen triple helix repeat containing 1 facilitates hepatitis B virus-associated hepatocellular carcinoma progression by regulating multiple cellular factors and signal cascades. Mol Carcinog. (2015) 54:1554–66. doi: 10.1002/mc.22229

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Liu Z, Li LX, Yang ZX, Luo WR, Li X, Yang HL, et al. Increased expression of MMP9 is correlated with poor prognosis of nasopharyngeal carcinoma. BMC Cancer. (2010) 10:270. doi: 10.1186/1471-2407-10-270

PubMed Abstract | CrossRef Full Text | Google Scholar

56. He WL, Zhang H, Wang YF, Zhou YB, Luo YF, Cui YM, et al. CTHRC1 induces non-small cell lung cancer (NSCLC) invasion through upregulating MMP-7/MMP-9. BMC Cancer. (2018) 18:400. doi: 10.1186/s12885-018-4317-6

PubMed Abstract | CrossRef Full Text | Google Scholar

57. King M, Arnold JS, Shanske A, Morrow BE. T-genes and limb bud development. Am J Med Genet A. (2006) 140:1407–13. doi: 10.1002/ajmg.a.31250

CrossRef Full Text | Google Scholar

58. Tickle C. How the embryo makes a limb: determination, polarity and identity. J Anat. (2015) 227:418–30. doi: 10.1111/joa.12361

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Rodriguez-Esteban C, Tsukui T, Yonei S, Magallon J, Tamura K, Izpisua Belmonte JC. The T-box genes Tbx4 and Tbx5 regulate limb outgrowth and identity. Nature. (1999) 398:814–8. doi: 10.1038/19769

PubMed Abstract | CrossRef Full Text | Google Scholar

60. Takeuchi JK, Koshiba-Takeuchi K, Matsumoto K, Vogel-Hopker A, Naitoh-Matsuo M, Ogura K, et al. Tbx5 and Tbx4 genes determine the wing/leg identity of limb buds. Nature. (1999) 398:810–4. doi: 10.1038/19762

PubMed Abstract | CrossRef Full Text | Google Scholar

61. Horie M, Miyashita N, Mikami Y, Noguchi S, Yamauchi Y, Suzukawa M, et al. TBX4 is involved in the super-enhancer-driven transcriptional programs underlying features specific to lung fibroblasts. Am J Physiol Lung Cell Mol Physiol. (2018) 314:L177–L91. doi: 10.1152/ajplung.00193.2017

PubMed Abstract | CrossRef Full Text | Google Scholar

62. Lawrance MF, Mansfield KL, Sutton E, Savage AE. Molecular evolution of fibropapilloma-associated herpesviruses infecting juvenile green and loggerhead sea turtles. Virology. (2018) 521:190–7. doi: 10.1016/j.virol.2018.06.012

PubMed Abstract | CrossRef Full Text | Google Scholar

63. Yamamoto S, Nishimura O, Misaki K, Nishita M, Minami Y, Yonemura S, et al. Cthrc1 selectively activates the planar cell polarity pathway of Wnt signaling by stabilizing the Wnt-receptor complex. Dev Cell. (2008) 15:23–36. doi: 10.1016/j.devcel.2008.05.007

PubMed Abstract | CrossRef Full Text | Google Scholar

64. Zhang R, Lu H, Lyu YY, Yang XM, Zhu LY, Yang GD, et al. E6/E7-P53-POU2F1-CTHRC1 axis promotes cervical cancer metastasis and activates Wnt/PCP pathway. Sci Rep. (2017) 7:44744. doi: 10.1038/srep44744

PubMed Abstract | CrossRef Full Text | Google Scholar

65. Wang Q, Ding H, He Y, Li X, Cheng Y, Xu Q, et al. NLRC5 mediates cell proliferation, migration, and invasion by regulating the Wnt/beta-catenin signalling pathway in clear cell renal cell carcinoma. Cancer Lett. (2019) 444:9–19. doi: 10.1016/j.canlet.2018.11.024

PubMed Abstract | CrossRef Full Text | Google Scholar

66. Meissner TB, Li A, Biswas A, Lee KH, Liu YJ, Bayir E, et al. NLR family member NLRC5 is a transcriptional regulator of MHC class I genes. Proc Natl Acad Sci U S A. (2010) 107:13794–9. doi: 10.1073/pnas.1008684107

PubMed Abstract | CrossRef Full Text | Google Scholar

67. Rodriguez GM, Bobbala D, Serrano D, Mayhue M, Champagne A, Saucier C, et al. NLRC5 elicits antitumor immunity by enhancing processing and presentation of tumor antigens to CD8(+) T lymphocytes. Oncoimmunology. (2016) 5:e1151593. doi: 10.1080/2162402X.2016.1151593

PubMed Abstract | CrossRef Full Text | Google Scholar

68. Yoshihama S, Roszik J, Downs I, Meissner TB, Vijayan S, Chapuy B, et al. NLRC5/MHC class I transactivator is a target for immune evasion in cancer. Proc Natl Acad Sci U S A. (2016) 113:5999–6004. doi: 10.1073/pnas.1602069113

PubMed Abstract | CrossRef Full Text | Google Scholar

69. Ong CEB, Patchett AL, Darby JM, Chen J, Liu G-S, Lyons AB, et al. NLRC5 regulates expression of MHC-I and provides a target for anti-tumor immunity in transmissible cancers. bioRxiv [Preprint]. (2020). doi: 10.1101/2020.09.06.274720

CrossRef Full Text | Google Scholar

70. Zhang J, Cai H, Sun L, Zhan P, Chen M, Zhang F, et al. LGR5, a novel functional glioma stem cell marker, promotes EMT by activating the Wnt/beta-catenin pathway and predicts poor survival of glioma patients. J Exp Clin Cancer Res. (2018) 37:225. doi: 10.1186/s13046-018-0864-6

PubMed Abstract | CrossRef Full Text | Google Scholar

71. Hou MF, Chen PM, Chu PY. LGR5 overexpression confers poor relapse-free survival in breast cancer patients. BMC Cancer. (2018) 18:219. doi: 10.1186/s12885-018-4018-1

PubMed Abstract | CrossRef Full Text | Google Scholar

72. Rot S, Kaune T, Taubert H, Greither T, Kotrba J, Guttler A, et al. Prognostic impact of mRNA levels of LGR5 transcript variants in OSCC patients. BMC Cancer. (2019) 19:155. doi: 10.1186/s12885-019-5327-8

PubMed Abstract | CrossRef Full Text | Google Scholar

73. Dranoff G. Cytokines in cancer pathogenesis and cancer therapy. Nat Rev Cancer. (2004) 4:11–22. doi: 10.1038/nrc1252

PubMed Abstract | CrossRef Full Text | Google Scholar

74. Torphy RJ, Schulick RD, Zhu Y. Newly emerging immune checkpoints: promises for future cancer therapy. Int J Mol Sci. (2017) 18:2642. doi: 10.3390/ijms18122642

PubMed Abstract | CrossRef Full Text | Google Scholar

75. Flies AS Wild comparative immunology C. rewilding immunology. Science. (2020) 369:37–8. doi: 10.1126/science.abb8664

CrossRef Full Text | Google Scholar

76. Jacobson ER, Mansell JL, Sundberg JP, Hajjar L, Reichmann ME, Ehrhart LM, et al. Cutaneous fibropapillomas of green turtles (Chelonia mydas). J Comp Pathol. (1989) 101:39–52. doi: 10.1016/0021-9975(89)90075-3

PubMed Abstract | CrossRef Full Text | Google Scholar

77. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnetjournal. (2011) 17:10–2. doi: 10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

78. Wang Z, Pascual-Anaya J, Zadissa A, Li WQ, Niimura Y, Huang ZY, et al. The draft genomes of soft-shell turtle and green sea turtle yield insights into the development and evolution of the turtle-specific body plan. Nat Genet. (2013) 45:701–6. doi: 10.1038/ng.2615

PubMed Abstract | CrossRef Full Text | Google Scholar

79. Kim D, Langmead B, Salzberg SL. HISAT: a fast spliced aligner with low memory requirements. Nat Methods. (2015) 12:357–60. doi: 10.1038/nmeth.3317

PubMed Abstract | CrossRef Full Text | Google Scholar

80. Anders S, Pyl PT, Huber W. HTSeq-a Python framework to work with high-throughput sequencing data. Bioinformatics. (2014) 31:166–9. doi: 10.1093/bioinformatics/btu638

PubMed Abstract | CrossRef Full Text | Google Scholar

81. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. (2014) 15:550. doi: 10.1186/s13059-014-0550-8

PubMed Abstract | CrossRef Full Text | Google Scholar

82. Risso D, Ngai J, Speed TP, Dudoit S. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat Biotechnol. (2014) 32:896–902. doi: 10.1038/nbt.2931

PubMed Abstract | CrossRef Full Text | Google Scholar

83. Blighe K, Rana S, Lewis M. Enhanced Volcano: Publication-Ready Volcano Plots With Enhanced Colouring and Labeling. R package version (2019) 1.

84. Altschul SF, Gish W, Miller W, Myers EW, Lipman DJ. Basic local alignment search tool. J Mol Biol. (1990) 215:403–10. doi: 10.1016/S0022-2836(05)80360-2

CrossRef Full Text | Google Scholar

85. Altschul SF, Madden TL, Schaffer AA, Zhang J, Zhang Z, Miller W, et al. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res. (1997) 25:3389–402. doi: 10.1093/nar/25.17.3389

PubMed Abstract | CrossRef Full Text | Google Scholar

86. Camacho C, Coulouris G, Avagyan V, Ma N, Papadopoulos J, Bealer K, et al. BLAST+: architecture and applications. BMC Bioinformatics. (2009) 10:421. doi: 10.1186/1471-2105-10-421

CrossRef Full Text | Google Scholar

87. Mudunuri U, Che A, Yi M, Stephens RM. bioDBnet: the biological database network. Bioinformatics. (2009) 25:555–6. doi: 10.1093/bioinformatics/btn654

CrossRef Full Text | Google Scholar

88. Yu GC, Wang LG, Han YY, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. Omics. (2012) 16:284–7. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

89. Ackermann M, Koriabine M, Hartmann-Fritsch F, de Jong PJ, Lewis TD, Schetle N, et al. The genome of Chelonid herpesvirus 5 harbors atypical genes. PLoS One. (2012) 7:e46623. doi: 10.1371/journal.pone.0046623

PubMed Abstract | CrossRef Full Text | Google Scholar

90. Blackburn NB, Leandro AC, Nahvi N, Devlin MA, Leandro M, Escobedo IM, et al. Transcriptomic profiling of fibropapillomatosis in green sea turtles Chelonia mydas from South Texas. bioRxiv [Preprint]. (2020). doi: 10.1101/2020.10.29.360834

CrossRef Full Text | Google Scholar