Molecular consequences of SARS-CoV-2 liver tropism

1In this study, we examined three patient cohorts to characterize severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2)-associated liver injury and infection. First, a clinical cohort (Hamburg; cohort 1, n = 99) was used to examine changes in liver function tests (LFTs) in patients admitted to hospital due to coronavirus disease 2019 (COVID-19) and those who were diagnosed with COVID-19 during hospitalization (Fig. 1a); 42% of these patients were at least 65 years old, 70% were males, only 23% had at least 3 coexisting conditions (Source Data Fig. 1). Among patients admitted due to COVID-19 (n = 72), raised aspartate aminotransferase (AST) and alanine aminotransferase (ALT) were shown in 63 and 39% of cases (Fig. 1b), respectively despite the low frequency of liver disease (1.4%), confirming previous studies–. Next, we analysed a subgroup of patients that were admitted to hospital due to alternative diagnoses (n = 27) and suffered from nosocomial SARS-CoV-2 infection. Both AST and ALT were significantly increased after the diagnosis of COVID-19 (Fig. 1c). While the percentage of patients with raised ALT almost doubled before and after COVID-19 diagnosis (33 versus 59%), the percentage of patients with elevated AST almost tripled (22 versus 67%) (Fig. 1d). Next, we evaluated a large validation cohort of patients admitted to hospital due to COVID-19 (Michigan; cohort 2, n = 1,219) (Fig. 1e); 44% of patients were at least 65 years old, 57% were males and only 16% had at least 3 coexisting conditions (Supplementary Table 1). Importantly, 57% (n = 699) of patients showed elevations of AST and 37% (n = 452) of ALT at admission (Fig. 1e). Furthermore, LFT elevations at admission and during the second week of hospitalization were associated with mortality (Fig. 1f), raising questions about their role in disease severity. Together, these observations clearly highlight hepatic injury as an important clinical feature of patients with COVID-19 requiring hospitalization.

2It has been postulated that LFT elevations in hospitalized patients with COVID-19 may result from systemic inflammation or severe cellular stress (for example, hypoxia), as generally observed in critically ill patients. However, an autopsy study reported ultrastructural evidence of SARS-CoV-2 (ref. ), and a second study reported histopathological findings in a liver biopsy of a patient with abnormal liver enzymes, including no obvious inflammation in the portal area, with normal description of the interlobular bile duct, interlobular vein, interlobular artery and hepatocytes with minimal inflammatory cell infiltration. Given that liver biopsies are not routinely performed in patients with COVID-19 with altered LFTs, we evaluated a third cohort (cohort 3, n = 45 autopsy cases) in search of direct evidence of liver infection (Supplementary Table 2); 73% (n = 33) were older than 65, 69% (n = 31) had at least 3 coexisting conditions and 62% (n = 28) were males, matching the demographic characteristics linked to severe COVID-19. SARS-CoV-2 RNA was detected using quantitative PCR with reverse transcription (RT–qPCR) targeting the E gene in 69% of cases (n = 31) (Fig. 2a), which was more frequently associated with older age, male sex and multiple coexisting conditions, as shown in previous reports,,.

3Pathology assessment only revealed signs of shock, probably due to severe systemic disease, and adipose changes without classical signs of hepatitis or significant immune infiltration (n = 18; Extended Data Fig. 1). However, a combination of indirect immunofluorescence and high-resolution confocal microscopy identified hepatic expression of angiotensin-converting enzyme 2 (ACE2), the main cellular receptor of SARS-CoV-2 (Fig. 2b and Extended Data Fig. 2; n = 5 liver samples), and SARS-CoV-2 spike protein (Fig. 2c and Extended Data Fig. 3, n = 3 examples). In addition, SARS-CoV-2 RNA was also found in hepatic cells with in situ hybridization (Supplementary Fig. 1; n = 4, that is, n = 2 SARS-CoV-2 PCR+ and n = 2 controls (Ctrls)), confirming the spatial location of SARS-CoV-2 in the human liver using multiple techniques that complement the findings of previous studies,. To provide additional context, we first compared the levels of SARS-CoV-2 RNA copies per cell between airway samples (previously reported in Puelles et al.) and a subset of autopsy livers, showing a similar range but a lower median viral RNA load in hepatic specimens (Supplementary Fig. 2a; n = 18 lung, n = 16 pharynx and n = 18 liver), and we quantified the number of SARS-CoV-2 spike-positive cells in a small subset of our autopsy samples (n = 6 liver samples, that is, n = 3 SARS-CoV-2 PCR+ and n = 3 Ctrls with quantification provided in Supplementary Fig. 2b), both of which may explain the lack of overt cytopathic changes (that is, hepatitis). Notably, we isolated infectious SARS-CoV-2 from 2 out of 3 autopsy livers and corresponding lungs, confirming the infectivity of postmortem liver tissue samples, leading to a 1,000 time increase in viral RNA after 48 h of cell infection in vitro (Fig. 2d), as previously reported in the kidney. Together, our data provide spatial and functional evidence of SARS-CoV-2 hepatic tropism, at least in cases of fatal COVID-19.

