Cells and viruses: The SARS-CoV-2 virus was isolated from a nasopharyngeal sample of a patient in Sweden and the isolated virus was confirmed as SARS-CoV-2 by sequencing (Genbank accession number MT093571). The human hepatocyte-derived cellular carcinoma cell line Huh7 (obtained from Marburg Virology Lab, Philipps-Universität Marburg, Marburg, Germany matching the STR reference profile of Huh7 [4 (
link)]), African Green monkey cell line Vero-E6 (ATCC
® CRL-1586
™) and 16HBE (human bronchial epithelial cell line, obtained from Lena Palmberg, Karolinska Institute) were used.
Antibodies and drugs: Akt (rabbit, Abcam Cat#ab8805), Akt (S473) (rabbit, Abcam Cat#ab81283), mTOR (rabbit, Abcam, Cat#ab32028), mTOR (S2448) (rabbit, Abcam Cat#ab109268), S6K (rabbit, Abcam Cat#ab32529), S6K (T389+T412) (rabbit, Abcam Cat#ab60948), eIF4EBP1(rabbit, Abcam, Cat#Ab32024), eI4EBP1 (T37) (rabbit, Abcam, Cat#ab75767), ENO-1 (rabbit, Abcam, Cat#ab155102) and HIF-1a (Clone 54) (mouse, BD Biosciences, Cat#610959), β-Actin (mouse, Sigma Aldrich, Cat#A5441) SARS-CoV-2 spike S2 (mouse, GeneTex, Cat#GTX632604) and SARS-CoV-2 nucleocapsid (rabbit, Bioserve, Cat#BSV-COV-AB-04). The drugs Wortmannin, MK-2206, Torin-1, BI-D1870, PX-478 and disulfiram was purchased from Selleckchem, US, while SAHA (vorinostat) was purchased from Sigma-Aldrich, US and Rapamycin from Abcam, US.
Infection and cytotoxicity: The infectivity dose of the virus was either determined by plaque-forming assay (for omics studies) or by determining TCID
50 in Vero-E6 cells. Infection was performed by incubating the cells with virus for one hour at 37°C, 5%CO
2 in DMEM supplemented with 5% heat-inactivated fetal bovine serum (FBS) followed by removal of virus and replenishing with fresh medium. The virus-mediated cytotoxicity was determined using Viral ToxGlo assay (Promega, US). The virus titer in the supernatant was determined by qPCR targeting either the
E-gene or
N-gene using Takara PrimeDirect probe, RT-qPCR mix (Takara Bio Inc, Japan).
SARS-CoV-2 infection of Huh7 cells for omics: Huh7 cells were plated in 6-well plates (2.5 × 10
5 cells/well) in DMEM (Thermo Fisher Scientific, US) supplemented with 10% heat-inactivated FBS (Thermo Fisher, US). At 90%–95% cell confluence the medium was removed, cells washed carefully with PBS and thereafter either cultured in medium only (uninfected control) or infected with SARS-CoV-2 at a multiplicity of infection (MOI) of 1 added in a total volume of 0.5 mL. After one hour of incubation (37°C, 5%CO
2) the inoculum was removed, cells washed with PBS and 2 mL DMEM supplemented with 5% heat-inactivated FBS was added to each well. Samples were collected at three different time points, 24, 48 and 72 hours post infection (hpi). Samples were collected for proteomics and RNAseq.
Total RNA extraction and Quantification of viral RNA: The cells (uninfected, 24hpi, 48hpi and 72hpi) were collected by adding Trizol
™ (Thermo Fisher Scientific, US) directly to the wells. RNA was extracted from SARS-CoV-2 infected and uninfected Huh7 cells and from supernatent using the Direct-zol RNA Miniprep (Zymo Research, US) and quantitative real-time polymerase chain reaction (qRT-PCR) was conducted using TaqMan Fast Virus 1-Step Master Mix (Thermofisher Scientific, US) with primers and probe specific for the SARS-CoV-2
E gene following guidelines by the World Health Organization (
https://www.who.int/docs/default-source/coronaviruse/wuhan-virus-assayv1991527e5122341d99287a1b17c111902.pdf) as described previously [5 (
link)].
Transcriptomics analysis (Illumina RNAseq): The samples were sequenced using Illumina NextSeq550 in single-end mode with read length of 75 bases. The raw sequence data were first subjected to quality check using FastQC tool kit version 0.11.8. Illumina adapter sequences and low-quality bases were removed from the raw reads using the tool Trim Galore version 0.6.1. Phred score of 30 was used as cut-off to remove low-quality bases. Quality of the data was again checked after pre-processing to assure high-quality data for further analysis. The pre-processed reads were then aligned against human reference genome version 38 Ensembl release 96. Short read aligner STAR version 2.7.3a was used for the alignment. STAR was executed by setting the parameter soloStrand to Reverse to perform strand specific alignment and rest of the required parameters were set to default. The alignment result was written in sorted by co-ordinate bam format. After the alignment gene level read count data was generated for each sample using the module featureCounts from the software subread version 2.0.0. Read counting was performed by setting attribute type in the annotation to gene_id and strand specificity to reverse. Human reference gene annotation version 38 Ensembl release 96 in gtf format was used for the read counting. Normalization factors were calculated using the R package edgeR [6 (
link)] from read counts matrix to scale the raw library sizes. Low expression genes with maximum counts per million (CPM) values under 1 per sample were removed from the sample. As recommended in RNAseq, data were transformed to CPM and variance weight was calculated using voom function. Square root of residual standard deviation against log2 CPMs was plotted to verify transformation quality.
Protein extraction and in-solution digestion: The cells (uninfected, 24hpi, 48hpi and 72hpi) were lyzed in lysis buffer (5% glycerol, 10 mM Tris, 150 mM NaCl, 10% SDS and protease inhibitor), NuPAGE
™ LDS sample buffer (ThermoFisher Scientific,US) was added and the samples were boiled at 99°C for 10 min. Aliquots of cell lysates (150 µL) were transferred to sample tubes and incubated at 37°C for 5 min at 550 rpm on a block heater and sonicated in water bath for 5 min. Each sample was reduced by adding 7 µL of 0.5 M dithiothreitol (DTT) at 37°C for 30 min and alkylated with 14 µL of 0.5 M iodoacetamide for 30 min at room temperature (RT) in the dark. Following the addition of 2 µL of concentrated phosphoric acid and 1211 µL of binding buffer, protein capturing was performed according to the manufacturer’s protocol using S-Trap™ Micro spin columns (Protifi, Huntington, NY). After washing with 150 µL of binding buffer four times the samples were subjected to proteolytic digestion using 1.2 µg trypsin (sequencing grade, Promega) for 2 h at 47°C. Then 40 µL of 50 mM TEAB was added following acidification with 40 µL of 0.2% formic acid (FA) and elution with 40 µL of 50% acetonitrile (AcN)/0.2% FA and the eluents were dried using a Vacufuge vacuum concentrator (Eppendorf, US). The resulted peptides were cleaned up in a HyperSep filter plate with bed volume of 40 µL (Thermo Fisher Scientific, Rockford, IL). Briefly, the plate was washed with 80% AcN/0.1% FA and equilibrated with 0.1% FA. Samples were filtered in the plate and washed with 0.1% FA. Peptides were eluted with 30% AcN/0.1% FA and 80% AcN/0.1% FA and dried in a vacuum concentrator prior to tandem mass tag (TMT) labeling.
TMT-Pro labeling: Dry samples were dissolved in 30 µL of 100 mM triethylammonium-bicarbonate (TEAB), pH 8, and 100 µg of TMT-Pro reagents (Thermo Scientific, US) in 15 µL of dry acetonitrile (AcN) were added. Samples were scrambled and incubated at RT at 550 rpm for two hours. The labeling reaction was stopped by adding 5 µL of 5% hydroxylamine and incubated at RT with 550 rpm for 15 min. Individual samples were combined to one analytical sample and dried ina vacuum concentrator.
High pH reversed phase LC fractionation and RPLC-MS/MS analysis: The TMTPro-labeled tryptic peptides were dissolved in 90 µL of 20 mM ammonium hydroxide and were separated on an XBridge Peptide BEH C18 column (2.1 mm inner diameter × 250 mm, 3.5 μm particle size, 300 Å pore size, Waters, Ireland) previously equilibrated with buffer A (20 mM NH
4OH) using a linear gradient of 1–23.5% buffer B (20 mM NH
4OH in AcN, pH 10.0) in 42 minutes, 23.5%–54% B in four minutes and 54%–63% B in two minutes at a flow rate of 200 µL/min. The chromatographic performance was monitored by sampling eluate with a UV detector (Ultimate 3000 UPLC, Thermo Scientific, US) monitoring at 214 nm. Fractions were collected at 30 second intervals into a 96-well plate and combined into twelve samples concatenating eight fractions representing the peak peptide elution. Each combined fraction sample (800 µL) was dried in a vacuum concentrator and the peptides were resuspended in 2% AcN/0.1% FA prior to LC-MS/MS analysis.
Approximately, 2 µg samples were injected in an Ultimate 3000 nano LC on-line coupled to an Orbitrap Fusion Lumos mass spectrometer (MS) (Thermo Scientific, San José, CA). The chromatographic separation of the peptides was achieved using a 50 cm long C18 Easy spray column (Thermo Scientific,US) at 55°C, with the following gradient: 4%–26% of solvent B (2% AcN/0.1% FA) in 120 min, 26%–95% in five minutes, and 95% of solvent B for five minutes at a flow rate of 300 nL/min. The MS acquisition method was comprised of one survey full mass spectrum ranging from
m/z 350 to 1700, acquired with a resolution of
R = 120,000 (at
m/z 200) targeting 4 × 10
5 ions and 50 ms maximum injection time (max IT), followed by data-dependent HCD fragmentations of precursor ions with a charge state 2+ to 7+ for 2 s, using 60 s dynamic exclusion. The tandem mass spectra were acquired with a resolution of
R = 50,000, targeting 5 × 10
4 ions and 86 ms max IT, setting isolation width to
m/z 1.4 and normalized collision energy to 35% setting first mass at
m/z 100.
Peptide identification and preprocessing: The raw files were imported to Proteome Discoverer v2.4 (Thermo Scientific) and searched against the
Homo sapiens SwissProt (2020_01 release with 20,595 entries) and the pre-leased SARS-CoV-2 UniProt (completed with 14 SARS-CoV-2 sequences of COVID-19 UniProtKB release 2020_04_06) protein databases with Mascot v 2.5.1 search engine (MatrixScience Ltd., UK). Parameters were chosen to allow two missed cleavage sites for trypsin while the mass tolerance of precursor and HCD fragment ions was 10 ppm and 0.05 Da, respectively. Carbamidomethylation of cysteine (+57.021 Da) was specified as a fixed modification, whereas TMTPro at peptide N-terminus and lysine, oxidation of methionine (+15.995 Da), deamidation of asparagine and glutamine were defined as variable modifications. For quantification both unique and razor peptides were requested. Protein raw data abundance was first filtered for empty rows with
in house script and quantile-normalize using R package NormalyzerDE [7 (
link)]. Principal component analysis (PCA) was applied to explore sample-to- sample relationships. One proteomics samples from the uninfected control were excluded as it turned out to be outlier.
Statistical analysis: Proteomics and transformed transcriptomics data were tested for normality using histograms with normal distribution superimposed. Differential expression through linear model was performed using R package LIMMA [8 (
link)]. LIMMA supports multifactor designed experiments in microarray, transcriptomics and proteomics. Its features are designed to support small number of arrays. The three infected replicates at 24hpi, 48hpi and 72hpi hours respectively were selected in order to perform an equi-spaced univariate time series analysis. In limma design matrix, separated coefficients were associated with time and replicates in order to extract the difference as a contrast. Moderated paired-t-test using limma with adjustment for replicates was applied. For pairwise comparisons, single factorial design was implemented to fit model with a coefficient for each of our four factors: uninfected, 24hpi, 48hpi and 72hpi. Comparisons were extracted as contrasts. In both analysis, significant differential genes and proteins were selected based on
p values after Benjamini-Hochberg (BH) adjustment. Genes with alpha value inferior to 0.05 were considered significant.
Bioinformatics Analysis: The transcriptomics and proteomics analysis were performed using all the protein-coding genes and proteins and a gene set of viral processes, response and diseases respectively. The viral response gene set is a catalogue of genes that is known to be involved in viral processes, response and diseases. The catalogue was enriched by mining biological process category of gene-ontology terms, Reactome pathways and gene sets associated with various viral diseases. Gene Ontology terms were selected by keeping, “response to virus (GO:0009615)” as parent term. All child terms of GO:0009615 were selected based on ontology term relationship “is a” and “regulates”. The pathway “Antiviral mechanism by IFN-stimulated genes” and two other events it participates were selected from Reactome database. Gene sets related to 42 virus-associated diseases and six virus-related diseases were selected from “Rare_Diseases_AutoRIF_Gene_Lists” library provided by gene set enrichment tool Enrichr [9 (
link)]. The viral response gene set contains total of 1517 protein coding genes. After filtering antiviral genes, up and downregulated proteins and transcripts were submitted separately to gene set enrichment analysis (GSEA) using gseapy v0.9.17. R package gplots v3.03 was used to generate heatmaps to display terms associated adjusted
p values contrasts over conditions.
Network and community analyses: Association analyses were performed by computing pairwise Spearman rank correlations for all features after removing null variant or genes with very low expression (RPKM < 1). Correlations were considered statistically significant at false discovery rate (FDR) < 0.01. Positive correlations were selected and used to build a weighted graph where Spearman ρ was used as weights. All network analyses were performed in igraph [10 ]. For all networks, diameter, average path lengths, clustering coefficients, and degree distributions were compared with those attained for similarly-sized random networks (Erdős-Rényi models, [11 ]). Degree centrality was computed for all networks and normalized for network size. Communities were identified by modularity maximization through the Leiden algorithm [12 (
link)]. Community centrality was computed by averaging node centrality and used to identify the most central communities in each network by degree comparison. Gene set enrichment analysis was performed on each community (
n > 30) through Enrichr for KEGG Human 2019 where backgrounds were selected based on the node number of each network. Community similarity was computed through hypergeometric testing of overlap between statistically significant KEGG terms for each transcriptomic vs proteomic pair of communities. Throughout, all statistical tests were considered at an FDR < 0.05 unless otherwise stated. All analyses were performed in Python 3.7.
Protein–protein interactions among human proteins were derived from Human Reference Interactome (HuRI). Interactions between human proteins and SARS-CoV-2 viral proteins were obtained from Human Protein Atlas (HPA). Protein interaction network is created using Cytoscape version 3.6.1 [13 ]. Edge-weighted spring-embedded layout was used for the network. R package gplots 3.03 was used to generate heatmaps to display terms associated
p values contrasts over conditions. Sankey Plot illustrates most important contribution genes to flow pathways. It was plotted using R package ggalluvial version 0.11.1 [14 ]. Scatter plots produced using ggplot2 represent the bivariate relationship between proteins and time.
Western Blot: Following 24 hpi and 48hpi with different doses of SARS-CoV-2 infection, the cells were lysed in 2x NuPage LDS sample buffer (Thermo Scientific, US) followed by boiling at 95°C for ten minutes to inactivate the virus. The protein concentration was evaluated by Pierce
™ 660nm Protein Assay kit (Thermo Scientific, US). Evaluation of protein expression was performed by running 20 μg of total protein lysate on NuPage Bis Tris 4%–12%, gels or NuPage Tris-Acetate 3%–8% gels (Invitrogen, Carlsbad, CA, USA). Proteins were transferred using iBlot dry transfer system (Invitrogen, Carlsbad, CA, USA) and blocked for one hour using 5% milk or bovine serum albumin (BSA) in 0.1% PBSt (0.1% Tween-20). Subsequent antibody incubation was performed at 4°C overnight or for one hour at room temperature for β-Actin. Membranes were washed using 0.1% PBSt and secondary antibody was incubated for one hour at room temperature using Dako Polyconal Goat Anti-Rabbit or Anti-Mouse Immunoglobulins/HRP (Aglient Technologies, Santa Clara, CA, USA). Membranes were washed using 0.1% PBSt and proteins were detected using ECL or ECL Select (GE Healthcare, Chicago, IL, USA) on ChemiDoc XRS+ System (Bio-Rad Laboratories, Hercules, CA, USA). The western blot analysis was performed by using antibodies targeting Akt, p-Akt, mTOR, p-mTOR, S6K, p-S6K, 4E-BP1, p-4E-BP1and HIF-1α. Viral RNA was quantifed from cell supernatent in all the time points as a confirmation of the infection by using Takara PrimeDirect probe, RT-qPCR mix (Takara Bio. Inc, Japan).
Drug treatment and virus infectivity: Inhibitors and modulators of PI3K and mTOR signaling pathways, namely Wortmannin, MK-2206, Torin-1, Rapamycin, BI-D1870, PX-478, Disulfiram and SAHA (vorinostat) were reconstituted in DMSO and cytotoxicity at different concentrations and time points (24h and 48 h) was determined in Huh7 cells using alamarBlue (Invitrogen, US). Changes in expression of different components of the pathway at effective and non-toxic concentration for 24 h and 48 h are shown in supplementary Figure S1. To determine the effect of these drugs on virus replication, Huh7 cells were pretreated with the drugs for twelve hours at a concentration that was suitable for 48 h incubation in DMEM supplemented with 2%FBS. The cells were infected with SARS-CoV-2 at MOI of 0.1 in presence of the drugs for one hour. The virus was removed, and the cells were further incubated for 24 h in presence of the drugs in DMEM supplemented with 2% FBS. DMSO was used as a control. The virus infectivity was determined by qPCR in both the supernatant and in the cells.
RNAScope® assay: RNAScope
® assay was performed targeting HIF-1α using the probe RNAscope
® Probe-Hs-HIF1A-C2 (ACD Bio, US) (Cat# 605221-C2) as described by us previously [15 (
link)].
Data and Code Availability. The raw RNAseq data can be obtained from the SRA using the project id. PRJNA627100. Proteomics data can be obtained from
https://zenodo.org/record/3754719#.XqgnSy2B3OQ. All the codes are available at github:
https://github.com/neogilab/COVID19
Appelberg S., Gupta S., Svensson Akusjärvi S., Ambikan A.T., Mikaeloff F., Saccon E., Végvári Á., Benfeitas R., Sperk M., Ståhlberg M., Krishnan S., Singh K., Penninger J.M., Mirazimi A, & Neogi U. (2020). Dysregulation in Akt/mTOR/HIF-1 signaling identified by proteo-transcriptomics of SARS-CoV-2 infected cells. Emerging Microbes & Infections, 9(1), 1748-1760.