A translational genomics approach identifies IL10RB as the top candidate gene target for COVID-19 susceptibility

Transcriptome-wide association studyTranscriptomic imputation model construction

Transcriptomic imputation models were constructed as previously described11 for peripheral tissues of the GTEx v87 (excluding brain and tissues with n ≤ 73) and STARNET8 cohorts (Supplementary Table 1; Fig. 1a). The genetic datasets of the GTEx and STARNET cohorts were uniformly processed for quality control (QC) steps before genotype imputation. We restricted our analysis to samples with European ancestry as previously described11. Genotypes were imputed using the University of Michigan server44 with the Haplotype Reference Consortium (HRC) reference panel45. Gene expression information was derived from RNA-seq gene-level counts, which were adjusted for known and hidden confounds, followed by quantile normalization. For GTEx, we used publicly available, quality-controlled gene expression datasets from the GTEx consortium (http://www.gtexportal.org/). RNA-seq data for STARNET were obtained in the form of residualized gene counts from a previously published study8. For the construction of the transcriptomic imputation models, we used elastic net-based methods; when epigenetic annotation information46 was available for a given tissue, we employed the EpiXcan11 method to maximize power; when not available, we used the PrediXcan25 method.

COVID-19 phenotypes GWAS summary statistics

Summary statistics for all 7 COVID-19 phenotypes (A1, A2, B1, B2, C1, C2, and D1; Supplementary Table 2) were obtained from the COVID-19 Host Genetics initiative4 (Release 4; 2020-10-20; https://www.covid19hg.org/results/r4/).

Multi-tissue transcriptome-wide association study (TWAS)

We performed the gene-trait association analysis as previously described11. Briefly, we applied the S-PrediXcan method26 to integrate the COVID-19 GWAS summary statistics and the transcriptomic imputation models constructed above to obtain gene-level association results.

Gene set enrichment analysis for TWAS results

To investigate whether the genes associated with a given trait exhibit enrichment for biological pathways, we used gene sets from MsigDB 5.147 and filtered out nonprotein coding genes, genes located at MHC, as well as genes whose expression could not be reliably imputed. Statistical significance was evaluated with one-sided Fisher’s exact test and the adjusted p values obtained by the Benjamini–Hochberg method48. We also performed separately LD-aware TWAS pathway enrichment analysis with JEPEGMIX2-P9 v01.1.0 with SNP and gene set annotations v0.3.0.

Genetically regulated gene expression (GReX-) based gene targeting approach

The gene targeting approach integrates genetically regulated gene expression (GReX) information (using the TWAS gene-trait-tissue association results) with a perturbagen signature library6 (Fig. 1).

Perturbagen Library used

We used the LINCS Phase II L1000 dataset (GSE70138) perturbagen reference library6; specifically the shRNA signatures (gene expression changes after knocking down a gene). All inferred genes (AIG; n = 12,328) were considered. Only “gold” signatures were considered.

Imputed transcriptomes used

We only considered GReX from 17 EpiXcan tissue models of the B2 phenotype that had significant TWAS results (FDR adjustment was applied to all COVID-19 phenotypes and tissues; Supplementary Table 1; the steps of the method are also outlined in Supplementary Fig. 5).

Signature antagonism of trait GReX

Each signature from the shRNA signature library (e.g., IL10RB shRNA treatment for 96 h in MCF7 cells) was assessed and ranked for its ability to reverse the trait-associated imputed transcriptomes using a previously published CDR method10.

Summarization of the effect of signatures across tissues

Signatures were grouped by peturbagen (shRNA), and we first tested whether the signatures for a specific perturbagen were more likely to be ranked higher or lower (Mann–Whitney U test). Then, we obtained a perturbagen-specific GReX antagonism pseudo-z-score which is defined as the negative Hodges-Lehmann estimator (of the median difference between that specific shRNA vs. the other shRNAs) divided by the standard deviation of the ranks of the shRNAs as follows: \(}\;z_}} = - \frac} - }\;}_}}}}}\;}\;}\;}\;}\;}}}\). A positive pseudo-z-score is interpreted as a potential therapeutic candidate, whereas a negative pseudo-z-score would suggest that the shRNA is not antagonizing the imputed transcriptome and is thus likely to exacerbate the phenotype. FDR is estimated with the Benjamini–Hochberg procedure48. All shRNAs were considered.