4Next, we aimed to determine the molecular changes associated with SARS-CoV-2 liver tropism. Thus, liver tissues from a subset of COVID-19 autopsies were selected for transcriptomic and proteomic analyses (n = 11; Fig. 2e). Briefly, 3 comparison groups were defined a priori: (1) patients with COVID-19 and SARS-CoV-2-positive livers (n = 5 for transcriptomics and n = 6 for proteomics); (2) patients with COVID-19 and SARS-CoV-2-negative livers (n = 5 for transcriptomics and proteomics) to account for the systemic effect of COVID-19; and (3) patients with non-COVID-19-related deaths (n = 5 for proteomics and n = 7 for transcriptomics, that is, Ctrls) to account for pleiotropic effects resulting from autopsy-related collection and storage. The groups were similar in terms of age, sex and number of coexisting conditions (Supplementary Table 3). First, transcriptomic profiling confirmed the expression of known SARS-CoV-2 entry receptors and facilitators in the liver autopsy samples (Extended Data Fig. 4), including ACE2, transmembrane protease serine 2 (TMPRSS2), procathepsin L and Ras-related protein Rab-7a,. Importantly, a direct comparison to a publicly available dataset contextualized the expression of these facilitators in relation to other organs, where liver samples showed similar expression levels as the lung, only surpassed by the kidney and small intestine, suggesting that hepatic tissue is susceptible to SARS-CoV-2 infection. Then, we detected SARS-CoV-2 subgenomic RNA within our own dataset and identified two samples out of six with values in the range of the infected respiratory tract, tracing SARS-CoV-2 to the liver by another independent methodology (Extended Data Fig. 5).

5Transcriptional changes associated with SARS-CoV-2 hepatic tropism revealed a significant upregulation of type I and II interferon (IFN) responses (Fig. 2f and Extended Data Fig. 6). Gene Ontology (GO) analysis confirmed upregulation of IFN and viral responses and increases in lipid and phospholipid metabolic processes (Extended Data Fig. 7), supporting a recent study linking perturbations of lipid metabolism with SARS-CoV-2 infection. Importantly, these findings were also consistent at a single patient level (Fig. 2g) and within single gene sets (Fig. 2h,i). While SARS-CoV-2 has developed strategies to evade IFN responses, leading to relatively weak type I and III IFN induction in host cells, a recent study showed a high expression level of IFN-stimulated genes in patients with a high SARS-CoV-2 viral load and high levels of pro-inflammatory cytokines but relatively intact lung morphology. These observations are in agreement with our findings, supporting upregulated IFN responses in tissues with traces of SARS-CoV-2 RNA and without classical cytopathic changes. In addition, we also identified that components of IFN-related JAK-STAT signalling (Supplementary Fig. 3), and other central cellular processes involved in lipid/cholesterol synthesis, including terpenoid and steroid biosynthesis were all modulated in patients with SARS-CoV-2 liver tropism. Altogether, these data provide a comprehensive framework to understand the molecular consequences of SARS-CoV-2-mediated liver injury.

