Single-cell RNA sequencing uncovers the nuclear decoy lincRNA PIRAT as a regulator of systemic monocyte immunity during COVID-19

Significance SARS-CoV-2–infected patients often display characteristic changes in the production of immune mediators that trigger life-threatening courses of COVID-19. The underlying molecular mechanisms are not yet fully understood. Here, we used single-cell RNA sequencing to investigate the involvement of the emerging class of long regulatory RNA in COVID-19. Our data reveal that a previously unknown regulatory RNA in the nucleus of immune cells is altered after SARS-CoV-2 infection. The degradation of this RNA removes a natural brake on the production of critical immune mediators that can promote the development of severe COVID-19. We believe that therapeutic intervention in this nuclear RNA circuit could counteract the overproduction of disease-causing immune mediators and protect against severe COVID-19.


Fig. S1. Identification of myeloid and lymphoid lincRNAs signatures. A) mRNA and lincRNA expression levels (median + inner quartiles) in indicated tissues (Human Bodymap). B) Consensus
Path DB KEGG pathway analysis (top left) and induced network analysis (lower left) of mRNAs from A. Pathway size = circle size. Candidates contained: color-coded (light to dark red). Overlap between pathways = line thickness. Overlap of genes contained in the top 5 pathways is illustrated in the right panel (each grey cell stands for one gene). Genes annotated in at least 3 signaling pathways are indicated. C) and D) PCA analysis with leukocyte-enriched lincRNAs and clustering of mRNAs from A (cell type markers indicated). E) Overlay histogram FACS plots showing successful enrichment of the indicated cell populations by MACS (left peak: unstained control; right peak: cells stained for indicated surface antigen). F) RNA-Seq based row Z-scores of selected myeloid (top 3) and lymphoid (bottom 3) lincRNA markers (data from Fig. 1B). G) and H) qRT-PCR validation of lincRNAs from indicated purified cell types, relative to human brain reference tissue. Horizontal bar indicates base-line (black) and 2-fold deviation from base-line (grey). Box plots and individual replicate values from four independent experiments are shown. I) Summary of lincRNA expression patterns in the studied leukocyte populations.       , Table S7). Fold-changes relative to unstimulated control cells. B) Ruxolitinib and BAY-11-7082 sensitivity of CXCL2, CXCR4 and LUCAT1 (monocytes; PAMP = 4 h LPS + polyI:C; inhibitor pre-stimulation: 2 h). Fold-changes relative to unstimulated vehicle control. C) Rescue of CXCL2, NAMPT and CXCR4 dysregulation in LUCAT1 deficient THP1 cells (cell line from Fig. 3) upon 2 h Ruxolitinib treatment. A-C: One-way ANOVA, 3 independent experiments. * = P ≤ 0.05; ** P ≤ 0.01; n.s. = no significant difference.     Fig. 5A. B) Consensus Path DB KEGG and Reactome pathway analysis of genes regulated in COVID-19 and upon PIRAT expression-manipulation in THP1 cells (overlap of Venn diagram in Fig. 5A). Pathway terms, pertaining genes, pathway sources and p-values are shown. Grey fill color indicates that a given gene is included in the respective pathway. C) Same as B, but for genes regulated in COVID-19 only (Fig. 5A, left circle, without overlap). D) Same as B, but for genes regulated upon PIRAT expression manipulation only (Fig. 5A, right circle, without overlap). E) qRT-PCR analysis of CHI3L1 expression in wild-type and PIRAT-deficient (+/-and -/-), as well as in PIRAT overexpressing (OE) THP1 cells. Mean, individual replicate values and standard deviation based on three independent experiments are shown. One-way ANOVA. F) Representative dot plots illustrating the gating strategy used for FACS-quantification of CD11c (ITGAX) positive THP1 cells (Fig. 5F). G) qRT-PCR analysis of expression changes of S100A8 and S100A9 upon PU.1compared to control-CRISPRi (THP1 cells. ≥ 5 experimental replicates. Replicate values, mean and standard deviation are shown. Two-tailed Student's t-test.). Predictions were done using the ENCODE and ChEA Consensus TFs tab at Enrichr (https://maayanlab.cloud/Enrichr/). Top 10 transcription factors are shown (for the PIRAT controlled gene set only 4 transcription factors were identified). Transcription factor predictions supported by experimental data in the present study are highlighted. B) qRT-PCR analysis of the expression of S100A8 and S100A9 in control (C) and LUCAT1 knockdown (KD) THP1 monocytes, stimulated with LPS + polyI:C (PAMP) for 4 h or left untreated. Results obtained with second knockdown cell line, generated with guideRNA 2 (gRNA2, Table S7). Fold-changes relative to unstimulated control cells. C) Restoration of S100A8/9 expression in LUCAT1 KD cells upon 2 h Ruxolitinib and 4 h LPS + polyI:C treatment. B-C: Three independent replicates, One-way ANOVA. * = P ≤ 0.05; ** = P ≤ 0.01.    Fig. 6G and H, but with additional primer pairs. B) Enumeration of single nucleotide variants discriminating cloned and Sanger-sequenced PCR amplicons of the PIRAT binding sites amplified using the same primers as in panel A (elution fraction, "LINC"), but with Advantage 2 proof-reading polymerase. % indicates percentage of variant among all Sanger sequences. # indicates the total number a given variant was detected in experiments performed with ChIRP DNA from monocytes isolated from blood samples of three different donors (D1, D2, D3). C) Aligned Sanger sequences from analysis performed in B, with enumerated nucleotide variants highlighted. D) ENSEMBL Genome Browser view of all obtained BLAST hits for Sanger sequences from C (human GRCh38 genome).  Table S1. Myeloid and lymphoid cell specific lincRNAs identified in the current study (RPKM and standard-deviation are shown).

Cell culture and human biomaterial
Human peripheral blood mononuclear cells from healthy donors (control patient cells see below) were isolated from buffy coats (transfusion medicine department, UKGM Giessen). Buffy coats were de-identified prior to further use. Leukocyte populations were purified from buffy coats using Lymphoprep gradient medium (Stemcell Technologies)  Patients suffering from SARS-CoV-2-infection were recruited after hospitalization. In addition, healthy subjects were recruited (Table S2 and  Obtained cells were analysed immediately. Further BALF (Fig. 7F) was obtained from patients at the Department of Infectious Diseases and Respiratory Medicine, Charité, Berlin. All patients underwent bronchoscopy including BAL on clinical indication and had provided oral and written informed consent. The study was approved by the local ethics committee (EA2/086/16). BAL was performed by instillation of 150 ml pre-warmed sterile 0.9% NaCl solution. In patients with focal abnormalities in chest imaging, BAL was performed in the corresponding pulmonary segment; in patients without radiological abnormalities or diffuse infiltrates, BAL was performed in the right middle lobe or lingula. Diagnosis of infection was made by a board-certified pulmonologist based on chest imaging, clinical signs of infection, culture and laboratory results, BALF cellular analysis and response to therapy. For the infection group, patients with non-mycobacterial infection were selected. Control patients showed no apparent lung disease and underwent bronchoscopy and BAL as part of rule-out diagnostics due to idiopathic coughing, for exclusion of pulmonary involvement of systemic disease (e.g. rheumatoid arthritis) or for exclusion of pulmonary tuberculosis. No obvious abnormalities in chest imaging and BALF composition were detected in these patients. Patient characteristics are listed in Table S6. All studies and procedures to obtain human specimen were conducted according to the Declaration of Helsinki.

qRT-PCR
RNA was isolated using TRIzol (Ambion), treated with DNaseI (Thermo Fisher) in the presence of recombinant RNase inhibitor (Promega) and concentration was determined (Nanodrop 2000 spectrometer, Thermo Scientific). cDNA was generated (High-Capacity cDNA Reverse Transcription Kit, Thermo Fisher) and quantitative PCR was performed (PowerUP SYBR Green Master Mix, Thermo Fisher) using a QuantStudio 3 instrument. For subcellular fractionation and CoIP analysis the Power SYBR RNA-to-Ct 1-Step Kit (Thermo Fisher) was used. Relative expression was calculated based on CT values, using the 2 -ΔΔCT method (3), where applicable relative to U6 snRNA. Primers are listed in Table S7.

Subcloning and sequencing of REXO1LP amplicons
DNA from the elution fractions of PIRAT ChIRP experiments was used as a template for PCR reactions with the same primers as in Fig. S15A and using Advantage 2 proof-reading PCR polymerase (Takara). Experiments were conducted with ChIRP DNA from monocytes from three different blood donors. PCR products were subjected to gel-purification and sub-cloned using the Strataclone TA PCR cloning kit (Agilent). Insert sequences were determined by Sanger sequencing (Seqlab GmbH) and aligned using Multalin (http://multalin.toulouse.inra.fr/multalin/). Nucleotide variants in the sequenced inserts were counted and assigned to the three different donors.

RACE-PCR
RACE-PCR was performed using the SMARTer 5'/3' RACE kit (Clontech) according to the manufacturer's instructions. Template poly(A) RNA was purified using oligo-d(T) coupled dynabeads (Thermo Fisher). RACE-PCR primers are listed in Table S7. RACE products were subjected to gel-purification and sub-cloned using the Strataclone UA PCR cloning kit (Agilent). Insert sequences were determined by Sanger sequencing (Seqlab GmbH).

Copy number enumeration
PIRAT cDNA (Fig. S8) was amplified using Phusion PCR polymerase (Thermo Fisher), according to the manufacturer's instructions, and primers listed in Table S7. The forward primer (OBS-1898) contained a T7 RNA polymerase consensus binding site. The PCR amplicon was extracted from an agarose gel (Macherey-Nagel™ NucleoSpin™ Gel and PCR Clean-up Kit) and 150 ng were used as a template for RNA in vitro transcription (for 4 h), using the MEGAscript T7 transcription kit (Thermo Fisher), according to the manufacturer's instructions. Synthesized RNA was cleaned up using the Monarch RNA Cleanup column kit (NEB). RNA concentration was determined using the Agilent Bioanalyzer system, with an RNA Nano chip. RNA integrity was additionally controlled by running a sample on a MOPS / 1.2 % agarose / 1 % formaldehyde gel (10x MOPS buffer: 50 mM MOPS, 50 mM Na-acetate, 10 mM EDTA, pH 7.0), with Millennium RNA size marker (Thermo Fisher). PIRAT copy number was determined by qRT-PCR, as described in the main manuscript text.

Subcellular fractionation
Cells were lysed (10 mM Tris, pH 8, 140 mM NaCl, 1.5 mM MgCl2, 0.5 % Igepal, 2 mM vanadyl ribonucleoside complex), incubated on ice for 5 min and centrifuged (1000 x g, 4 °C, 3 min). The supernatant (cytosolic fraction) was transferred to a new tube, centrifuged (3 min, maximum speed) and transferred to a new tube for RNA-extraction. The pellet (nuclear fraction) was washed two times with lysis buffer and once with lysis buffer containing 0.5 % deoxycholic acid (centrifugations at 4 °C and 1000 x g), followed by RNA-extraction.

RNA-FISH
Tissues were derived as described (4) Table S8. Western blot full-scan is shown in Fig. S16A.

Flow cytometry
Cells were identified by plotting the respective fluorescence channel against backgroundfluorescence or the side-scatter. The gating strategy is illustrated in Fig. S11F. For surface marker staining, 2 µl of fluorophore-coupled primary antibody were added to cells in 100 µl PBS containing 1 % FCS, followed by incubation on ice for 30 min. Cells were washed and resuspended in PBS containing 0.5 % FCS and subjected to FACS analysis (Guava EasyCyte, Millipore).

UV crosslinking & Co-immunoprecipitation
For co-immunoprecipitation (CoIP), 10 7 cells were UV-crosslinked (300 mJ / cm 2 ) in petri dishes, on an ice bath. The CoIP procedure published by Tawk et al. (5) was used with minor modifications. For protein purification protein G dynabeads (Thermo Fisher), coupled with 2.5 µg of antibody (Table S8)  The cDNA reaction mix was prepared as indicated in the manufacturer's protocol and mixed with the beads. The mixture was incubated in a thermomixer (37 °C, 1200 rpm, 20 minutes). The supernatant was removed and replaced by the Exonuclease I mix prepared according to the manufacturer's protocol. The bead suspension was placed on the thermomixer (37°C, 1200 rpm, 30 minutes, followed by 80 °C without shaking for 20 minutes). The suspension was then briefly placed on ice and the supernatant was removed. Finally, beads were resuspended in Bead Resuspension Buffer (Cat. No. 650000066).
Single cell mRNA, multiplex sample Tag, and AbSeq libraries were prepared using the BD Rhapsody™ Single-Cell Analysis system (Cat. No. 633774) according to manufacturer's recommendation (Doc ID: 214508). Briefly, the Bead Resuspension Buffer was removed from the beads and replaced by the PCR1 reaction mix containing the primers specific for the AbSeq and mutiplex sample tags, and genes of the Human Immune Response Panel (Cat. No. 633750) supplemented with custom-made primers for additional genes (see NCBI GEO GSE142503). Beads were placed in the thermal cycler for 11 cycles of the PCR program indicated in the protocol. The supernatant was retained and the PCR products for Abseq and multiplex sample tags, as well as the mRNA PCR product were separated and purified by double-sided size selection using AMPure XP magnetic beads (Cat. No. A63880; Beckman Coulter, Krefeld, Germany). A fraction of the Abseq/multiplex sample tag PCR 1 product, as well as the mRNA PCR 1 product were further amplified with a second PCR of 10 cycles and subsequent purification using AMPure XP beads, resulting in multiplex sample tag and mRNA PCR 2 product. Finally, the Abseq/multiplex sample tag PCR 1 product for the Abseq library, and both PCR 2 products for each the multiplex sample tag and mRNA libraries were amplified by the final index PCR for 7 cycles each with subsequent purification afterwards. Concentrations of the index PCR products were determined using the Qubit Fluorometer and the Qubit dsDNA HS Kit (Cat. No. Q32851; Thermo Fisher Scientific) and quality control was performed on the Agilent 2100 Bioanalyzer with the High Sensitivity DNA Kit (Cat. No. 5067-4626; Agilent, Waldbronn, Germany). Mixed libraries were sequenced on a NextSeq550 with 2 x 75 bp paired-end reads.
After pre-processing of BD Rhapsody scRNA-seq data, read counts were loaded into the R (v3.6.3) environment and further analyzed using the Seurat package (v3.1.4). The following quality criteria were used to include cells for the downstream analysis: at least 25 genes were expressed, and at least 1,000 but no more than 70,000 transcripts were detected per cell.
Following the Seurat workflow, the read counts were normalized and scaled by NormalizeData and ScaleData functions of the Seurat package, respectively. Principal component analysis (PCA) was performed by RunPCA using top 2,000 variable features that were selected using the default selection method ("vst") in Seurat. Next, based on the first 15 PCs, cell clusters were identified with the Louvain algorithm at resolution of 0.4. Finally, in a two-dimensional space, a UMAP was generated to visualize the identified cell clusters.
To identify marker genes of each cell cluster, differentially expressed (DE) genes were tested by FindAllMarkers functions in Seurat using the default test (Wilcoxon Rank Sum test). Significantly differentially expressed genes were determined by 1) log-fold changes > 0.3, 2) expressed in at least a fraction of 0.2 cells in each tested population, and 3) adjusted p value < 0.05 (Bonferroni correction). DE analyses were used to identify cluster marker genes by comparing the expression of upregulated genes in cells between one cluster and the rest of cells.
Cell clusters were firstly assigned using the SingleR (v1.0.6) package based on four reference dataset which are provided in the package, including BlueprintEncodeData, DatabaseImmuneCellExpressionData, HumanPrimaryCellAtlasData, and MonacoImmuneData. Then, the assigned cell cluster annotations were double-checked by comparing the cell type specifically expressed marker genes from public resources. To dissect the different profiles between COVID-19 patients and controls, publicly reported COVID-19 related genes were selected and their expression profiles in patients and controls were