Coevolutionary Analysis Implicates Toll-Like Receptor 9 in Papillomavirus Restriction

Abstract

Upon infection, DNA viruses can be sensed by pattern recognition receptors (PRRs), leading to the activation of type I and III interferons to block infection. Therefore, viruses must inhibit these signaling pathways, avoid being detected, or both. Papillomavirus virions are trafficked from early endosomes to the Golgi apparatus and wait for the onset of mitosis to complete nuclear entry. This unique subcellular trafficking strategy avoids detection by cytoplasmic PRRs, a property that may contribute to the establishment of infection. However, as the capsid uncoats within acidic endosomal compartments, the viral DNA may be exposed to detection by Toll-like receptor 9 (TLR9). In this study, we characterized two new papillomaviruses from bats and used molecular archeology to demonstrate that their genomes altered their nucleotide compositions to avoid detection by TLR9, providing evidence that TLR9 acts as a PRR during papillomavirus infection. Furthermore, we showed that TLR9, like other components of the innate immune system, is under evolutionary selection in bats, providing the first direct evidence for coevolution between papillomaviruses and their hosts. Finally, we demonstrated that the cancer-associated human papillomaviruses show a reduction in CpG dinucleotides within a TLR9 recognition complex. IMPORTANCE Viruses must avoid detection by the innate immune system. In this study, we characterized two new papillomaviruses from bats and used molecular archeology to demonstrate that their genomes altered their nucleotide compositions to avoid detection by TLR9, providing evidence that TLR9 acts as a PRR during papillomavirus infection. Furthermore, we demonstrated that TLR9, like other components of the innate immune system, is under evolutionary selection in bats, providing the first direct evidence for coevolution between papillomaviruses and their hosts.

Keywords: Mexican free-tailed bat; Papillomaviridae; TLR9; evolutionary biology; innate immunity; papillomavirus; speciation.

FIG 1

Evolutionary relationship of novel bat papillomaviruses. Shown is a maximum-likelihood phylogenetic tree inferred using concatenated E1, E2, and L1 protein sequences. Papillomaviruses associated with Chiroptera are highlighted: Yangochiroptera (orange) and Yinpterochiroptera (green). Papillomavirus genera are collapsed (the number of types within each genus is indicated in parentheses). Bootstrap-generated branch support values are given using symbols and color gradients. Host species are shown using Sonoran Desert-dwelling animals. The red arrow indicates the subtree used for further analyses as described throughout the paper. Inset plots the percentage pairwise identities across the L1 nucleotide sequences.

FIG 2

Coevolution of papillomaviruses. (A) Optimized tanglegram between subtree based on concatenated E1-E2-L1 maximum likelihood phylogenetic tree (Fig. 1) and associated host species. We downloaded the host species tree from www.timetree.org. Papillomaviruses are linked to their host phylogenies. Papillomaviruses associated with Chiroptera are highlighted: Yangochiroptera (orange) and Yinpterochiroptera (blue). The underlined virus/host pairs were used for the analysis in panel C. (B) Procrustean Approach to Cophylogeny analysis based on the interaction network and phylogenies shown in panel A supports the notion that papillomaviruses coevolved with their hosts. The observed best‐fit Procrustean superimposition (red dotted line) lies outside the 95% confidence interval (shaded area of the curves) of the distribution of network randomizations in the null model. (C) As in panel B but using a subset of the interaction network and phylogenies (underlined virus-host pairs in panel A).

FIG 3

CpG dinucleotide sequences are significantly depleted in papillomavirus genomes. The observed-versus-expected (O/E) ratios of each dinucleotide in the papillomavirus genome sequences shown in Fig. 2 were calculated using a custom wrapper around the CompSeq program from the Emboss software suite. The red line at 1.0 indicates the ratio where a dinucleotide is seen as often as expected by chance.

FIG 4