6Global proteomic profiling in hepatic tissues of the same autopsy cohort showed substantial regulation of IFN, viral and in particular SARS-CoV-2 responses (Fig. 3a and Source Data Fig. 3). Protein ontology analysis revealed upregulation of IFN signalling and downregulation of basic biological processes, including oxidation reduction, oxidative phosphorylation and cellular respiration (Fig. 3b). Type I IFN responses were confirmed at the single patient level (Fig. 3c). Gene/protein sets for hallmark ‘interferon-α (IFN-α) response’ (R = 0.65, P < 0.05; Fig. 3d) and the GO Biological Process ‘defence response to virus’ (R = 0.71, P < 0.05; Fig. 3e), as well as single gene/protein analysis for IFN-inducible genes (R = 0.74–0.87, P < 0.05; Fig. 3f), showed strong associations between transcriptomic and proteomic profiles, highlighting a coherent modulation of SARS-CoV-2-related pathways at the transcriptional and translational level. Furthermore, transcription factor activity profiling (Fig. 3g) provided individual (Fig. 3h) and global (Fig. 3i) signatures that combined transcription factor upregulation (for example, signal transducer and activator of transcription 2 (STAT2), a signal transducer and activator of transcription that mediates signalling by type I IFNs) and downregulation (that is, nuclear receptor corepressor 2 (N-CoR2), which aids histone deacetylases to modify chromatin structure, thus preventing basal transcriptional activity). Canonical pathway analysis confirmed upregulation of immune-mediated, IFN and JAK-STAT-related responses (Extended Data Fig. 8a), as well as downregulation of glucose and carbohydrate metabolism in association with liver tropism (Extended Data Fig. 8b). Interestingly, modulation of JAK-STAT signalling has been proposed to target SARS-CoV-2 viral entry and replication in a study that dissected SARS-CoV-2 infectivity using organotypic three-dimensional cultures of primary human liver cells.

7The clinical and molecular data from our study suggest that SARS-CoV-2-mediated liver injury is comparable to previously characterized hepatotropic viruses. For this reason, we compared publicly available datasets containing bulk RNA sequencing (RNA-seq) data from liver tissue samples from patients with hepatitis C virus (HCV), human immunodeficiency virus (HIV) and hepatitis B virus (HBV) via gene set enrichment analysis (GSEA) (Fig. 4a). While our literature search was carefully performed, it was not possible to find studies that could be perfectly matched our study demographics (for example, our population was notably older than the other three studies as per the demographic trends of fatal COVID-19; Supplementary Table 4). Interestingly, and despite of multiple methodological limitations, all three infections affecting the liver, including HBV (Fig. 4b), HCV (Fig. 4c) and HIV (Fig. 4d), showed overlapping gene expression patterns with our proposed SARS-CoV-2 signature. Specific examples include upregulation of IFN signalling gene sets (Fig. 4e), especially types I, II and III IFN responses (Fig. 4f), and downregulation of viral gene expression and related cellular functions (that is, translational initiation and peptide metabolic processes) (Fig. 4g). Pathway analysis also revealed remarkable overlap between SARS-CoV-2 and HCV (13 of 20 pathways for upregulated; Fig. 4h; 8 of 20 pathways for downregulated; Fig. 4i). In addition, we identified similar transcription factor use in livers from patients with severe COVID-19 and patients with liver steatosis and HIV, inactive carriers of HBV and chronic HCV (Fig. 4j). Taken together, this comparative analysis highlights a shared molecular signature between SARS-CoV-2 and multiple viruses with well-characterized patterns of liver injury.

