Infectious vaccine-derived rubella viruses emerge, persist, and evolve in cutaneous granulomas of children with primary immunodeficiencies

Abstract

Rubella viruses (RV) have been found in an association with granulomas in children with primary immune deficiencies (PID). Here, we report the recovery and characterization of infectious immunodeficiency-related vaccine-derived rubella viruses (iVDRV) from diagnostic skin biopsies of four patients. Sequence evolution within PID hosts was studied by comparison of the complete genomic sequences of the iVDRVs with the genome of the vaccine virus RA27/3. The degree of divergence of each iVDRV correlated with the duration of persistence indicating continuous intrahost evolution. The evolution rates for synonymous and nonsynonymous substitutions were estimated to be 5.7 x 10-3 subs/site/year and 8.9 x 10-4 subs/site/year, respectively. Mutational spectra and signatures indicated a major role for APOBEC cytidine deaminases and a secondary role for ADAR adenosine deaminases in generating diversity of iVDRVs. The distributions of mutations across the genes and 3D hotspots for amino acid substitutions in the E1 glycoprotein identified regions that may be under positive selective pressure. Quasispecies diversity was higher in granulomas than in recovered infectious iVDRVs. Growth properties of iVDRVs were assessed in WI-38 fibroblast cultures. None of the iVDRV isolates showed complete reversion to wild type phenotype but the replicative and persistence characteristics of iVDRVs were different from those of the RA27/3 vaccine strain, making predictions of iVDRV transmissibility and teratogenicity difficult. However, detection of iVDRV RNA in nasopharyngeal specimen and poor neutralization of some iVDRV strains by sera from vaccinated persons suggests possible public health risks associated with iVDRV carriers. Detection of IgM antibody to RV in sera of two out of three patients may be a marker of virus persistence, potentially useful for identifying patients with iVDRV before development of lesions. Studies of the evolutionary dynamics of iVDRV during persistence will contribute to development of infection control strategies and antiviral therapies.

Fig 1. Phylogenetic tree of iVDRV.

The genetic relationships between the consensus genome sequences from each original granuloma sample and the whole genomes of the WHO reference viruses were inferred using the Maximum Likelihood method in MEGA7. All taxa are labeled with WHO names with iVDRV sequences marked with red dots. The genetic distances were computed using the Maximum Composite Likelihood method. The scale bar indicates the number of base substitutions per site. RA27/3 and iVDRVs represent a separate branch on the tree with RA27/3 being basal.

Fig 2. Diversity of iVDRV sequences.