Gene prioritization approach

For prioritization, we estimated, for each gene, the p-value corresponding to the joint statistic of the two approaches \(\left( }} = \underline } + }\;z_}}} \right)\).

IL10RB and IFNAR2 GReX association with COVID-19 severity and other phenotypes in the Million Veteran ProgramCohort

Within the broader cohort of the Million Veteran Program14, for the COVID-19 severity analysis we used all COVID-19-positive individuals as of March 11, 2021 (nEUR = 14,262, nAFR = 5828, nHIS = 2870, nASN = 266; EUR: European; AFR: African; HIS: Hispanic; ASN: Asian) from MVP release 4. For the phenome-wide association study (PheWAS), we used individuals of European Ancestry (n = 296,407) from MVP release 3. Ancestries were defined by the HARE method49. Genotypes used for imputation were filtered by Minor Allele Frequency (>0.01), Variant level Missingness (<0.02), as well as imputation R2 (>0.9). We considered the MVP severity cohort an independent cohort from the GWAS since less than 7% (1519) of its participants were included in the COVID-19-related hospitalization GWAS (“B2_ALL_eur_leave_23andme”; release 4), comprising less than 7% and 0.2% of the GWAS’s cases and total individuals, respectively. For the individual imputation, we used the EpiXcan tissue model of blood from the STARNET cohort for the following reasons: (1) as tissue, blood is relevant (immune cells) and accessible—allowing for testing and validation, and (2) as an imputation model, the blood (STARNET) is the most powerful model (identifying the most FDR-significant gene-trait associations; Supplementary Fig. 4), allows the concurrent study of both IL10RB and IFNAR2 (Fig. 2a), and is based on collection from beating heart donors (in contrast to the GTEx model which is based mostly on postmortem blood). These analyses were conducted under the protocol “VA Million Veteran Program COVID-19 Science Initiative”, which was approved by the Veterans Affairs Central Institutional Review Board and by the Research & Development Committee at the James J. Peters VA Medical Center.

Phenotypes

There are four COVID-19 severity levels: mild, moderate, severe, and death (see Supplementary Table 12 for more information on the phenotypic definition and counts).

Transcriptomic imputation

GReX for blood IL10RB and IFNAR2 was calculated with the EpiXcan11 Blood (STARNET) transcriptomic imputation model (Supplementary Table 1). For TWAS, we only considered SNPs with imputation R2 ≥ 0.3. Ancestry-specific principal component analysis was performed using the EIGENSOFT50 v6 software as previously described51.

GReX association with COVID-19 severity

Associations of GReX and COVID-19 severity (Supplementary Table 12) were independently performed in each ancestry group. All associations were performed on scaled GReX with the following covariates: Elixhauser comorbidity index15 for 2 years, sex, age, age squared (age2), and top 10 ancestry-specific principal components. To estimate the effect of GReX on COVID-19-related death, we used a logistic regression analysis (binomial distribution) where death was defined as “1”, while mild, moderate, and severe cases were defined as “0”. To estimate the effect of GReX on COVID-19-related outcome severity, we performed an ordinal logistic regression where COVID-19 severity was ordered as follows: mild, moderate, severe, and death. The ancestry-specific associations were meta-analyzed with a fixed-effect model using the inverse-variance method to estimate the effect of IL10RB and IFNAR2 GReX in the total population.

GReX PheWAS

Phecodes52 assigned to clinical encounters up to 2018 (predating the COVID-19 pandemic) were grouped into categories using Phecode Map v1.2 with manual curation for some uncategorized phecodes (as provided in Supplementary Data 10). All phecodes with at least one count in more than 25 individuals in the cohort were considered for further analysis. Association testing was performed between scaled GReX and counts of each phecode with a negative binomial distribution—this is the appropriate distribution to capture the full range of phecode counts since variance was higher than the mean phecode count in 99.95% of the phecodes evaluated (1840/1841; the only exception was “Other disorders of purine and pyrimidine metabolism”). The following covariates were used: total number of phecodes per individual, length of the record, sex, age, and top 10 ancestry-specific principal components. Phecodes with nonconvergent regression models were dropped. Significance was tested at the 0.05 false discovery rate (FDR) level.