8Since the topic of SARS-CoV-2 hepatic tropism has been highly debated, we aimed to tackle this question from multiple perspectives, including clinical (two large independent patient cohorts that showcase COVID-19-associated liver injury), histopathological (identification of SARS-CoV-2 at the RNA and protein levels with high spatial resolution), virological (RT–qPCR in human autopsy samples and isolation of SARS-CoV-2 from human autopsy liver specimens), molecular (transcriptomic, proteomic and transcription factor activity use profiling, indicating the presence of SARS-CoV-2 in association with significant tissue responses) and bioinformatic (detection of subgenomic RNA and data validation with external public datasets of virus-mediated hepatic injury). Despite all these layers of evidence supporting SARS-CoV-2 hepatic tropism, the exact mechanism of infection is unclear. Our findings show that the human liver is susceptible to SARS-CoV-2 infection since cell entry receptors and facilitators are clearly expressed in hepatic cells. However, we observed a mismatch between the expression of the ACE2 protein, mostly observed in Kupffer cells, and the location of the SARS-CoV-2 spike protein, which was observed in Kupffer and parenchymal cells (for example, hepatocytes). In this context, it is worth considering a recent study that proposed high-density lipoprotein scavenger receptor class B member 1 (SRB1) as a cell entry facilitator for SARS-CoV-2. SRB1 acts as a critical receptor that affects HCV entry, an infection that shares some molecular features with SARS-CoV-2 liver tropism. While Wei et al. reported strong protein expression of SRB1 in the human liver, we confirmed messenger RNA expression of SCARB1 (SRB1 gene) in multiple human organs, showing the highest expression levels in the liver followed by the brain, small intestine, lung, kidney and muscle (Extended Data Fig. 9). In addition, we confirmed SRB1 protein expression on human hepatocytes and in association with the expression of SARS-CoV-2 spike protein (Extended Data Fig. 10). It is plausible that even if the expression of ACE2 in hepatocytes is relatively low at both RNA and protein levels, high levels of SRB1 may enhance the potential for SARS-CoV-2 entry.

9It is worth remembering that our observations represent a late ‘snapshot’ of a severe disease state and consequent multi-organ dysfunction (that is, severe COVID-19 leading to death), which provides a limited view affected by disease course, coexisting conditions, therapeutic interventions and postmortem processes. All these factors may influence our view of receptor expression patterns (for example, potential up/downregulation during infection), definition of cell-specific tropism (that is, severe injury to infected cells and local protein clearance) and even molecular signatures, which likely reflect profound cellular changes in neighbouring cells that have survived an active phase of SARS-CoV-2 tropism and could explain the common features between patients with severe COVID-19 and hepatotropic viruses (that is, HCV). Furthermore, our findings should be carefully considered in the context of lack of access to full clinical data in our large clinical cohorts, limiting our ability to rule out other causes of liver injury (for example, drug-mediated hepatoxicity), and relatively small sample sizes in our multi-omics experiments, which were mainly focused on patients who died due to severe COVID-19, raising questions about generalizability to mild and moderate cases. Despite this, a recent clinical study followed 443 non-hospitalized individuals for approximately 9.6 months after their first positive SARS-CoV-2 test, revealing signs of multi-organ injury, which allows us to question the potential for similar molecular changes in mild and moderate forms of COVID-19.

10In summary, our findings provide multimodal evidence of SARS-CoV-2 human liver tropism and a molecular signature, featuring a strong upregulation of IFN responses, JAK-STAT signalling and liver-specific metabolic modulation. This unique human dataset can be used as a rational framework for future studies aiming to characterize the potential consequences of SARS-CoV-2 extrapulmonary manifestations and identify personalized therapeutic strategies.

11Autopsies were performed at the Institute of Legal Medicine of the University Medical Center Hamburg-Eppendorf. From every liver specimen collected by the Institute of Legal Medicine, multiple randomly chosen small samples were available for different analyses. The ethics committee of the Hamburg Chamber of Physicians was informed about the study (nos. 2020-10353-BO-ff and PV7311). Ctrls included cases of sudden, non-infectious deaths. The postmortem interval was on average 6 d. Informed consent was obtained from a next of kin or legal representative for autopsy and tissue sampling. No compensation was paid.

12The study protocol for clinical data collection (patient cohorts) was approved by the institutional review board (IRB) of the University of Michigan (no. HUM00178971) and Hamburg (no. WF-052/20). The IRB approved a waiver of informed consent for this observational study.

