Lentiviruses were produced using a standard Lipofectamine 2000™transfection protocol, using one million HEK293 cells seeded in a 25 cm² flask in DMEM medium supplemented with 10% bovine serum. A mixture of four plasmids (3 µg pMDLg/pRRE (addgene ”12,251”), 1.4 µg pRSV-Rev (addgene ”12,253”), 2 µg pVSV-G (addgene ”12,259”) and 2.5 µg of the plasmid containing the expression cassette to package the pooled CROP-seq guides) was transfected. Forty-eight h later, the medium was collected, centrifuged for 5 min at 3000 rpm, and 2.5 mL supernatant containing the viral particles was collected and used to infect cells or aliquoted and stored at −80°C. Large scale preparations of lentiviruses were produced at the Vectorology facility, PVM, Biocampus (CNRS UMS3426), Montpellier, France.

The A549 lung adenocarcinoma cell line was infected with a lentivirus produced from the plasmid lenti-dCas9-KRAB-MeCP2 (a gift from Andrea Califano, addgene 122,205) allowing the expression of a dCas9-KRAB-MeCP2 fusion protein and a gene conferring resistance to blasticidin. Infected cells were then grown in the presence of 10 μg/mL of blasticidin (Sigma). Selection of A549-KRAB-MeCP2 cells was complete within 3–5 days. Bulk blasticidin resistant cells were amplified and cloned for the CRISPRi scRNA-seq experiments. The best clone was selected according to the expression level of dCas9-KRAB-MeCP2 mRNA and to the most effective inhibition of NLUCAT1 using the NLUCAT1 sg3 RNA.

The CROPseq-Guide-Puro plasmid (Datlinger et al., 2017) (a gift from C Bock, Addgene plasmid 86,708) was digested using the restriction enzyme BsmBI (NEB R0580) for 2 h at 50°C. The relevant fragments (around 8 kB) were gel-purified using the Qiagen Gel purification kit and stored at −20°C in 20-fmol aliquots. Guides against the targeted genes (see Supplemental Table S1 for selected sequences) were cloned using the Gibson assembly method (NEBuilder HiFi DNA Assembly Master Mix, NEB E2621). Aliquoted, BsmBI-digested plasmid was mixed with 0.55 µL guide oligonucleotide (200 nM) in 10 µL total volume, combined with 10 µL 2X NEBuilder HiFi Assembling Master mix and the mixture was incubated at 50°C for 20 min. Then, 8 µL of NEBuilder Assembling mixture was incubated with 100 µL of Stable competent E coli. The mixture was heat-shocked at 42°C for 45 s and transferred to ice for 2 min. SOC medium (900 µL) was added to the Stabl2-NEBuilder mixture and the mix was incubated at 37°C for 1 h. Transformed bacterial cells (350 µL) were plated onto LB agarose plates containing ampicillin (100 μg/mL) and incubated overnight at 37°C. Individual colonies were picked and grown overnight in 5 mL of Terrific Broth medium containing 150 μg/mL ampicillin and low-endotoxin, small scale preparation of plasmid DNA were performed using the ToxOut EndoFree Plasmid Mini Kit from BioVision (K1326-250). All plasmids were verified by Sanger sequencing with the primer 5′-TTG​GGC​ACT​GAC​AAT​TCC​GT-3’.

A549-dCas9-KRAB-MeCP2 cells were infected with lentiviruses obtained from individual CROPseq-Guide-Puro plasmids, encoding individual guides. Infected cells were then grown in the presence of 1 μg/mL of puromycin (Sigma). A week later, total RNAs were purified from A549-KRAB-MeCP2 cells infected with guide encoding lentiviruses and RT-qPCR (primers sequences presented in Supplemental Table S2) were performed to measure expression of the targeted genes. RT-qPCR was performed using Fast SYBR Green Master Mix (Thermo Fisher Scientific) and ABI 7900HT real-time PCR machine. A validated guide was defined as a guide providing at least 75% inhibition of targeted gene expression compared to a control guide.

A549-dCas9-KRAB-MeCP2 cells were transduced with different amounts of the viral stock containing the library of pooled, selected gRNA. After 6 h, the virus-containing medium was replaced by fresh complete culture medium. Puromycin selection (1 μg/mL) was started at 48 h post-transduction, and 2 days later, the plate with about 30% surviving cells was selected, corresponding roughly to a MOI = 0.3. The cells were then amplified under puromycin selection for 5 days. The cells were then plated and further cultured in normoxia or in hypoxic condition (1% O2) for 3 h, 6 h or 24 h. Cells were trypsinized counted and assessed for cell viability using the Countess 3 FL (Fisher Scientific). Samples were then stained for multiplexing using cell hashing Stoeckius et al. (2018), according to the Cell Hashing Total-Seq-ATM protocol (Biolegend), using 4 distinct Hash Tag Oligonucleotides-conjugated mAbs (TotalSeq™-B0255, B0256, B0257 and B0258). Briefly, for each condition, 1.106 cells were resuspended in 100 µL of PBS, 2% BSA, 0.01% Tween and incubated with 10 µL Fc Blocking reagent for 10 min at 4°C then stained with 0.5 µg of cell hashing antibody for 20 min at 4°C. After washing with PBS, 2% BSA, 0.01% Tween, samples were counted and merged at the same proportion, spun 5 min 350 x g at 4°C and resuspended in PBS supplemented with 0.04% of bovine serum albumin at final concentration of 500 cells/μL. Samples were then adjusted to the same concentration, mixed in PBS supplemented with 0.04% of bovine serum albumin and pooled sample were immediately loaded onto10X Genomics Chromium device to perform the single cell capture.

