Tissue tropism and functional adaptation of the SARS-CoV-2 spike protein in a fatal case of COVID-19

Abstract

Systemic spread of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) to extrapulmonary tissues has been observed following acute infections. Autopsy studies further indicate tissue-specific virus diversity, including in immune-privileged sites. Questions remain on the viral dynamics leading to the tissue tropism of SARS-CoV-2, including evolutionary trajectories and functional adaptations that could impact persistence and transmission. In this study, we characterized SARS-CoV-2 genomes from 27 distinct tissues collected from an autopsy case where the patient had a primary immune deficiency. We identified tissue-specific virus genotypes, in some instances coexisting within the same sites, with mutations primarily in the receptor-binding domain of the spike protein. Protein simulations and isolation of infectious virus indicate combinations of spike substitutions that would lead to increased protein stability and stronger binding of the virus to host cells. This highlights the importance of studying patients with weakened immune responses where potential tissue reservoirs provide an environment permissive for SARS-CoV-2 evolution and diversification.IMPORTANCEPersistent severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infections in immunocompromised individuals are considered a potential source of new viral variants. Beyond the respiratory tract, the virus can spread within days to organs like the brain, heart, and kidneys, where distinct tissue microenvironments may further drive viral evolution and the emergence of new mutations. In this study, we compared the genetic diversity of SARS-CoV-2 genomic RNA isolated from 27 distinct tissue sites collected from an individual with a weakened immune system. By linking viral population dynamics across these tissue sites, we defined the extent of compartmentalization during multi-organ spread, highlighting how non-respiratory tissues can impact SARS-CoV-2 diversification.

Keywords: coronavirus; tissue tropism; viral evolution; viral intrahost diversity.

Fig 1

Multiple SARS-CoV-2 genotypes are identified across tissues in an autopsy case. (A) Maximum-likelihood phylogenetic tree of 27 consensus sequences from the autopsy case and 6,000 randomly selected AY.119 consensus sequences circulating in the USA between 1 August 2021 and 31 January 2022. Color indicates if the sample was collected in Florida (FL, orange), New York (NY, blue), a state outside of NY or FL (no point), or our autopsy case (green). Outlier sequences are excluded from the tree. (B) The proportion of 6,000 randomly sampled AY.119 sequences from GISAID collected in the USA from 1 August 2021 to 31 January 2022. Color indicates if the sample was collected and sequenced in New York (NY, purple), Florida (FL, orange), or a different state (gray). (C) The TreeTime output of the root-to-tip distance versus the collection date. Outliers and autopsy samples were ignored (red points). (D) Maximum-likelihood phylogenetic tree of the SARS-CoV-2 virus consensus sequences collected and sequenced from 27 different tissue sites. The x-axis represents the estimated number of substitutions from the shared ancestral node.

Fig 2

Consensus diversity is primarily in the spike coding region. (A) Consensus mutations along the SARS-CoV-2 genome. Mutations are determined by comparing each nucleotide position to the AY.119 consensus reference (Materials and Methods). Each mutation is labeled with the coding region: AY.119 nucleotide–nucleotide position–variant nucleotide. The number of samples with the variant and the amino acid (aa) information (AY.119 aa–aa position–variant aa) is provided in parentheses. Both synonymous and nonsynonymous substitutions are labeled. Open reading frames and coding regions are represented by different colored boxes; only regions with variants are labeled. UTR, untranslated region; ORF, open reading frame; 3a, open reading frame 3a (ORF3a). (B) The estimated spike haplotype (H) frequency (y-axis) calculated using CliqueSNV for each virus sequenced from each tissue site (x-axis). Only spike haplotypes predicted to be present at ≥10% are provided. Color indicates the spike haplotype (H1–H10) identified in each tissue site, with the amino acid residues for each haplotype outlined to the right. Only nonsynonymous amino acid substitutions are provided. Asterisks (*) mark haplotypes found in multiple samples. See also Table 1 and Fig. S2 and S3.

Fig 3