13Tissue samples were systematically sampled during the autopsy procedure in 45 individuals. Automated nucleic acid extraction was performed according to the manufacturer’s recommendations with whole process control (control kit; Roche), with a final elution volume of 100 μl. For virus quantification, a previously published assay was adopted with modifications using chimeric 2′-O-methyl RNA bases at the penultimate base of both primers (mG and mC) to reduce primer dimer formation. The forward primer 5′-ACAGGTACGTTAATAGTTAATAGCmGT-3′ (400 nM end concentration), 5′-TATTGCAGCAGTACGCACAmCA-3′ (400 nM end concentration) and probe 5′-Fam-ACACTAGCC/ZEN/ATCCTTACTGCGCTTCG-Iowa Black FQ-3′ (100 nM end concentration) were used. Primer and probes were obtained from Integrated DNA Technologies. One-step RT–PCR (25 μl volume) was run on the LightCycler 480 system (Roche) using the one-step RNA control kit as master mix (Roche) and 5 μl of eluate. The Ct value for the target SARS-CoV-2 RNA (FAM) was determined using the second derivative maximum method. To quantify the standard in vitro-transcribed RNA, the E gene of SARS-CoV-2 was used. The standard was obtained via the European Virus Archive. The linear range of the assay was between 1 × 103 and 1 × 109 copies ml−1. β-Globin qPCR was performed with the commercial TaqMan primer set (catalogue no. 401846; Thermo Fisher Scientific) and Roche DNA control kit. The PCR was run on the LightCycler 480 system. The amount of DNA was normalized using a human DNA standard (KR0454). SARS-CoV-2 RNA levels in tissues were normalized to β-globin DNA.

14Liver tissue was fixed in formaldehyde. Then, 2-μm slides were cut and haematoxylin and eosin, periodic acid–Schiff and Masson–Goldner stainings were performed. Expert pathology review was conducted by two independent and experienced pathologists.

15ACE2 and SARS-CoV-2 spike were detected in formalin-fixed paraffin-embedded sections using a protocol for indirect immunofluorescence and confocal microscopy,. We used primary antibodies against ACE2 (dilution: 1:200, catalogue no. AF933; R&D Systems) and SARS-CoV/SARS-CoV-2 (COVID-19) (dilution 1:200, catalogue no. GTX632604; GeneTex), which were validated in previous studies,. We also used a primary antibody against SRB1 (dilution 1:200, catalogue no. ab217318; Abcam).

16For quantification, and considering the size and degree of autolysis, we performed targeted sampling of five fields of view based on the presence of at least one SARS-CoV-2 spike+ cell with a random location within the field; for Ctrls, five random fields were chosen from non-autolytic sites. Then, we quantified the number of SARS-CoV-2 spike+ cells and the total number of complete nuclei per field, which allowed us to calculate a percentage of positive cells per field. Statistics were performed with Prism v.9.2.0 (GraphPad Software).

17In situ hybridization was carried out to detect virus RNA of SARS-CoV-2 on paraffin sections utilizing the RNAscope 2.5 HD detection kit (catalogue no. 322310; Advanced Cell Diagnostics) according to the manufacturer’s instructions. Briefly, tissue sections were deparaffinized in xylene followed by target retrieval at 95 °C for 10 min. Subsequently, internal peroxidase activity was quenched by hydrogen peroxide incubation for 10 min followed by permeabilization using protease plus treatment at 40 °C for 30 min. The SARS-CoV-2-specific RNAscope probe V-nCoV2019-S (catalogue no. 848561) was hybridized at 40 °C for 2 h. RNAscope probes specific for either the human ubiquitin C mRNA (catalogue no. 310041) or the bacterial dihydrodipicolinate reductase mRNA (catalogue no. 310043) were used as positive or negative Ctrl, respectively. The RNAscope signal was developed with 3,3′-diaminobenzidine and nuclei were counterstained with haematoxylin.

18Liver tissues were homogenized and 250 μl of the homogenized tissue solution were used to infect Vero cells (CRL-1586; ATCC). Growth was confirmed by RT–qPCR of cell culture supernatants.

19Sample selection was based on collection date and RNA quality. RNA was extracted with the QIAGEN RNeasy Micro Kit according to the manufacturer’s instructions followed by Agilent Bioanalyzer sample quality control, library preparation with Lexogen CORALL Total RNA and ribosomal RNA depletion. Single-end RNA-seq was done using 75-base pair NextSeq v.2.5 with >30 million reads. Raw reads were trimmed with Trim Galore! v.0.4.3 (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) and aligned with RNA STAR: Galaxy Version 2.7.2b. Mapped reads were counted with HTSeq: Galaxy Version 0.9.1..