CpG content is significantly lower in papillomaviruses associated with Yangochiroptera than in related viruses. (A) A maximum likelihood phylogenetic tree is shown comparing the O/E ratios of CpG dinucleotides. Viruses infecting Yangochiroptera (red), Yinpterochiroptera (green), and related hosts (gray) are indicated. (B) Mean (±standard deviation) CpG O/E ratios for each group of viruses were compared using one-way ANOVA with Tukey’s post hoc test. (C) CpG O/E ratios are compared to total GC content for each viral genome in panel A.

FIG 5

Yangochiroptera papillomaviruses have a restricted codon usage. (A) Codon usage tables for each virus in Fig. 2 were compared using the “codcmp” program from the Emboss software suite. Root-mean-square deviation (RMSD) values for each pairwise comparison are plotted as box-and-whisker plots with the outliers (colored circles) identified using Tukey’s method. Individual values are shown as single black dots. (B) The codon adaptation indices for a subset of viruses (see Materials and Methods) in Fig. 2 are plotted as box-and-whisker plots, with the outliers (colored circles) identified using Tukey’s method. Individual values are shown as single black dots. (C) RSCU values for the indicated amino acid/codons were calculated and plotted as box-and-whisker plots, with the outliers (colored circles) identified using Tukey’s method. RSCU values for each amino acid were compared using two-way ANOVA with Tukey’s post hoc test. *, P < 0.05; **, P < 0.01; ***, P < 0.001. (D) The amino acid composition was calculated as described in Materials and Methods. Mean values ± standard deviations are plotted.

FIG 6

Yangochiroptera TLR9 is evolving under diversifying selection. (A) Structure of horse TLR9 in complex with agonistic DNA (PDB code 3WPC) (94). Amino acids of interest are highlighted. (B) Maximum likelihood phylogenetic tree of mammalian TLR9 sequences clusters Yangochiroptera and Yinpterochiroptera separately from the mammalian TLR9. Red branches display evidence of episodic diversifying selection as identified by aBSREL (92). Alignments show sequences of interest. The sequence logo is based on the alignment of 29 nonchiropteran TLR9 sequences. Numbering is based on the mouse TLR9. The analyses identified residues indicated by “$” as being selected using FEL. Residues indicated by “*” were previously identified as evolving under diversifying selection (66), while those indicated by “#” were identified as functionally important through site-directed mutagenesis (94).

FIG 7

Yangochiroptera papillomaviruses depleted a TLR9 recognition motif from their genomes. (A) The O/E ratios of each nCGn tetramer in the Yangochiroptera papillomavirus genomes sequences were calculated using a custom wrapper around the CompSeq program from the Emboss software suite. Mean values ± standard deviations are plotted. (B) The O/E ratios of each nCGn tetramer in the different groups were calculated as for panel A. The proportion of these ratios provides a normalized view of tetramer depletion across papillomavirus genomes shown in Fig. 2. (C) The Yangochiroptera versus Yinpterochiroptera nCGn proportions (as in panel B) are plotted as brown dots and compared to 1,000 randomly shuffled sequences (green violin) plots. Only ACGT is statistically underrepresented in the Yangochiroptera.

FIG 8

Human papillomaviruses depleted a TLR9 recognition motif from their genomes. (A) The O/E ratios for statistically underrepresented nCGn tetramers are plotted as box-and-whisker plots, with the outliers (circles) identified using Tukey’s method. All known HPVs within the Alphapapillomavirus genus are shown. We also separately analyzed the mucosal and cutaneous HPVs in this genus. (B) Statistical depletion of nCGn tetramers was based on comparison of distributions of the actual to randomized sequences. The distributions for TCGA and TCGT (n = 65 alphapapillomavirus genomes) are shown. (C) The O/E ratios for statistically underrepresented nnCGnn hexamers are plotted. Individual values are plotted as dots. Numbers near the x axis indicate the number of genomes in which the hexamer was not observed. (D) Statistical depletion of nnCGnn hexamers was based on comparison of distributions of the actual (red curves) to randomized (yellow curves) sequences. The distributions for the 10 hexamers identified in panel C (n = 65 alphapapillomavirus genomes) are shown.