After single-cell capture on the 10X Genomics Chromium device (3′ V3), libraries were prepared as recommended, following the Chromium Next GEM Single Cell 3’ Reagent Feature Barcoding V3.1 kit (10X Genomics) and a targeted gRNA amplification Hill et al. (2018) with respectively 6, 8 and 10 PCR cycles. Libraries were then quantified, pooled (80% RNA libraries, 10% gRNA libraries and 10% hashing libraries) and sequenced on an Illumina NextSeq 2,000. Alignment of reads from the single cell RNA-seq library and unique molecular identifiers (UMI) counting, as well as oligonucleotides tags (HTOs) counting, were performed with 10X Genomics Cell Ranger tool (v6.0.2). Reads of the gRNA library were counted with CITE-seq-Count (v1.4.2). Cells without gRNA counts were discared. Counts matrices of total UMI, HTOs, and gRNA were thus integrated on a single object using Seurat R package (v4.1.0), from which the data were processed for analysis. On the total of 19,663 cells, 817 cells without gRNA counts were discared. HTOs and gRNA were demultiplexed with HTODemux() and MULTIseqDemux(autoThresh = TRUE) functions respectively, in order to assign treatment and received gRNA for each cell. On the remaining 18,846 cells, only cells identified as “Singlet” after demultiplexing of HTO counts were conserved (14,276 cells). The repartition of cells assigned as ”Doublet” (high expression of at list 2 different gRNA), ”Negative” (no detected gRNA) and ”Singlet” (a unique detected gRNA) in all conditions is showed in (Supplemental Table S3). Finally, after transforming the data of the subset of ”Singlet” cells using SCTransform(), computing PCA, and performing KNN clustering, 2 clusters of low UMI content and high mitochondrial content cells (3087 cells) were eliminated for the rest of the analysis.

In order to gain new insights into the molecular functions of 6 hypoxia-regulated lncRNAs in LUAD cells, we performed a single-cell CRISPRi transcriptome screening based on the CROP-Seq approach (Figure 2A). We transduced A549 cells expressing double transcriptionnal repressor dCas9-KRAB-MeCP2 with a mini-library containing 12 validated gRNA targeting CYTOR (also known as LINC00152), LUCAT1, MALAT1, NEAT1, SNHG12 and SNHG21 as well as the two key regulators of the hypoxic response, HIF1A and HIF2/EPAS1 (Table 1). Each guide was individually validated by qPCR in A549 cells, showing a 75%–95% inhibition of the target compared with control cells (Figure 2B). Two additional guides, with no effect on the genome, were used as negative controls. In order to mimic the hypoxic environment in which tumors develop in vivo, we equally divided the transduced dCas9-KRAB-MeCP2 A549 cells in 4 samples that we then cultured in normoxia or in hypoxia during 3, 6 or 24 h. Cells from each sample were labeled with a specific barcoded antibody (HTOs), pooled, and simultaneously sequenced using droplet based scRNA-seq (10X Genomics Chromium). The received gRNA and the culture condition were subsequently assigned for each cell by demultiplexing both gRNA and HTOs counts respectively.

Overall, we found a balanced representation for each treatment and for each gRNA among the sequenced cells, except for the cells targeted by ”SNHG12-sg3” which were depleted in all conditions (Figure 2C, Supplemental Figure SA and Supplemental Table S3). Moreover, the expression of this particular gRNA was lowly detected in those cells, confirming previous observations that this gRNA induced cell death and that only cells with low expression survive. Inhibition of target gene expression in the presence of their corresponding gRNA were validated in all 4 conditions, as well as their progressive increase (CYTOR, LUCAT1, NEAT1, SNHG12) or decrease (HIF1A and SNHG21) during hypoxia exposure (Figure 2D).

In order to validate the SSAE approach, we first evaluated the effect of feature and cell selection on the accuracy of the model, and then the biological relevance of the detected transcriptomic perturbations induced by the knock-down of the two main regulators of the hypoxic response, HIF1A and HIF2.