Most combinations of spike substitutions are associated with increased protein stability and higher binding energy to the ACE2 receptor. (A) Structural mapping of individual spike substitutions on the 3D AY.119 spike-ACE2 complex. Substitutions are represented as the AY.119 amino acid, amino acid position, and variant amino acid. The spike domains are represented by color. NTD, N-terminal domain (periwinkle); RBD, receptor-binding domain (brown); SD1, sub-domain 1 (blue); and SD2, sub-domain 2 (blue). (B, C) Distributions of averaged root mean square deviations (RMSD) (y-axis) for the (B) complete spike-ACE2 complex or the (C) NTD and RBD of the spike across all spike sequences (x-axis). Average RMSD values were calculated for 100 individual frames at every nanosecond of an 80 ns molecular dynamics simulation time (n = 8,000 per spike haplotype). Global RMSD values (B) for all mutant haplotypes were significantly different from the AY.119 reference (adj. P-value <0.0001, Kruskal-Wallis). Color indicates the haplotype and matches the haplotype colors outlined in Fig. 2. All boxplots represent the median (middle line), first and third quartiles (box), 1.5 * interquartile range (whiskers), and outliers (points). (D) The mean and standard deviation of the change in binding free energy (kcal/mol) (x-axis) between the entire spike structure and ACE2 for each spike sequence (y-axis) using molecular mechanics Poisson-Boltzmann surface area. Non-significant differences between AY.119 and each mutant spike haplotype are noted with ns, while significant differences are marked with * and indicate a P-value <0.0001 using a one-way ANOVA test. (E) A heatmap of the decomposed energetic interactions for each residue (y-axis) where a consensus variant was found in our data set. AY.119 reference residues include R19, K417, V445, G446, Y453, L455, G476, S477, K478, Q493, and P561. Only variant residues are outlined for each mutant spike haplotype along the x-axis. If not provided, the residue is the reference AY.119 residue. White boxes indicate the residue is significantly involved in the total binding energy. See also Fig. S4.

Fig 4

Combinations of spike substitutions impact the binding and internalization of infectious isolates. (A) The log₁₀ fold change in SARS-CoV-2 RdRp gene expression for each plaque-purified isolate (1–5) with respect to mock infections for binding (left) and internalization (right) assays at an MOI of 0.1 (top) and 1 (bottom) in VeroE6 cell lines constitutively expressing human ACE2-TMPRSS2 or TMPRSS2 only (n = 3 for each isolate and assay). Data were normalized using GAPDH gene expression, and changes in RdRp gene expression were calculated and compared to those of the mock controls. P-values were calculated by comparing each isolate’s binding and internalization measurements to those of isolate 2 or by comparing measurements of each isolate within ACE2-TMPRSS2 and TMPRSS2 cells. Only P-values <0.01 are shown in the figure (**P < 0.01, ***P < 0.001). Color indicates the different infectious plaque-purified isolates. The mutant spike sequence is provided for each isolate and includes the AY.119 amino acid residue, amino acid position, and variant amino acid for each variant site. The most similar haplotype sequence found in Fig. 2 is provided in parentheses next to each isolate sequence. See also Tables S2 and S3. (B) Distributions of averaged RMSD (y-axis) for the complete spike-ACE2 complex across the spike isolates (1–4 only) and the AY.119 reference (x-axis). Average RMSD values were calculated for 100 individual frames at every nanosecond of an 80 ns molecular dynamics simulation time (n = 8,000). Boxplots represent the median (middle line), first and third quartiles (box), 1.5 * interquartile range (whiskers), and outliers (points). Non-significant differences between AY.119 and each isolate’s spike haplotype are noted with ns, while significant differences (P-value <0.0001, one-way ANOVA test) are marked with *. (C) The mean and standard deviation of the change in binding free energy (kcal/mol) (x-axis) between the entire spike structure and ACE2 for each spike sequence (y-axis) using molecular mechanics Poisson-Boltzmann surface area.

Fig 5

Tissue compartmentalization and mixing of minor variant populations. (A) A network visualizing the hierarchical clustering of pairwise BCI values. Links represent a BCI less than the mean BCI for the data set (0.66). Red lines represent connections that span clusters, and black lines indicate connections within a cluster. Each node is colored by the general tissue location. Each cluster is numbered (1–3), and the mean and standard deviation of the intra- and inter-cluster comparisons are provided. BCI values assume independence between variants and were calculated using positions where a variant was present as a minor variant in at least one sample. (B) Boxplots of the number of minor (≥5%–49.99%) SNVs (y-axis) for the Bray-Curtis clusters. Each point represents a single tissue site. Color represents the general tissue location and is maintained for the figure. (C) The nonsynonymous (nonsyn) and synonymous (syn) minor variant divergence for each nucleotide site across the genome, normalized by the total number of available nonsynonymous and synonymous sites, and excluding the untranslated and non-coding regions of the genome. Data are by Bray-Curtis cluster. Each line connects tissue site values. Asterisks represent significance (P-values: ***<0.001) between noted comparisons using the Mann-Whitney U-test. (D) The nucleotide variants found in more than two tissue sites that define or are shared by the three Bray-Curtis clusters. Fill indicates the relative frequency of the given variant (x-axis) within each tissue site (y-axis). Variants along the x-axis are ordered by their cluster association and the total number of samples in which they were identified. Asterisks indicate whether a mutation is unique to cluster 1 (*) or cluster 2 (**). NSP, non-structural protein; ORF, open reading frame. See also Fig. S6.