a. Nucleotide and b. amino acid substitutions in iVDRV genomes relative to parental RA27/3 virus (GenBank #FJ211588) are depicted as vertical lines. The total number of substitutions are shown for each iVDRV on the right side. The consensus in b. represents amino acids identical in >50% sequences with the shared substitutions shown in red above the consensus; the indicated positions in the substitutions were the positions in the corresponding RV proteins. The positions of the coding sequences (CDS) and proteins are indicated by dark yellow and green pointed bars, respectively. The reference sequence is indicated by light yellow shading. The figure was prepared with Geneious software.

Fig 3. iVDRV quasispecies in tissue and viral isolate.

a. Nucleotide and amino acid substitutions in consensus RVi sequences relative to paired consensus RVs sequences shown by vertical lines for each case. The reference RVs sequences are indicated by light yellow shading. The number of substitutions are shown for each RVi on the right. The regions in the genomic sequence which encode proteins are indicated by the green pointed bars. b. Distribution of pairwise genetic distances between individual quasispecies within primary granuloma sample and the virus isolate from the CA patient. Each bar in the binned histogram represents the number of comparisons in each distance class. Note, the distance categories are not identical due to rounding. Underlying data can be found in the S4 Data file. c. Neighbor-joining tree (non-rooted) for quasispecies within the granuloma sample (n = 39, blue circles) and virus isolate (n = 30, red circles) from the CA case. The genetic distances were computed using the Maximum Composite Likelihood method. The scale bar indicates the number of base substitutions per site.

Fig 4. Analysis of iVDRV evolution patterns.

a. Relationship between the number of synonymous and nonsynonymous substitutions in the consensus genomic sequences directly from the specimens (RVs) and time elapsed after vaccination. b. Sequence divergence of the individual genes by a patient. dN/dS ratios were calculated with the SNAP tool. c. Plots of cumulative behavior of the average number of synonymous and nonsynonymous substitutions per a codon (cumulative substitution index) across each gene were calculated with the SNAP tool. Underlying data for Fig 4C can be found in the S5 Data file. The following domains are indicated: MT—methyltransferase, HVR -hypervariable region, X—ADP-ribose binding, P—protease, FL—fusion loops I and II, NT—neutralizing epitopes, TM—transmembrane domains. Note the NT domain is actually composed of four separate epitopes.

Fig 5. Different iVDRV lineages in different body sites (arm skin and NP).

a. Divergence of the iVDRV sequence in the NP fluids relative to iVDRV in the skin granuloma. Nucleotide and amino acid substitutions in the NP-derived sequence relative to granuloma-derived sequence (shaded in light yellow) are depicted as vertical lines. The number of substitutions are shown at the right side. The positions of encoded proteins in the genomic sequence are indicated by green pointed bars. b. Neighbor-Joining phylogenetic tree showing the genetic relationships between the whole genome sequences of RA27/3, LA-NP and LA-GR viruses. The genetic distances were computed using the Maximum Composite Likelihood method. The scale bar indicates the number of base substitutions per site.

Fig 6. RNA editing signatures in iVDRVs.

a. Spectra of single-nucleotide substitutions in iVDRV RVs genomes. Underlying data for Fig 6A–6C can be found in the S6 Data file. b. ADAR1 and ADAR2 prediction scores for adenines. The scores were calculated for 60 [A to_G] mutations in the positive strand and for 141 [A_to_G] mutations in the negative strand of the viral genome, as well as for the adenines that mutated to bases other than guanines along with adenines that did not mutate at all ([A–[A_to_G]]) (Note: Since, by convention, all mutations in S6 Data file are listed as changes in the positive strand, the negative strand [A to G] mutations are presented in this Table as [U to C] changes in the positive strand). The y axis denotes the prediction scores of hADAR1 and hADAR2 activity expressed as log2. The black horizontal line in the graphs denotes the median value of the prediction scores. P-values are shown above the scatter plots. c. Enrichments with APOBEC editing motifs in iVDRV RVs genomes. Mutated nucleotides are shown in capital letters within trinucleotide mutation motifs. Enrichment values were calculated for C to U mutations induced in the RV positive strand as described in Methods. P-values of one-sided Fisher’s exact test calculated as described in [37] can be found in the tab “Fisher_test_Fig6C” of S6 Data file.

Fig 7. Location of amino acid substitutions in 3D structure of E1 glycoprotein.

a. E1 monomer with all the identified substitutions mapped. The models of mutations were based on the coordinates from PDB entry 4adg. Mutational hotspots are shown encircled by a red line, the hotspot I including the residues I32, D34, K158, Q351, P415 and the hotspot II including the residues I50, V57, F84, V87, E118. The mutation A24->V is common in all iVDRV strains studied. The mutation F84->L was found in 5 of 6 strains and is located in the vicinity of the fusion loops. The color codes used: the variant residues found both in RVs and RVi (white spacefill), only in RVs (cyan spacefill), only in RVi (red spacefill); neutralizing epitopes NT1 (residues 225–235, brown), NT2 (245–251, orange), NT3 (260–266, violet) and NT4 (274–285, yellow); fusion loops (FL-I, residues 89–98, blue; FL-II, residues 130–138, magenta). b, c. Mutated amino acid residues on the surface of E1 monomer (b) and trimer (c) in iVDRVs (solvent accessibility >0.5, probe size 1.4A). d. The NT epitopes in E1 monomer. e. Possible escape mutations in the E1 neutralizing epitopes NT2, NT3 and NT4.

Fig 8. iVDRV growth properties.

a. Virus yields and percentages of RV-positive cells after infection of WI-38 with iVDRV, wtRV and RA27/3 strains at MOI = 5 (2 dpi) and MOI = 0.1 (3 dpi). Virus titers in the media were determined by titration on Vero cells, the number of infected cells was estimated by immunostaining for E1 protein. Data are presented as a mean +/- s.d. (n = 3, each experiment was performed in duplicate). b. Foci of infection of iVDRV isolates on Vero cells in comparison to wt and vaccine foci revealed by immunostaining for E1 at 6 dpi. c. Representative images from two independent experiments showing the results of the RNA-FISH for positive-strand RV RNA in WI-38 mock infected or infected with RA27/3 or RV-Dz at MOI = 5 (2 dpi). Nuclei were counterstained with DAPI. d. iVDRV persistence in WI-38. Growth curves of the indicated strains were shown for MOI = 5 and 0.1. Media were collected every 3–6 days and the extracellular viruses were titered on Vero cells. The representative results of two independent experiments each done in duplicate are shown. Limit of the assay detection (1x10² ffu/ml) is depicted by the dashed line. e. Phase contrast images of mock infected or infected (MOI = 5) cells at 36 dpi. Note cytopathic effects of RI6318 and LA3331 and the lack of live cells in RA27/3-infected wells. The adherent cells from one of two duplicate wells at 36 dpi were counted; the percentage of remaining adherent cells in each well was calculated relative to the mock infected well (yellow text). The cells in the second duplicate well were immunostained for E1.

Fig 9. Neutralization of iVDRV isolates by sera from immunized immunologically-competent individuals and PID patients.

Comparison of neutralization titers in 10 sera from healthy MMR vaccinees against iVDRVs and RA27/3. A neutralization titer was expressed as a log2 reciprocal of the serum dilution that protected 50% of the input virus. The cutoff of the assay (NT = 10) was depicted by the red dashed line.