We then applied the SSAE method to classify cells treated with the 6 gRNA targeting hypoxia-regulated lncRNAs and cultured in the 4 conditions. Globally, the SSAE was able to classify perturbed and control cells with a good overall accuracy around 80%, except for SNHG12 and SNHG21 (Table 5; Figure 4). However, despite their promising accuracies, we did not detect any other stable perturbations than the target gene for both MALAT1 and NEAT1 targeted cells, as indicated by the obtained high means and standard deviations of the computed ranks, while those two lncRNAs were previously associated with various gene regulation functions Dong et al. (2018); Amodio et al. (2018) (Figures 6A,B). The SSAE outcomes were different for the classification of LUCAT1-targeted cells. Indeed, the transcriptomic inhibition of LUCAT1 resulted in a stable upregulation for PCOLCE2 and ISCA1 in normoxia, HDHD2 after 3 h, and 6 h of hypoxia (Figure 6C). ISCA1 and HDHD2 encode for metal ion binding proteins, whereas TFCP2 is a known oncogene. After 24 h of hypoxia, a completely different perturbation signature was found, with at least 6 stably modulated genes, including the upregulation of KDM5C, TMEM175 and NIT1, as well as the downregulation of ATP6AP1, PEX1 and PHF20. ATP6AP1 and PEX1 are respectively components of the V-ATPase and the peroxisomal ATPase complexes, while TMEM175 is a proton channel also involved in pH regulation. KDMC5 and PHF20 are both involved in chromatin remodeling and transcriptomic regulation, while NIT1 is associated to tumor suppressor functions. This particular signature allowed the classification of 90,4% of targeted cells with an accuracy of 85,67% (Table 5; Figure 4). These results indicate that LUCAT1 inhibition may induce hypoxic condition-dependent transcriptomic modulations that potentially impact tumor survival and gene regulatory processes during prolonged exposure to hypoxic conditions, completing our previous observations Moreno Leon et al. (2019).

For LINC00152, the combined inhibition of CYTOR/LINC00152 with MIR4435-2HG (Figure 7A), whose sequences are highly homologous (99% in the 220 bp region including the most efficient gRNA), was sufficient to select half of the targeted cells with an accuracy above 85% regardless of the condition (Table 4; Figure 4).

For the SNHG12 and SNHG21 datasets, the first round of SAE selected only around 10% of perturbed cells (Table 5; Figure 4). Thus we could not run the SAE for the second round because of a too low number of cells for the 4 fold cross validation. For those 2 genes, we just reported the average accuracies obtained after the first round of the SSAE (Table 5). As SNHG21 expression is relatively low in LUAD cells and is decreased by hypoxic stress, the extent of its inhibition was therefore weaker and not sufficient to distinguish targeted from control cells. Combined with the lack of transcriptomic effect induced by its knock-down, it explains the poor classification results and the randomness of features selected for cells targeted by this particular gene under all conditions (Figure 7B).

Looking at the SSAE classification outcomes for SNHG12-targeted cells, only about 15% of them were classified as perturbed in normoxia and after 3 h of hypoxia, with a poor accuracy (Table 5; Figure 4). The number of selected cells was even worse for a longer exposure to hypoxia.

Nevertheless, the ranked list of top discriminant features between the few perturbed cells and control cells obtained for the first two time points showed a notable perturbation signature. In normoxia, it was only composed of BAG1 upregulation, whereas after 3 h of hypoxia exposure, this signature was completed by GAS5 (snoRNAs-containing lncRNA gene), ARRB2, ATF5, and ETHE1 upregulations (Figure 7C). These 5 genes are all known anti-apoptotic factors. We hypothesized that this anti-apoptotic signature was expressed by a subset of LUAD cells that were actively escaping the cytotoxic effect we systematically observed for the most efficient of the two gRNA selected for targeting SNHG12, SNHG12-sg3 (Table 1). Indeed, most of the cells classified as perturbed were specific to this particular gRNA (Figure 4). As this signature progressively attenuated over time under hypoxic conditions, we speculate that the activation of this anti-apoptotic response may be inhibited by hypoxic stress, or that hypoxia may protect against the cytotoxic effect of this guide. These results demonstrate the precision of the SAE-based approach to detect a short signature, even restricted to a small subset of cells.

We compared the classification performance and the biological relevance of extracted features between the SSAE, with and without cell selection of the most responsive cells, SAE Le et al. (2018) and Random Forests using 400 estimators and the Gini importance (GI) for feature ranking. We performed this comparison for 2 representative datasets, namely, HIF2-targeted cells versus control cells and LUCAT1-targeted cells versus control cells following 24 h of hypoxia exposure.

For the first dataset, as HIF2 inhibition induced a strong perturbation signature, the 10 first selected features between all methods were highly similar, even between SSAE/SAE and Random Forest (Figure 8A). However, using the SSAE with cell selection outperforms Random Forest with an increase of 25.25% of accuracy.

For LUCAT1 dataset, few overlaps were found between the first 10 selected features with each method. The inhibition of the target gene LUCAT1 was the only feature commonly detected (Figure 8B). ATP6AP1 and PEX1 were the only two overlapping genes between SSAE and SAE obtained signatures, while KDMC5, NIT1, PHF20 and CCDC142 were specific to SSAE regardless of the cell selection step. While the signature of Random Forests was very specific, note that the author of RF proposes two measures for feature ranking, the variable importance (VI) and Gini importance (GI): Altmann et al. (2010) showed that if predictors are real with multimodal Gaussian distributions, both measures are biased. Moreover, since using the SSAE with cell selection outperforms RF by 31.07%, SAE by 27% and SSAE without cell selection by 1.67%, it is reasonable to claim that the SSAE perturbation signature is the most relevant.