This study was conducted according to the principles expressed in the Declaration of Helsinki. Ethics approval was obtained from the Ethics Committee of Northwestern and Central Switzerland (Project-ID 2020-00629). For all patients, either personal and/or family consent was obtained for autopsy and sample collection.
Patients and sample collection
The study is based on the analysis of 16 out of 21 consecutive COVID-19 autopsies performed between March 9th and April 14th 2020 at the Institute of Pathology Liestal and Institute of Medical Genetics and Pathology Basel, Switzerland. Clinical features including symptoms, course of disease, comorbidities, laboratory results and therapy are listed in Extended Data Table 1a. Detailed autopsy findings for each patient were recently published, and the identifiers (with the prefix “C”) for each COVID-19 patient are consistent with the description of this Swiss COVID-19 autopsy cohort3. In this study, we analysed formalin fixed and paraffin embedded (FFPE) lung tissue of distinct areas of the lungs of 16 of these COVID-19 patients. All 16 COVID-19 patients had positive nasopharyngeal swabs collected while alive. In all COVID-19 patients, diagnosis was confirmed by detection of SARS-CoV-2 in postmortal lung tissues. 5/16 patients were additionally tested by postmortal nasopharyngeal swabs which were positive for SARS-CoV-2 in all 5 cases.
As a control cohort, we selected 6 autopsies performed between January 2019 and March 2020 at the Institute of Pathology Liestal (“normal” patients N1 – N6). These control patients died of other, non-infectious causes and had a similar age, gender and cardiovascular risk profile. Patients with infections were excluded from this control cohort. Another control cohort consisted of 4 autopsies of patients suffering from various infections mainly with bacteria affecting the lung (patients with lung pathology, P1 – P4). Details for both control cohorts are listed in Extended Data Table 1b,c. SARS-CoV-2 was ruled out for each control patient by PCR-examination of lung tissue samples.
Nucleic acid extraction
RNA was extracted from up to six sections of FFPE tissue blocks using RecoverAll Total Nucleic Acid Isolation Kit (Cat No. AM1975, Thermo Fisher Scientific, Waltham, MA, USA). Extraction of DNA from up to 10 sections of FFPE tissue samples was automated by EZ1 Advanced XL (Qiagen, Hilden, Germany) using the EZ1 DNA Tissue Kit (Cat No. 953034, Qiagen, Hilden, Germany). Concentration of DNA and RNA were measured with Qubit 2.0 Fluorometer and Qubit dsDNA HS Assay or Qubit RNA HS Assay (Cat No. Q33230 & Q32852, Thermo Fisher Scientific, Waltham, MA, USA), respectively.
Quantification of SARS-CoV-2 in FFPE tissue samples
Post mortem viral load was individually measured in all lung tissue blocks from all patients included in this study. SARS-CoV-2 was detected in 15ng of human total RNA using the TaqMan 2019-nCoV Assay Kit v1 (Cat No. A47532, Thermo Fisher Scientific, Waltham, MA, USA), which targets three genomic regions (ORFab1, S Protein, N Protein) specific for SARS-CoV-2 and the human RNase P gene (RPPH1). The copy numbers of the SARS-CoV-2 viral genome was determined by utilizing the TaqMan 2019-nCoV Control Kit v1 (Cat No. A47533, Thermo Fisher Scientific, Waltham, MA, USA) and a comparative “ΔΔCт” method. The control kit contains a synthetic sample with a defined amount of target molecules for the human RPPH1 and the three SARS-CoV-2 assays, and was re-analyzed in parallel with patient samples. For each patient sample, this method resulted in individual copy numbers of the human RPPH1 and the three SARS-CoV-2 targets. Finally, the mean copy number of the SARS-CoV-2 targets was normalized to 1 x 106 RPPH1 transcripts.
Profiling of immune response by targeted RNAseq
The expression levels of 398 genes, including genes relevant in innate and adaptive immune response and housekeeping genes for normalization, were analyzed with the Oncomine Immune Response Research Assay (OIRRA, Cat No. A32881, Thermo Fisher Scientific, Waltham, MA, USA). The OIRRA is a targeted gene expression assay designed for the Ion™ next-generation sequencing (NGS) platform. This gene expression assay was originally designed to interrogate the tumor microenvironment to enable mechanistic studies and identification of predictive biomarkers for immunotherapy in cancer. The assay is optimized to measure the expression of genes involved in immune cell interactions and signaling, including genes expressed at low levels and involved in inflammatory signaling. The 398 genes covered by this assay are listed in Extended Data Table 2.
The NGS libraries were prepared as recommended by the supplier. In brief, 30ng of total RNA were used for reverse transcription (SuperScript VILO, Cat No. 11754250, Thermo Fisher Scientific, Waltham, MA, USA) and subsequent library preparation. The libraries were quantified (Ion Library TaqMan Quantitation Kit, Cat No. 4468802, Thermo Fisher Scientific, Waltham, MA, USA), equimolarly pooled and sequenced utilizing the Ion GeneStudio S5xl (Thermo Fisher Scientific, Waltham, MA, USA). De-multiplexing and gene expression level quantification were performed with the standard setting of the ImmuneResponseRNA plugin (version 220.127.116.11) within the Torrent Suite (version 5.12.1), provided as part of the OIRRA by Thermo Fisher Scientific, Waltham, MA, USA.
Detection of co-infections by whole genome sequencing
To identify potential pathogens accompanying an infection with SARS-CoV-2, we analyzed the DNA of the same tissue samples used for detection and profiling of the SARS-CoV-2-specific immune response. First, 250ng of genomic DNA was enzymatically sheared (15 minutes at 37°C) and barcoded using the Ion Xpress Plus Fragment Library Kit (Cat No. 4471269, Thermo Fisher Scientific, Waltham, MA, USA). Subsequently, the libraries were quantified (Ion Library TaqMan Quantitation Kit, Cat No. 4468802, Thermo Fisher Scientific, Waltham, MA, USA) and up to three libraries were pooled at equimolar levels for analysis with Ion GeneStudio S5xl (Thermo Fisher Scientific, Waltham, MA, USA). Sequencing data for each sample was analysed using the CLC genomics workbench (version 20.0.3, Qiagen, Hilden, Germany) in combination with the microbial genomics module (version 20.0.1, Qiagen, Hilden, Germany): The raw reads were trimmed by quality (Mott algorithm with limit 0.05 and a maximum of 2 ambiguous bases per read) and mapped to the human genome (GRCh37 hg19, match score: 1, mismatch cost: 2, indel opening cost: 6, indel extension cost: 1). Unmapped reads were analysed by taxonomic profiling to identify reads of viral or bacterial origin. The profiling utilized an index of 11’540 viral genomes with a minimum length of 1’000 bp and 2’715 bacterial reference genomes with a minimum length of 500’000 bp, retrieved from the NCBI Reference Sequence Database (date of download: 2020-04-02).
Immunohistochemical analyses for CD3, CD4, CD8, CD15, CD20, CD68, CD123, CD163, PD-1, MPO, p53, Ki67, C3d and C5b-9 were performed on all lung tissue blocks used in this study. Antibodies, staining protocols and conditions are detailed in Extended Data Table 4.
Qualitative and semiquantitative assessment of histopathological lung damage and neutrophilic infiltration
Hematoxylin and eosin (H&E) and Elastica van Gieson (EvG) stained sections of all lung tissues used in this study were independently evaluated by two experienced and board certified pathologists (VZ and KDM) (Extended Data Table 5). Both pathologists evaluated the presence of diffuse alveolar damage (DAD), and if present, its stage, intra-alveolar edema and hemorrhage. In addition, both pathologists evaluated the severity of histopathological changes in COVID-19 lungs (1 = mild / discrete alterations, 2 = moderate, 3 = severe changes) based on resemblance between normal and pathologically altered lung tissues. Parameters that were taken into account included reduction of alveolar air-filled spaces, typical histologic features of DAD with hyaline membrane formation, infiltration of lymphocytes, monocytes and neutrophils into interstitial and alveolar spaces, type 2 pneumocyte hyperplasia, desquamation of pneumocytes, histologic features of organizing pneumonia including intra-alveolar fibrin deposition and fibrosis (acute fibrinous and organizing pneumonia, AFOP)17,18. The number of neutrophils per lung tissue section was estimated on H&E stained sections and by immunohistochemical stains for CD15 and MPO using a three tiered system (1 = few or no neutrophils, 2 = moderate number of neutrophils, 3 = high number of neutrophils). Assessment of the two pathologists was concordant in the vast majority of cases. Discrepant cases were reviewed by a third pathologist (NW) to reach consent.
Digital image analysis
Slides were digitalized on a 3DHistech™ P1000 slide scanner at 400x magnification (3DHISTECH Ltd. Budapest, Hungary). Digital slide review and quality control was performed by a board-certified pathologist (VHK). Tissue regions with staining artefacts, folds or other technical artefacts were excluded from analysis. A deep neural network (DNN) algorithm (Simoyan and Zisserman VGG, HALO AI™ on HALO™ 3.0.311.167, Indica Labs, Corrales, NM) was trained using pathologist annotations to automatically localize and measure the area of each lung tissue sample on the digital slides. Background regions and glass were excluded from analysis. Mark-up images for tissue classification were generated and classification accuracy was confirmed through pathology review. For cell-level analysis, color deconvolution for DAB, AP and hematoxylin channels was performed and nuclear segmentation was optimized using cell-morphometric parameters. Marker-positive cells in stromal and epithelial regions were quantified. For CD3, CD4, CD8, CD20, CD68, CD123, CD163 and PD1, staining detection was optimized for the cytoplasmic / membranous compartment and marker expression was measured on a continuous scale at single cell resolution. For assessment of CD8/PD1 double stains, color deconvolution was optimized for separation of DAB (PD1) and AP (CD8) staining products. Internal controls (non-immune cells) and external controls (tonsil) were used to calibrate the detection limits and cross-validated by visual review. For each tissue sample, the total area of lung tissue in mm2, the absolute number of marker‐positive cells, cell morphometric parameters and staining intensity were recorded.
Identification of SARS-CoV-2 immune response pattern
Gene expression analysis
Samples were included in the study based on quality of libraries and alignment performance. Applied inclusion criteria are: >1 million of mapped reads, good concentration of libraries, average read length >100bp, > 300 target genes with more than 10 reads. One sample with > 1 Mio reads was excluded from the study because of shorter read length and a low library concentration. Notably this sample had the longest time between death and autopsy (72h) before analysis. Differential expression analysis was performed using the edgeR package comparing normal lung samples, COVID-19 samples and samples from patients with other infections. Genes were selected for downstream analyses by fdr <0.05 and |logFC| >1 for clustering analysis. Clustering analysis was performed using k-means algorithm and complete linkage. Ideal number of clusters (n=3) was chosen based on 30 different algorithms19 and the final clustering derives from the consensus of 2000 iterations. Expression of gene signatures was calculated as median of log2(cpm + 1) of selected genes.
Functional enrichment analysis
Biological processes enrichment was performed using the enrichGO function of the package clusterProfiler20 setting all the genes included in the assay as universe.
All the analyses and graphical representations were performed using the R statistical environment software21 and the following packages: ggplot222, circlize23, ComplexHeatmap24, ggfortify25, reshape226 and factoextra27. Correlation between transcripts and viral counts was performed using Pearson's correlation. Association between continuous and categorical data were tested using Wilcoxon rank sum test.
Box-plots elements indicate the median (center line), upper and lower quartiles (box limits) and show all the data points. Whiskers extend to the most extreme value included in 1.5x interquartile range.