20Raw counts from several healthy tissues were downloaded directly from the GTEX portal (www.gtexportal.org) on 4 March 2020.

21Genes with a non-null count in at least two samples were selected. Raw counts were processed with the limma R package v.3.42.2 (ref. ) and differentially regulated genes between PCR-positive and PCR-negative groups were identified after library size normalization. Both groups were compared to GTEX liver. A Benjamini–Hochberg-adjusted P < 0.05 was considered significant.

22GSEA was performed using the GAGE R package v.2.36.0 (ref. ) with the Molecular Signatures Database (MSigDB) gene set v.7.1, comparing PCR-positive to PCR-negative groups. Both groups were also compared with GTEX liver. A Benjamini–Hochberg-adjusted P < 0.05 was considered significant.

23First, each sample was normalized to healthy liver tissue to depreciate the background noise and emphasize the true signal. Therefore, the average GTEX liver transcripts per kilobase million (TPM) intensity was subtracted for every single gene. In each sample, genes were ranked according to their normalized TPM value and used as input for single-sample enrichment analysis using the fgsea R package v.1.12.0 (ref. ) with MSigDB gene set v.7.1 (ref. ) and ConsensusPathDB v.34 (ref. ). In all gene sets, the enrichment score was calculated based on a random walk over the ranked list of genes. The significance of the enrichment scores was assessed with 1,000 permutations. A Benjamini–Hochberg-adjusted P < 0.05 was considered significant.

24Kyoto Encyclopedia of Genes and Genomes pathways were created with PathView v.1.34.0 (ref. ) with log2 positive versus negative fold change.

25The relative fraction of single-guide RNAs (sgRNAs) was determined by identifying the leader sgRNA fusion containing sequencing reads normalized to read counts containing the genomic RNA leader sequence. Pharyngeal swab data served as a positive control for the detection of sgRNA due to the high amount of replicating virus in this specimen. Swab sample data were reproduced from our own data.

