Canine OMM samples (n = 8) were obtained following surgical resection as treatment at Animal Medical Center, Yamaguchi University (Table 1). The bleed of all dogs used in this study was Dachshund. The patient's owners were informed before sample collection. Total RNA was isolated from fresh OMM from eight dogs using RNA Premium Kit (FastGene #FG-81050). We collected mRNA by NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB #E7490). Strand-specific libraries were generated by NEBNext Ultra Directional RNA Library Prep Kit for Illumina (NEB #E7420) and indexed by NEBNext Multiplex Oligos for Illumina (NEB #E7335). All libraries were sequenced on Illumina NextSeq 500. We also obtained transcriptome sequencing data of canine OMM and oral healthy tissues (8 and 3 samples, respectively) from NCBI sequence read archives (https://www.ncbi.nlm.nih.gov/sra) (20) to confirm our RNA-seq analyses. Accession numbers and read mapping results were summarized in Supplementary Table 1.

All paired-end FASTQ files were trimmed, and their sequence qualities were checked using Trim Galore (version 0.6.4) (https://www.bioinformatics.babraham.ac.uk/projects/trim_galore/). Then, the filtered short reads were aligned to the canine reference genome (un-masked canFam3 genome from UCSC genome browser: https://hgdownload.soe.ucsc.edu/goldenPath/canFam3/bigZips/) using HISAT2 (version 2.1.0) (21). For each RNA-seq sample, we conducted transcript assemblies using StringTie2 (version 2.1.1) (21) with “–rf” options to specify strandedness. Resultant GTF files obtained from the eight canine OMM samples and the canine RefSeq GTF file (canFam3, downloaded from UCSC genome browser on June 15, 2021) were merged using StringTie2 with “–merge” option to generate a merged GTF file. The “-m 500 -f 0.05” options were combined to improve the specificity of the transcript assembly. Each melanoma transcript identical to a RefSeq transcript was allocated to the RefSeq nomenclature. The resultant transcript assemblies are available in Supplementary File 1 (a compressed GTF file).

A BED file of canine ERV-ORFs was downloaded from gEVE (version 1.1) (http://geve.med.u-tokai.ac.jp/) (14). In gEVE, pol genes that were thought to be derived from long interspersed nuclear elements (LINEs) were included (14). To remove LINE-derived sequences, ERV-ORFs which were annotated as “LINE” by RepeatMasker and/or “YP_073558.1” or “NP_048132.1” by BLASTP against the NCBI Viral Genome Database were removed. Transcripts in the merged GTF file overlapped with ERV-ORFs were identified using BEDtools intersect with “-s” option to consider strandedness (22). Forty-three ERV-derived genes detected at this step were listed in Supplementary Table 2. All ERV-derived genes were manually checked using Integrative Genomics Viewer (version 2.4.9) (23). Transcripts per kilobase million (TPM) was calculated using StringTie2 with “-e –rf” options using the merged GTF file. Eleven ERV-derived genes with the mean TPM > 1 were listed including Canis-env2 (Table 2). Differentially expressed genes were identified and visualized with an MA plot and a volcano plot using DESeq2 (version 1.26.0) (24) in R (version 3.6.1).

Total RNA was extracted from cultured OMM cells from four animals (Supplementary Table 3), which were maintained in Dulbecco's modified Eagle medium (Sigma-Aldrich #D5796) supplemented with 10% heat-inactivated fetal bovine serum (Gibco #10270), penicillin (100 IU/ml) and streptomycin (100 μg/ml) (Nacalai tesuque # 09367-34). To synthesize the cDNA, RNA was reverse transcribed using Verso cDNA Synthesis kit (Thermo Fisher Scientific #AB1453B). For negative controls, samples without reverse transcriptase were prepared. PCR was performed using KOD One (TOYOBO #KMM-101). Primer pairs used for NMG-1 and NMG-4 were 5′-ACCCGCTGACTATGACTCAGGAAC-3′ (forward) and 5′-AACAGTCTTTGCCTCTGCTGTCAG-3′ (reverse), and 5′-AGGCACTCCTCCACGCCACTACTAG-3′ (forward) and 5′-TCAAGTCTGTGCTTATGTAGGGACCAG-3′ (reverse), respectively. For loading control, GAPDH was amplified with a primer pair, 5′-AAGGTCGGAGTGAACGGATTTGG-3′ (forward) and 5′-CGGTTGCTGTAGCCAAATTCATTG-3′ (reverse). After preheating at 94°C for 2 min, thermal cycling reaction (98°C for 10 s and 68°C for 10 s) was repeated 35 times (for NMG-1 and NMG-4) or 30 times (for GAPDH). The sequences of the amplicons were checked by the Sanger sequencing (Fasmac Co., Ltd.).

To collect Canis-env2 orthologs in Carnivora genomes, coding sequences of Canis-env2 in the dog was searched using Blat (25) in the following represented species in the UCSC genome browser (https://genome.ucsc.edu): Southern Sea otter, enhLutNer1 (Jun. 2019); Ferret, musFur1 (Apr. 2011); Hawaiian monk seal, neoSch1 (Jun. 2017); panda, ailMel1 (Dec. 2009); and Cat, felCat9 (Nov. 2017). The hits with the highest similarity scores were confirmed to be located between BTN2A1 and BTN1A1. To conduct the following evolutionary analyses of Canis-env2, we used MEGA-X (26). We first aligned nucleotide sequences of Canis-env2 genes by MUSCLE program with “Align codons” option (27), and pairwise p-distance values of amino acids, which are proportions of amino acid sites at which the two sequences to be compared are different, were calculated. The numbers of non-synonymous substitutions (dN) and synonymous substitutions (dS) per site were estimated by Nei-Gojobori method with Jukes-Canter correlation (28).

We conducted RNA-seq analyses using canine oral malignant melanoma (OMM) cases (n = 8). To identify ERV-derived genes from these sequence data, we constructed a pipeline to identify ERV-ORFs expressed in canine melanoma (Figure 1A). First, we constructed de novo transcript assemblies with reference to the RefSeq gene coordinates using StringTie2. In this procedure, transcripts expressed from the overlapped locus were given a single gene name. In total, we obtained 106,209 transcripts from the 8 OMM RNA-seq samples, and these transcripts were organized into 34,292 genes, of which 3,035 genes were not reported in the reference gene annotation (Figure 1B). Next, we extracted genes whose transcripts were overlapped with ERV-ORFs in the gEVE database (http://geve.med.u-tokai.ac.jp/) (14) (Figure 1C). As a result, we obtained 43 genes, of which five genes were newly identified (Figure 1D, Supplementary Table 2).

Next, we examined the structures of these 43 genes overlapped with ERV-ORFs. In nine genes, ERV-ORFs were identical to full or partial coding sequences of host genes: ASPRV1, CAR1, CNBP, CTSD, GIN1, PEG10, RTL1, SUGP2, and TAF1 (Figure 2). In these genes, four genes—ASPRV1, CAR1, PEG10, and RTL1 — were known to originate from ERVs or LTR-type retrotransposons and mammalian or Carnivora-specific genes. ASPRV1 codes an aspartic protease conserved in mammals and has roles in skin maintenance (29). PEG10 and RTL1 are gag-pol genes derived from Ty3/Gypsy retrotransposon and were involved in placental development (30, 31). CAR1, also known as syncytin-Car1, is a retroviral env gene conserved in Carnivora (32). Syncytin-Car1 protein is involved in the cell fusion of syncytiotrophoblast during placentation (32). We used these four definite retroviral genes for further analysis. The remaining five genes—CNBP, CTSD, GIN1, SUGP2, and TAF1—consist of several exons, and one exon identified as containing ERV-ORFs (Figure 2A). Orthologs of these exons were found in at least chickens, mice, and humans and may be derived from ERVs before the divergence of birds and mammals. This possibility has been particularly suggested for GIN1, where its integrase domain shows a high similarity to the GIN element, an animal-specific DNA transposon that encodes a Gypsy/Ty3 retrotransposon-like integrase (33). On the other hand, to the best of our knowledge, no evidence has been reported to indicate that CNBP, CTSD, SUGP2, and TAF1 were derived from transposons. It is also plausible that these genes independently acquired protein motifs similar to ERVs. Since it is beyond the focus of this study to make conclusions, we did not regard these five genes as ERV-derived genes in this study. Two splicing variants of CRADD were overlapped with two ERV-ORFs, although they were in the 5′ UTR. Therefore, CRADD was also removed for further analyses. Others were uncharacterized genes in the RefSeq database (i.e., genes whose names start with LOC), including non-coding genes or newly identified non-RefSeq genes. To identify sufficiently expressing genes, we collected genes that were expressed in the eight canine OMM samples with an average TPM > 1. Finally, we identified 11 ERV-derived genes, including the newly identified melanoma genes 1 to 4 (NMG-1 to 4), which were expressed in canine oral melanoma (Table 2). Since NMG-1 to 4 overlapped with RepeatMasker track of an ERV group, CfERV1-int (Figure 2B; Supplementary Figure 1), these genes are presumed to be derived from ERVs.

Next, we identified the ERV-derived genes that showed higher expression levels in the OMM samples than those in healthy tissues. Although we did not conduct the transcriptome sequencing of control samples, we utilized publicly available RNA-seq data of eight canine OMM samples and three canine oral healthy samples, all of which were obtained from the SRA database (20) (Supplementary Table 1). Then, we found that four ERV-derived genes were significantly upregulated (PEG10, NMG-1, NMG-4, and LOC102155597), and two ERV-derived genes were down-regulated (ASPRV1 and LOC10215559) in OMM (adjusted p < 0.05) (Figures 3A,B, Supplementary Table 4). TPMs of these four upregulated genes varied among melanoma samples, but their expression levels were low in all healthy tissues (Figure 3C). PEG10 was reported as an important gene in several species of human cancers, such as hepatocellular carcinoma (34), breast cancer (35), lung cancer (36), and neuroendocrine prostate cancer (37). NMG-1 and NMG-4 are newly identified genes, and their expressions were experimentally verified by RT-PCR (Figure 3D). NMG-1 gene has a single transcript isoform that is composed of two exons, and its second exon was overlapped with four ERV-ORFs annotated by gEVE database (Figure 2B). NMG-4 gene has two transcripts, and both are composed of three exons. LOC102155597 was overlapped with an Env-like ERV-ORF of 447 amino acid length (Figure 2B). This env gene was previously designated as Canis-env2 (32). We found that Canis-env2 homologs located between BTN2A1 and BTN1A1 among six Carnivora species available in UCSC genome browser, and they were presumed to be orthologs (Figure 4A). We found that the coding sequences of Canis-env2 orthologs were under purifying selection, suggesting that the proteins are functionally important (Figure 4B). This Canis-Env2 protein retains the signal peptide, the furin cleavage site, and the immunosuppressive domain, but not the transmembrane (TM) domain (Figure 4C). Together, we successfully identified newly identified or uncharacterized retroviral coding genes expressed in canine OMM.