Gene expression profiling and EHR-based phenotyping in the Mount Sinai COVID-19 Biobank

Bio-specimens for this analysis were obtained from 568 individuals16; some of whom (n = 392) contributed more than one bio-specimen. The complete biobank dataset and analyses will be presented in Thompson et al.53.

RNA-seq gene expression profiling

RNA was extracted from whole blood samples and used to prepare RNA-seq libraries which were quality controlled and sequenced, as previously described53,54. In addition, we confirmed that no samples were mislabeled. RNA-seq reads were processed, quality controlled, and aligned to a reference genome as previously described53,54. After removing lowly expressed genes (keeping genes with counts per million >1 in at least half of the number of control subjects in the cohort), we normalized the raw count data of the 21,194 remaining genes using voom55 with dream56 weights from the variancePartition R package57. Samples that failed to pass all quality controls were removed from further analyses. Principal component analysis was performed using the prcomp R function to explore covariate effect on gene expression variance genome-wide. Batch effect was calculated on a per gene basis using technical replicates sequenced in all batches, as previously described53. The following additional covariates were included in the model: Subject ID, number of days since first blood sample, RNA Library Prep Plate, Sex, DV200 Percent, and PICARD metrics PCT_R2_TRANSCRIPT_STRAND_READS, PCT_INTRONIC_BASES, WIDTH_OF_95_PERCENT, MEDIAN_5PRIME_BIAS, MEDIAN_3PRIME_BIAS53. Cell type proportions were calculated with CIBERSORTx58 using the LM22 reference53,59. COVID-19 severity-associated cell types are identified as having nonzero coefficients in a linear mixed model lasso procedure (R package glmmlasso) predicting COVID-19 severity53. Finally, we used dream56 (a precision-weighted linear mixed model to consider repeated measures) for differential gene expression analysis while accounting for covariates identified by variancePartition57, above, as well as the proportions of COVID-19 severity associated with cell types identified above. A total of 6 DE signatures were generated (Severe end-organ damage Vs. Severe; Severe end-organ damage Vs. Moderate; Severe end-organ damage Vs. Control; Severe Vs. Moderate; Severe Vs. Control; Moderate Vs. Control). Multiple testing was controlled separately for each DE comparison accounting for the 21,114 genes tested using the false discovery rate (FDR) estimation method of Benjamini–Hochberg48.

COVID-19 severity scale

Phenotypic information was obtained by the EMR of the Mount Sinai Health System, which is reviewed by a screening team that includes practicing physicians. Each bio-specimen was associated, when possible—given the information in the EMR—with a COVID-19 severity measurement that corresponded to the time of collection. There are four levels of severity60: controls, moderate, severe, and severe end-organ damage, summarized in Supplementary Table 13.

Ethics statement

This study was approved by the Human Research Protection Program at the Icahn School of Medicine at Mount Sinai (STUDY-20-00341). All patients admitted to the Mount Sinai Health System were made aware of the research study by a notice included in their hospital intake packet. The notice outlined details of the specimen collection and planned research. Flyers announcing the study were also posted in the hospital, and a video was run on the in-room hospital video channel. Given the monumental hurdles of consenting sick and infectious patients in isolation rooms, the Human Research Protection Program allowed for sample collection, which occurred at the time of clinical collection and included at most an extra 5–10 cc of blood prior to obtaining research consent. Limited existing clinical data obtained from the medical record was collected and associated with the samples. As the research laboratory processing needed to begin proximal to sample collection, a portion of the data was generated prior to obtaining informed consent. During or after hospitalization, research participants and/or their legally authorized representative provided consent to the research study, including genetic profiling for research and data sharing on an individual level. In those circumstances where consent could not be obtained (13.8% of subjects, 0% of subjects who completed the post-discharge checklist), data already generated could continue to be used for analysis purposes only when not doing so would have compromised the scientific integrity of the work. The data were not identifiable to the researchers doing the analyses.

Manipulation of IL10RB and IFNAR2 expression in NGN2 cells and their effect on SARS-CoV-2 viral load and transcriptional profilesOverview