26For sample preparation, paraffin-embedded liver samples were deparaffinized, resuspended in 100 μl lysis buffer (50% trifluoroethanol (TFE), 300 mM Tris-HCl, pH 8.0) and processed. Briefly, samples were sonicated at 4 °C for 15 min (‘high’ settings on a Bioruptor; Diagenode), boiled at 95 °C for 90 min and sonicated again as described above. Cleared lysates were reduced/alkylated (5 mM dithiothreitol, 25 mM chloroacetamide) for 1 h in the dark and concentrated on a vacuum centrifuge (45 min, 60 °C). A total of 50 µg proteins were digested overnight in digestion buffer (10% TFE containing trypsin/LysC (1:50 protein/protein ratio). Digestion was stopped with 1% trifluoroacetic acid (TFA) and peptides were purified on stage tips with two SDB-RPS Empore filter discs (3M) and resuspended in 2% acetonitrile (ACN)/0.1% TFA to a final concentration of 250 ng µl−1.

27Samples were analysed on a nanoElute (plug-in v.1.1.0.27; Bruker) coupled to a trapped ion mobility spectrometry quadrupole time of flight (timsTOF Pro) (Bruker) equipped with a CaptiveSpray source. Peptides (500 ng) were injected into a Trap cartridge (5 mm × 300 μm, 5 μm C18; Thermo Fisher Scientific) and next separated on a 25 cm × 75 μm analytical column, 1.6 μm C18 beads with a packed emitter tip (IonOpticks). The column temperature was maintained at 50 °C using an integrated column oven (Sonation GmbH). The column was equilibrated using 4 column volumes before loading samples in 100% buffer A (99.9% Milli-Q water, 0.1% formic acid (FA)). Samples were separated at 400 nl min−1 using a linear gradient from 2 to 17% buffer B (99.9% ACN, 0.1% FA) over 60 min before ramping up to 25% (30 min), 37% (10 min) and 95% of buffer B (10 min) and sustained for 10 min (total separation method time, 120 min). The timsTOF Pro was operated in parallel accumulation-serial fragmentation (PASEF) mode using Compass Hystar v.5.0.36.0. Settings were as follows: mass range 100–1700 m/z, 1/K0 start 0.6 V⋅s/cm2End 1.6 V⋅s/cm2; ramp time 110.1 ms; lock duty cycle to 100%; capillary voltage 1,600 V; dry gas 3 l min−1; dry temperature 180 °C. The PASEF settings were: 10 tandem mass spectrometry (MS) scans (total cycle time, 1.27 s); charge range 0–5; active exclusion for 0.4 min; scheduling target intensity 10,000; intensity threshold 2,500; collision-induced dissociation energy 42 eV.

28Raw MS data were processed with the MaxQuant software v.1.6.17 using the built-in Andromeda search engine to search against the human proteome (UniprotKB, release 2019_10) containing forward and reverse sequences concatenated with the SARS-CoV-2 polyprotein with the individual viral open reading frames manually annotated, and the label-free quantitation algorithm. Additionally, the intensity-based absolute quantification (iBAQ) algorithm and match between runs option were used. In MaxQuant, carbamidomethylation was set as fixed and methionine oxidation and N-acetylation as variable modifications. Search peptide tolerance was set at 70 p.p.m. and the main search was set at 30 p.p.m. (other settings left as default). Experiment type was set as TIMS-DDA with no modification to the default settings. Search results were filtered with a false discovery rate of 0.01 for peptide and protein identification. The Perseus software v.1.6.10.4 was used to process the data further. Protein tables were filtered to eliminate the identifications from the reverse database and common contaminants. When analysing the MS data, only proteins identified on the basis of at least one peptide and a minimum of three quantitation events in at least one experimental group were considered. The iBAQ protein intensity values were normalized against the median intensity of each sample (using only peptides with recorded intensity values across all samples and biological replicates) and log-transformed; missing values were filled by imputation with random numbers drawn from a normal distribution calculated for each sample. Differential analysis and GSEA were performed as described above on mRNA.

29The mRNA and protein datasets were integrated and the correlation analysis was performed based on two different approaches: gene-wise and sample-wise.

30Spearman’s correlation between mRNA and protein across all samples was calculated for every single gene separately. For each gene, lower- and upper-tail P values were calculated from the empirical cumulative distribution function (ECDF). Significantly highly correlated genes (upper-tail P < 0.05) were selected to perform GSEA using Fisher’s exact test on the MSigDB. A Benjamini–Hochberg-adjusted P < 0.05 was considered significant.

31For every single gene, we calculated the average fold change between positive and negative. The gene set-specific correlation was assessed by first selecting the genes that belonged to a given gene set, then by calculating the mRNA fold change versus the protein fold change Spearman’s correlation of these genes. To estimate the significance of the correlation for one gene set, correlation values were also calculated from 10,000 random sets of genes, using the same number of genes as the given gene set. Finally, the ECDF was used to return lower- and upper-tail P values.

32Functional analysis was performed with the DoRothEA package (1.4.1) for transcription factor activity and the PROGENy package (1.14.0) for pathway activity. To test the difference in activities, a factorial experiment was designed for each modality using the limma package. In both functional types, the first factor was the infection status of the patient (COVID-19+/COVID-19−) and the second was if there were traces of the virus in their liver (PCR+/PCR−). Thus, in total there were three factor combinations: a healthy patient (COVID-19−/PCR−: Ctrl); an early-stage infection (COVID-19+/PCR−: negative); and a late-stage infection (COVID-19+/PCR+: positive). A linear model was fitted for each transcription factor and pathway activity to obtain coefficients for each factorial combination. Then, these coefficients were used to compute three contrasts: effect of COVID-19+/PCR− based on COVID-19−/PCR− samples (negative versus Ctrl); effects of COVID-19+/PCR− based on COVID-19−/PCR− samples (positive versus Ctrl); and the difference between these two comparisons (positive versus negative). The first and second contrasts were aimed at detecting systematic changes in activity when an early or late-stage COVID-19 infection was taking place. However, the third contrast was aimed at detecting the specific difference in activity between late-stage and early-stage infection. The sign of the obtained coefficients can be interpreted as an increase or decrease of the mean activity for a given comparison, each with an associated probability.

33Further information on research design is available in the Nature Research Reporting Summary linked to this article.