High-depth sequencing characterization of viral dynamics across tissues in fatal COVID-19 reveals compartmentalized infection

Abstract

SARS-CoV-2 distribution and circulation dynamics are not well understood due to challenges in assessing genomic data from tissue samples. We develop experimental and computational workflows for high-depth viral sequencing and high-resolution genomic analyses from formalin-fixed, paraffin-embedded tissues and apply them to 120 specimens from six subjects with fatal COVID-19. To varying degrees, viral RNA is present in extrapulmonary tissues from all subjects. The majority of the 180 viral variants identified within subjects are unique to individual tissue samples. We find more high-frequency (>10%) minor variants in subjects with a longer disease course, with one subject harboring ten such variants, exclusively in extrapulmonary tissues. One tissue-specific high-frequency variant was a nonsynonymous mutation in the furin-cleavage site of the spike protein. Our findings suggest adaptation and/or compartmentalized infection, illuminating the basis of extrapulmonary COVID-19 symptoms and potential for viral reservoirs, and have broad utility for investigating human pathogens.

Fig. 1. Overview of study design and sample selection.

A Overview of sample selection and analysis of the FFPE tissue specimens (B) Normalized SARS-CoV-2 viral load in lung samples from a cohort of 39 subjects from which samples in the study were selected (median lung viral load is highlighted for six subjects of interest). C Summary of selected individuals sequenced in this study, including clinical characteristics and viral strain information. Abbreviations: OSA (obstructive sleep apnea); HF (heart failure); HTN (hypertension); DM (diabetes mellitus); CAD (coronary artery disease); COPD (chronic obstructive pulmonary disease); ALL (acute lymphocytic leukemia); BM (bone marrow); CKD (chronic kidney disease); RA/SLE (rheumatoid arthritis/systemic lupus erythematosus); ISL (interstitial lung disease); MGUS (monoclonal gammopathy of undetermined significance); PVD (peripheral vascular disease). The ‘*’ designates uncertainty around the time between symptom onset and death. D Boxplot (top) and heatmap (bottom) each demonstrate normalized SARS-CoV-2 quantification across tissue samples available for four or more subjects, and testis. In the boxplot, boxes delineate quartiles and whiskers show the range of all samples available (3-6 subjects). In the heatmap, composite samples (those where 2 or more tissues were in the same FFPE block) are designated with an asterisk (Supplementary Fig. 1b and Methods), gray represents virus not detected, and white designates sample was not available. E Representative IHC sections from study samples show the presence of SARS-CoV-2 protein (brown) in multiple tissues, including pneumocytes in lung, ciliated respiratory epithelium in the trachea, cardiomyocytes in the heart, hepatocytes in the liver, small intestine epithelial cells, rete testis, and tubular epithelium in kidneys. Staining was performed once. Scale bars are 20-microns. These findings are in agreement with sequencing data.

Fig. 2. Refined sequencing methods and robust variant calling methodology enables confident analysis of variants from autopsy tissues.

A Schematization of enhanced sequencing methods. B Two independent libraries from four samples with high viral depth of coverage (>3,000x mean coverage) were downsampled to a range of mean coverage depths (90x to 1,250x), and variant profiles were identified in each condition, then compared to those detected in the highest coverage condition (>3,000x mean coverage). For each library from all four samples, at each coverage depth condition, precision and recall (top; points represent the mean, while error bars represent standard deviation), the number of variants identified (middle; boxes delineate quartiles, whiskers delineate range excluding outliers), and the frequency distribution of variants (bottom; points represent variants across all samples within a condition) were compared to the highest depth of coverage condition. Green represents conditions >500x mean depth of coverage, the threshold selected for high resolution genomic analyses. C For 12 samples, the same libraries were sequenced with and without hybrid capture enrichment, and were then downsampled to the same number of raw reads; mean depth of viral coverage was calculated and plotted for each sample. The order of magnitude “OM” of enrichment is annotated across the two-dimensional space. D Variants were identified for all samples in the sample set with at least 500x mean depth of coverage, including those sequenced with hybrid capture enrichment (purple) and without (blue). Frequency of variants identified in two independent libraries were compared (top), demonstrating high correlation (R² = 0.998). Number of variants identified (middle) and frequency of variants identified (bottom) was compared across samples, each showing poor correlation with viral load (R² = 0.229 and R² = 0.167, respectively).

Fig. 3. Complete SARS-CoV-2 genomes across subjects and tissues.

A A maximum likelihood phylogenetic tree of the unique viral genomes identified, including two unique consensus genomes from S02, demonstrates genetic divergence. The viral sequences are presented in the context of 729 sequences from a six month window centered on the time these infections occurred, with color representing SARS-CoV-2 Pango lineage (legend). Inset, the SNP that differentiates unique genomes from S02 (A24292G) is highlighted; allele frequencies across all seven tissues from which a genome assembled are annotated (for two samples, black outline designates >500x mean depth coverage for high-confidence variant quantification). B A schematic representation of complete viral genomes assembled across the six subjects. Colored shapes represent the assembly of a complete genome, each unique color-shape combination represents a unique genome, which is represented in the phylogenetic tree. A black outline designates >500x mean depth coverage, suitable for minor variant and transcriptomic analysis. C Principal component analysis of viral fragment abundance for 32 unique samples (those with >500x mean depth of viral coverage, some being duplicate libraries), colored subject, demonstrate that samples separate in PC1 by time between symptom onset and death. D Nucleocapsid (N) gene expression, normalized by the total number of viral reads, decreased with time between symptom onset and death. Legends indicate day between symptom onset and death, and whether a sample underwent targeted enrichment.

Fig. 4. SARS-CoV-2 minor variants across tissues.

A Number of variants and viral population diversity (Shannon entropy) for each sample with >500x depth of coverage, arranged by subject (ordered by time between symptom onset and death). In the boxplot, boxes delineate quartiles and whiskers show the range, excluding outliers (as determined by the interquartile range). S01 and S02 had only two samples for comparison and were not analyzed further for circulation or compartmentalization; the remainder of the subjects had 4–11 samples. The ‘*’ designates uncertainty around the time between symptom onset and death. B On the left, all variants are displayed as genome position vs frequency for each of subjects S03-S06. The size of the points reflects the number of tissues the variant was observed in; the standard deviation of the frequency at which the variant occurred across tissues is depicted by error bars. Red points reflect nonsynonymous changes whereas grey points reflect synonymous or noncoding variants. Dashed lines are present at 10 and 90% frequency; variants falling between dashed lines were considered high-frequency. On the right, high-frequency variants were quantified and diagrams were constructed to demonstrate genetic distance. Dashed line represents the number of SNPs differentiating the consensus genome from SARS-CoV-2 Wuhan reference strain (NC_045512.2).