NGN2-glutamatergic neurons were derived from hiPSC-NPCs of donor NSB260761 as previously described62. gRNAs were designed using the CRISPR-ERA (http://crispr-era.stanford.edu) web tool and cloned into a lentiviral transfer vector using Gibson assembly18,62; shRNAs were ordered as glycerol stocks from Sigma. Wild type or dCas9 expressing NPCs were infected with rtTA and NGN2 lentiviruses, as well as desired shRNA or CRISPRa lentiviruses, and differentiated for 7 days before SARS-CoV-2 infection with a multiplicity of infection (MOI) of 0.1 or mock infection for 24 h. After the completion of the experiment, RNA was isolated, quality metrics were obtained, and 200 ng of RNA was processed through a total RNA library prep using the KAPA RNA Hyper Prep Kit + RiboErase HMR kit (Roche, cat no: 8098140702) following the manufacturer’s instructions. The resulting libraries were sequenced on the NovaSeq 6000 S4 flow cell (Illumina), obtaining 2 × 150 bp reads at 60 M reads per sample. An additional 20 ng of RNA was run on a SARS-CoV-2 targeted primer panel using AmpliSeq Library Plus and cDNA Synthesis for Illumina kits (Illumina, Cat no: 20019103 and 20022654). All samples were normalized, pooled, and run on the NovaSeq 6000 S4 in a 2 × 150 run targeting 750 K reads per sample. The STAR aligner v2.5.2a63 was used to align reads to the GRCh38 genome (canonical chromosomes only) and Gencode v25 annotation. The module featureCounts64 from the Subread package v1.4.3-p165 was used to quantify genes. RSeQC v2.6.166 and Picard v1.7767 were used to generate QC metrics. Differential expression analysis was performed with limma68 using the first two components of multidimensional scaling and RIN as covariates. Competitive gene set testing using sets from Gene Ontology69 and the COVID-19 Drug and Gene Set Library21 was performed with camera70. SARS-CoV-2 quantification was performed by taxonomically classifying short-read data with taxMaps71. The AmpliSeq approach confirmed the presence or absence of the virus in our samples. A more detailed version of this section can be found in the Supplementary Methods section.

hiPSC-NPC culture and donor

hiPSC-NPCs of line NSB2607 (male, 15 years old, European descent)61 were cultured in hNPC media (DMEM/F12 (Life Technologies #10565), 1× N-2 (Life Technologies #17502-048), 1× B-27-RA (Life Technologies #12587-010), 20 ng/mL FGF2 (Life Technologies)) on Matrigel (Corning, #354230). hiPSC-NPCs at the full confluence (1–1.5 × 107 cells/well of a six-well plate) were dissociated with Accutase (Innovative Cell Technologies) for 5 min, spun down (5 min × 1000 × g), resuspended, and seeded onto Matrigel-coated plates at 3–5 × 106 cells/well. Media was replaced every 24 h for 4 to 7 days until the next passage.

SARS-CoV-2 virus propagation and infections

SARS-related coronavirus 2 (SARS-CoV-2), isolate USA-WA1/2020 (NR-52281) was deposited by the Center for Disease Control and Prevention and obtained through BEI Resources, NIAID, NIH. SARS-CoV-2 was propagated in Vero E6 cells in DMEM supplemented with 2% FBS, 4.5 g/L D-glucose, 4 mM L-glutamine, 10 mM Non-Essential Amino Acids, 1 mM Sodium Pyruvate and 10 mM HEPES. Virus stock was filtered by centrifugation using Amicon Ultra-15 Centrifugal filter unit (Sigma, Cat # UFC910096) and resuspended in viral propagation media. All infections were performed with either passage 3 or 4 SARS-CoV-2. Infectious titers of SARS-CoV-2 were determined by plaque assay in Vero E6 cells in Minimum Essential Media supplemented with 4mM L-glutamine, 0.2% BSA, 10 mM HEPES and 0.12% NaHCO3, and 0.7% Oxoid agar (Cat #OXLP0028B). All SARS-CoV-2 infections were performed in the CDC/USDA-approved BSL-3 facility of the Global Health and Emerging Pathogens Institute at the Icahn School of Medicine at Mount Sinai in accordance with institutional biosafety requirements.

gRNA design and cloning and shRNAs

gRNAs were designed using the CRISPR-ERA (http://crispr-era.stanford.edu) web tool. gRNAs were selected based on their specific locations at decreasing distances from the TSS as well as their lack of predicted off-targets and E scores (http://crispr-era.stanford.edu). For lentiviral cloning: synthesized oligonucleotides were phospho-annealed (37 °C for 30 min, 95 °C for 5 min, ramped-down to 25 °C at 5 °C per min), diluted 1:100, ligated into BsmB1-digested lentiGuide-Hygro-mTagBFP2 (addgene Plasmid #99374) and transformed into NEB10-beta E. coli, according to manufacturer’s instructions (NEB # C3019H). shRNAs were ordered as glycerol stocks from Sigma (IL10RB # SHCLNG-NM_000628; IFNAR2 # SHCLNG-NM_000874). Gibson Assembly of Vectors: Unless specified, all cloning reagents were from NEB, and plasmid backbones were from Addgene (https://www.addgene.org/). Primers were synthesized by Thermo Fisher Scientific. All fragments were assembled using NEBuilder HiFi DNA Assembly Master Mix (NEB, no. E2621X). All assemblies were transformed into either DH5a Extreme Efficiency Competent Cells (Allele Biotechnology, no. ABP-CE-CC02050) or Stbl3 Chemically Competent E. coli (Thermo Fisher Scientific, no. C737303). Positive clones were confirmed by restriction digest and Sanger sequencing (GENEWIZ). The following vectors have been deposited at Addgene: lenti-EF1a- dCas9-VP64-Puro, lenti-EF1a-dCas9-VPR-Puro, lenti-EF1a-dCas9-KRAB-Puro, lentiGuide-Hygro-mTagBFP2, lentiGuide-Hygro-eGFP, lentiGuide-Hygro-dTomato, lentiGuide-Hygro-iRFP670, and pLV-TetO-hNGN2-Neo. lentiGuide-dTomato and lentiGuide-mTagBFP2-Hygro lentiGuide-Puro (Addgene, no. 52963) were digested with Mlu1 and BsiWI. dTomato was amplified from AAV-hSyn1-GCaMP6f-P2A-NLS-dTomato (Addgene, no.51085). HygroR sequence was amplified from lentiMS2-P65-HSF1_Hygro (Addgene, no. 61426). mTagBFP2 was amplified from pBAD-mTagBFP2 (Addgene, no. 3463). The P2A self-cleaving peptide sequence was amplified using a reverse primer of HygroR and a forward primer of mTagBFP2. All sequences are provided in Supplementary Table 14.

Lentiviral dCas9 effectors

To engineer a lentiviral transfer vector that expresses dCas9: VP64-T2A-Puro (EF1a-NLS-dCas9(N863)-VP64-T2A-Puro-WPRE), dCas9:VP64-T2A-Blast (EF1a-NLS-dCas9(N863)-VP64-T2A-Blast-WPRE) (Addgene, no. 61,425) was digested with BsrGI and EcoRI. T2A-PuroR was amplified from pLV-TetO-hNGN2-P2A-eGFP-T2A-Puro (Addgene, no. 79823). Fragments were then assembled using NEBuilder HiFi DNA Assembly Master Mix (NEB, no. E2621). To engineer a lentiviral transfer vector that expresses dCas9:KRAB-Puro (EF1a-NLS-dCas9(N863)-KRAB-T2A-Puro-WPRE), dCas9:VP64-T2A-Blast (EF1a-NLS-dCas9(N863)-VP64-T2A-Blast-WPRE) (Addgene, no.61425) was first digested with BamHI and BsrGI. KRAB was then amplified from pHAGE-TRE-dCas9:KRAB (Addgene, no. 50917). Fragments were assembled using NEBuilder HiFi DNA Assembly Master Mix. dCas9:KRAB-Blast was digested with BsrGI and EcoRI, and T2A-PuroR was amplified from pLV-TetO-hNGN2- P2A-eGFP-T2A-Puro (Addgene, no. 79823). Fragments were then assembled using NEBuilder HiFi DNA Assembly Master Mix. To engineer a lentiviral transfer vector that expresses dCas9:VPR-Puro (EF1a- NLS-dCas9(D10A, D839A, H840A, and N863A)-VPR-T2A-Puro- WPRE), dCas9:VPR was first amplified from SP-dCas9-VPR (Addgene, no. 63798), and T2A-PuroR was amplified from pLV-TetO-hNGN2- P2A-eGFP-T2A-Puro (Addgene, no. 79823). dCas9:KRAB-T2A-Puro was digested with BsiWI and EcoRI. Fragments were then assembled using NEBuilder HiFi DNA Assembly Master Mix.

NGN2-glutamatergic neuron induction of shRNA and CRISPRa treated neurons18,62

On day -1 NPCs were dissociated with Accutase Cell Detachment Solution for 5 min at 37 °C, counted, and seeded at a density of 5 × 105 cells/well on Matrigel-coated 24-well plates in hNPC media (DMEM/F12 (Life Technologies #10565), 1× N-2 (Life Technologies #17502-048), 1× B-27-RA, 20 ng/mL FGF2 (Life Technologies)) on Matrigel (Corning, #354230). On day 0, cells were transduced with rtTA and NGN2 lentiviruses as well as desired shRNA or CRISPRa viruses in NPC media containing 10 μM Thiazovivin and spinfected (centrifuged for 1 h at 1000 × g). On day 1, media was replaced, and doxycycline was added with 1 ug/mL working concentration. On day 2, transduced hNPCs were treated with corresponding antibiotics to the lentiviruses (1 μg/mL puromycin for shRNA, 1 mg/mL G-418 for NGN2-Neo). On day 4, the medium was switched to Brainphys neuron medium containing 1 μg/mL dox. The medium was replaced every second day until SARS-CoV-2 (MOI of 0.1) or mock infection on day 7. The samples were harvested in Trizol (Invitrogen, Cat #15596026) 24 h later. RNA was isolated by phenol/chloroform extraction prior to purification using the RNeasy Mini Kit (Qiagen, Cat # 74106).

Sequencing platform

RNA samples were submitted to the New York Genome Center and, following an initial quality check, were normalized onto two different 96-well plates for a total RNA with RiboErase assay and a SARS-CoV-2 targeted assay. For the total RNA assay, 200 ng of RNA were normalized into a plate to be run through the KAPA RNA Hyper Prep Kit + RiboErase HMR (Roche, cat no: 8098140702). This total RNA prep followed the manufacturer’s protocol with minor adjustments for automation on the PerkinElmer sciclone. Briefly, the RNA first goes through an oligo hybridization and rRNA depletion, followed by 1st and 2nd strand synthesis. The cDNA was then adenylated, and unique dual indexed adapters ligated onto the ends. Finally, samples were cleaned up, enriched, and purified. The final library was quantified by picogreen and ran on a fragment analyzer to determine the final library size. Samples were normalized, pooled and run on a NovaSeq 6000 S4 in a 2 × 150 bp run format, targeting 60 M reads per sample. For the SARS-CoV-2 targeted assay, we used the AmpliSeq Library Plus and cDNA Synthesis for Illumina kits (Illumina, Cat no: 20019103 & 20022654). Briefly, 20 ng of RNA were reverse transcribed, the cDNA targets were then amplified with the Illumina SARS-CoV-2 research panel (Illumina, 20020496). The amplicons were partially digested, and AmpliSeq CD Indexes were ligated onto the amplicons. The library was then cleaned up and amplified. After amplification, there was a final clean-up, and the libraries were quantified, pooled, and run on a NovaSeq 6000 S4, obtaining 2 × 150 bp reads.

SARS-CoV-2 quantification

Short-read data were taxonomically classified using taxMaps71. As part of the taxMaps pipeline, reads were processed prior to mapping. Adapter sequences and low-quality (Q < 20) bases were trimmed out, and low complexity reads discarded. The remaining reads were then concurrently mapped against (1) the phiX174 reference genome (NC_001422.1); (2) the SARS-CoV-2 reference genome (NC_045512.2); and (3) a combined index encompassing the entire NCBI’s nt database, RefSeq archaeal, bacterial, fungal, protozoan and viral genomes, as well as a selection of RefSeq model organism genomes, including the human GRCh38 reference72, to produce the final classification. Given that some human sequences of ancestral origin (that constitute variation between individuals) are absent from the GRCh38 reference, a small percentage of human reads usually maps to other primate genomes and, consequently, was classified as such. To obtain more accurate estimates of the human content in these samples, all reads classified as “primate” were considered of human origin and reclassified accordingly. SARS-CoV-2 viral load was determined as the number of SARS-CoV-2 reads over the host (human) reads.

Competitive gene set testing

Competitive gene set testing using sets from Gene Ontology69 and the COVID-19 Drug and Gene Set Library21 was performed with camera70. First, we performed differential expression analysis with limma68 using the first two components of multidimensional scaling and RIN as covariates to identify the signature of SARS-CoV-2 infection in our cells while adjusting for other treatments. We then performed competitive gene set enrichment analysis for all gene ontology and betacoronavirus gene sets (n = 18,553). For gene ontology datasets, we kept all significantly enriched gene sets (FDR < 0.05) and kept those with a Jaccard index less than 0.2. For the betacoronavirus gene sets, we kept all the gene sets and filtered them based on a Jaccard index of 0.2. The combined SARS-CoV-2 gene set collection with the two datasets above (significantly pruned gene ontology and all pruned betacoronavirus) was used for all the following competitive gene set testing except as otherwise indicated (n = 296). Thus, in Fig. 4d, enrichment analysis is run across the whole exploratory dataset (n = 18,553) for SARS-CoV-2 infection (first row), whereas for all other conditions, we are only exploring the combined SARS-CoV-2 gene set collection (n = 296).

Manipulation of IL10RB and IFNAR2 expression in A549-ACE2 alveolar cells and their effect on SARS-CoV-2 viral loadOverview

ACE2-expressing A549 cells (A549-ACE2), a gift from Brad Rosenberg23, were either transfected with siRNA, or transduced with the TetOne inducible system prior to infection with SARS-CoV-2. For knockdown, A549-ACE2 cells were transfected with pooled siRNAs (Dharmacon) 48 h prior to SARS-CoV-2 infection. For overexpression, A549-ACE2 cells were transduced using the TetOne Inducible Expression System (Takara Bio) with lentivirus-containing TetOne-eGFP, -IL10RB, or -IFNAR2 for 48 h followed by puromycin selection. To induce gene expression, A549-ACE2 TetOne-eGFP, -IL10RB, or -IFNAR2 cells were treated with either 0 ng/mL or 100 ng/mL doxycycline for 24 h before SARS-CoV-2 infection. Cells were either mock-infected or infected with media containing SARS-CoV-2 for a multiplicity of infection (MOI) of 0.02. The cells were incubated for 48 h, and then the cells were harvested, and the virus was inactivated by Trizol (Invitrogen) or RIPA lysis buffer for safe removal from the BSL-3 facility. Lysates were stored at −80 °C until further analysis. RNA was isolated from Trizol, cleaned using the Qiagen RNeasy Mini Kit (Cat # 74106), and quantified by QuBit. A starting input of 500 ng was used to prepare cDNA via the High-Capacity cDNA Reverse Transcription Kit (Applied Biosystems, Cat # 4368813). Real-time qPCR was performed using TaqMan probes (Thermo Fisher Scientific) for the SARS-CoV-2 S (vi07918636_s1), IL10RB (hs00175123_m1), IFNAR2 (hs01022059_m1) and control gene, GAPDH (hs02786624_g1). qPCR reactions were performed on the Applied Biosystems QuantStudio 5 and analyzed using the -ΔΔCt method for fold change expression to validate genetic manipulation and quantify SARS-CoV-2 infection. A more detailed version of this section can be found in the Supplementary Methods section.

Plasmids

The pLVX.TetOne-2xstrept-eGFP plasmid was a gift from Nevan Krogan73. pLVX.TetOne-IL10RB and pLVX.TetOne-IFNAR2 where cloned using the parental pLVX-TetOne vector (Takara Bio, Inc.) and cDNA encoding human IL10RB or h

留言 (0)

沒有登入
gif