We perfected the lentiMPRA technique by combining it with specific chromatin immunoprecipitation and with inhibition of RNA interference (ChIP-lentiMPRA). The general scheme for obtaining highly representative libraries of active promoters and enhancers is shown in Figure 1a.

For immunoprecipitation, we used antibodies to H3K27ac, which is the major chromatin tag of active promoters and enhancers. This allowed us to obtain an input (primary) library of ChIP DNA presumably enriched with active regulatory elements. Then we selected a fraction of the ChIP DNA fragments that are capable of enhancing the expression of the reporter gene. For this, we amplified ChIP DNA by PCR and cloned it to the lentiviral vector pLVPGm.1-mP, which was specially designed for lentiMPRA in this work. This vector contains the gene of CopGFP fluorescent protein under the control of the cytomegalovirus minimal promoter and the cloning site for enhancer analysis. pLVPGm.1-mP showed high efficiency in cloning and selection of genome regulatory elements. During cloning, we used inhibition of RNA interference by the Nodamura virus protein B2. Such an approach has been shown to dramatically increase the titer of lentiviral vectors [34]; however, it has not been previously used in the construction of libraries. After cloning, the pool of lentiviral particles was obtained, and MIA PaCa-2 cells were transduced by the virus library. Transduced cells were propagated on puromycin to eliminate cells without the vector, and then GFP-positive cells were selected using FACS. DNA from the secondary (output) library was isolated, sequenced by Next Generation Sequencing (NGS), and analyzed.

In the first stage, we constructed a basic second-generation self-inactivating lentiviral vector pLVPGm.1-mP (Figure 1b, the annotated sequence is in the Supplementary File pLVPGm.1-mP.gb). Similar vectors have been successfully used previously [31,35].

The pLVPGm.1-mP vector contains a region for cloning the DNA fragments under study in close proximity to the minimal cytomegalovirus promoter (mP) controlling the reporter gene of the fluorescent protein CopGFP (Evrogen, Russia). The vector also includes the mPgk-1-PuroR cassette, which provides transduced cells resistance to puromycin. On the basis of pLVPGm.1-mP, a positive control vector pLVPGm.1-EmP was also obtained, which contains a cytomegalovirus enhancer inserted in the cloning region near the mP promoter.

To test the functionality of the obtained lentiviral vectors, we used MIA PaCa-2 pancreatic ductal adenocarcinoma cells. MIA PaCa-2 cells transduced with a basic vector without an enhancer or a control vector containing a cytomegalovirus enhancer were obtained (Figure 1c). After puromycin selection, the cells were analyzed by flow cytofluorimetry. The cell population transduced with the enhancer control vector contained significantly more GFP-positive cells with higher fluorescence than the population transduced with the basic vector without enhancer (Figure S1). To show that the observed differences are indeed related to the activity of the promoter-enhancer pair rather than to the different number of viral integrations, we analyzed the lentiviral DNA content in the studied populations. For this purpose, genomic DNA was isolated from the cells and analyzed by quantitative PCR using primers targeted to lentiviral DNA and to a segment of genomic DNA (see Materials and Methods). The analysis showed that both samples contained comparable amounts of lentiviral DNA, corresponding to 2–3 copies per cell (Figure S2). This proves that the minimal cytomegalovirus promoter is indeed activated by the control enhancer in MIA PaCa-2 cells.

To evaluate the performance of the proposed approach, we analyzed viral titer using the FACS-based titration on MIA PaCa-2 cells (see Materials and Methods) and evaluated the dependence of transduction efficiency (percentage of transduced cells) and average lentiviral DNA copy number in transduced populations on virus dilutions during transduction (Figure S2).

At a titer of 2.3 million transducing units (TU) per ml and 15 mL of viral medium volume, we obtained about 35 million functional viral particles (TU). We showed that transduction with undiluted viral stock of the control vector (MOI 7.2 TU per cell), 61% of GFP-positive cells in the population were detected with a relatively low number of integrated lentiviral DNA of 2.3 copies per cell (Figure S2). These data indicate the high performance of the proposed approach, which can ensure the selection of functional sequences from millions of DNA fragments. Subsequently, all transduction experiments were performed using undiluted viral stocks.

It was previously shown that antisense transcripts in lentiviral plasmids lead to the destruction of viral pre-mRNA in the packaging cells because of RNA interference and, consequently, to a decrease in viral titer [36]. Since enhancers and promoters often initiate transcription in the opposite direction or bidirectionally [3], we hypothesized that the production of viral particles by the pLVPGm.1-mP vector containing promoters and enhancers can also be reduced by RNA interference.

To test this assumption, we examined the production of viral particles in 293T cells transduced with an equimolar mixture of pLVPGm.1 plasmids containing the CopGFP gene under the control of a strong cytomegalovirus (CMV) promoter/enhancer in sense or antisense orientations (Figure S3). Indeed, the presence of a vector containing the CMV promoter in the antisense orientation reduces the total production of viral particles by a factor of 5.6 compared with a vector containing the CMV promoter in the forward orientation alone.

Addition of a plasmid expressing the RNA interference suppressor protein NovB2 (Nodamuravirus B2 protein) during co-transfection of the packaging cells resulted in almost complete restoration of the viral titer values, which is consistent with the results of other studies [36]. Therefore, in all experiments, a plasmid with the NovB2 gene was used for the production of viral particles in order to level out the decrease in lentiviral titer caused by RNA interference.

We performed chromatin immunoprecipitation (ChIP) of MIA PaCa-2 cells with antibodies to H3K27ac to obtain a primary DNA library enriched with promoters and enhancers. The library of ChIP DNA fragments was ligated to the adapter containing the library primer annealing site and the restriction site required for cloning (See Materials and Methods). After PCR amplification, the fragments (average length about 750 bp, Figure S4a) were cloned into the basic lentiviral plasmid pLVPGm.1-mP. Approximately 20 million colonies were obtained. PCR screening of 30 randomly selected colonies showed that they all contained inserts (Figure S4b). A pool of lentiviral plasmids pLVPGm.1-L-mP containing cloned ChIP DNA fragments was isolated, producing the E. coli library of 20 million clones with close to 100% clones with inserts. This indicates a high representation of ChIP DNA fragments in the pool of pLVPGm.1-L-mP plasmids.

Using quantitative PCR, we evaluated the enrichment of immunoprecipitated DNA with fragments of functional elements at different stages of library construction using promoter regions of the housekeeping genes PSMB2, PSMB5, and COPZ1 (Figure 2).

As shown in Figure 2, the content of promoter fragments in the immunoprecipitated DNA exceeds four to 20 times the content of “non-functional” sites located at the same loci and not overlapping with any known ENCODE cis-regulatory elements [37] (see Materials and Methods).

It is important that despite the numerous procedures performed during the cloning of the primary ChIP DNA library (adapter ligation, PCR-amplification, restriction endonuclease digestion, ligation, multiple DNA purifications, etc.), its enrichment with functional sites did not change qualitatively.

A pool of lentiviral plasmids pLVPGm.1-L-mP from the previous stage was used to produce viral particles in the 293T packaging cell line. After transduction and selection of transduced cells with puromycin, 12 × 10⁵ GFP-positive cells with the highest fluorescence intensity were selected using FACS. This experiment was repeated, and 25 × 10⁵ cells were obtained additionally. Genomic DNA was isolated from both populations, and preparative PCR amplification of the selected fragments was performed (see Materials and Methods). As a result, L1 and L2 libraries, presumably enriched with functional elements, were obtained from the two populations of selected cells.

To estimate the proportion of functional elements in the obtained libraries, the activity of 28 randomly selected DNA fragments from the L1 library was analyzed in the transient expression system. The fragments were cloned into the pGRm.1-mP plasmid (Figure 3a) 5′-end to the cytomegalovirus minimal promoter (mP) in the same orientations relative to the CopGFP gene as in the lentiviral DNA. The pGRm.1-mP plasmid contains an additional reporter gene for the red fluorescent protein mKate2, located head-to-head relative to the CopGFP gene. The second reporter gene allows the detection of dual activity, such as promoter and enhancer activity, simultaneously, as well as bidirectional promoter activity in the fragments under study. Transient transfection of MIA PaCa-2 cells was performed, and the proportion of GFP-positive cells was measured. Out of 28 DNA fragments tested, six (21%) significantly enhanced the minimal promoter activity, leading to an increase in the fraction of GFP-positive cells in the transfected populations (Figure 3b). The six active elements found were sequenced and mapped in the human genome (Table 1).

All six fragments are localized in known promoter regions and overlap with known DNase I hypersensitivity sites, ChIP-seq peaks of transcription factors [39], and ChIP-seq peaks of H3K27ac in MIA PaCa-2 cells [38]. Five fragments overlap with known transcription initiation regions (CAGE-peaks) [40]. One of the fragments (#13), mapped in the promoter region of the bidirectional gene pair ATP5MK and PDCD11 (Figure 3c), enhanced the expression of both reporter genes (Figure 3d).

In our transfection experiments, just about 21% of the selected fragments contained active functional elements, and the rest seemed to be false positives. To understand the nature of the false results better, we assessed the multiplicity of lentiviral DNA integration in the L1 population using quantitative PCR (see Materials and Methods). The lentiviral DNA content in the L1 population was 2.4 copies per cell on average, which is comparable with the results obtained previously for the control populations (Figure S2). This means that the presence of a significant number of false-positive fragments in the library can be explained by multiple integrations of the lentivirus. Another source of false positives in our experiments is the so-called “position effect” that occurs when an inactive DNA fragment in the integrated lentiviral construct is influenced by a nearby genomic regulatory element (e.g., an enhancer), which activates the promoter and, consequently, the expression of CopGFP. The position effect can be overcome by non-random integration of the lentiviral construct into a certain place in the genome that is known to be not influenced by regulatory elements or isolated from them, for example, by using insulators. Unfortunately, this approach does not seem to be satisfactory because of the low efficiency of homologous recombination. This approach also ceases to be universal and requires the creation of special cell lines of each type. In principle, the position effect can be partially overcome by parallel analysis of two or more independently derived libraries through random integration of lentiviral constructs into the genome.

The DNA fragments of the obtained libraries were sequenced using massive parallel paired-end sequencing. The resulting reads were mapped to the human reference genome (GRCh38/hg38). In total, 6244 and 5245 fragments were mapped for the L1 and L2 libraries, respectively. Analysis of the overlap of the mapped fragments with different functional regions showed that both libraries are equally enriched in functional sites (DNase I hypersensitivity sites, transcription factor binding sites, etc.) compared to the library of random fragments obtained in silico (Table S1). Therefore, fragments from the two libraries were combined into one; overlapping fragments were merged, resulting in a map containing 11,264 unique regions (mapped regions) (Table S2).

The obtained map was compared with the maps of regions with elevated H3K27ac modification levels in MIA PaCa-2 and Capan-2 cells, as reported by Diafferia et al. [38]. Of the 11,264 regions we mapped, 5372 (~48%) overlap with high H3K27ac regions in MIA PaCa-2 cells (Tables S1 and S2). In Capan-2 cells of the epithelial pancreatic cancer subtype, the proportion of overlapping regions is slightly lower, ~41% (4651) (Table S2). For randomly selected genome fragments, the extent of overlapping with H3K27ac ChIP-seq peaks is significantly lower and does not exceed 5.7% (Table S1). These data confirm that the library obtained is enriched with sequences that overlap with the H3K27ac regions and contains, among others, elements specific to MIA PaCa-2 cells. For example, in the TP63 transcription factor gene locus, which is expressed in MIA PaCa-2 cells but is inactive in Capan-2 cells, 12 of 14 mapped fragments overlap with the H3K27ac regions in MIA PaCa-2 cells, 9 of which being specific MIA PaCa-2 (Figure 4a). At the same time, we found no fragments overlapping with Capan-2-specific H3K27ac regions, while only three fragments overlapped with H3K27ac regions present in both cell types.

The mapped fragments tend to co-localize with regions with the highest content of the H3K27ac chromatin tag (Figure 4b). High H3K27ac content is known to be characteristic of super-enhancers, which are clusters of enhancers with high cell specificity. Super-enhancers ensure high transcriptional activity of the genes they regulate, many of which are key cell identity genes [41,42]. Indeed, 8.7% of the mapped regions overlap with known MIA PaCa-2 super-enhancers, compared to only 0.8% of randomly selected genome fragments (random regions, Table S1). The mapped regions are also enriched in other features of functional chromatin traits (DNase I hypersensitivity sites, transcription factor binding sites, and CAGE peaks). For example, sequences containing transcription initiation regions (CAGE-peaks) [40] are seven times more frequently occurring among our mapped regions than among randomly selected regions of the genome (Table S1). A detailed analysis of the distribution of the mapped and randomly sampled regions in the genome showed that, indeed, the found regions are concentrated in the promoter regions of genes within ±2 kbp of their 5′-ends (Figure S5).

As noted above, a total of 11,264 selected fragments were mapped in the genome. Our goal was to search for functionally active elements of the genome, so we further considered only the H3K27ac enriched regions. We found that the mapped regions are concentrated near gene promoters (Figure S5). We divided them into two groups: proximal (P-regions hereafter) located within ±2 kb of the 5′-end of genes and distal (D-regions) located 2–500 kb from the 5′-end of genes (Table S2).

A total of 2537 proximal and 2801 distal regions were identified. Using the GREAT service [43], we determined the association of the mapped regions with genes. GREAT associates genomic regions with genes in whose regulatory domains it lies based on data of chromosomal conformation capture approaches as well as computational methods.

Using Expression Atlas data [44], we showed that the genes associated with proximal (P genes, 1888 genes) or distal (D genes, 2155 genes) regions (Table S4) are expressed at high levels in MIA PaCa-2 cells (Figure 5a and Table S4). This means that the regions we mapped tend to localize near the most active genes and that the libraries we obtained are enriched with elements that are probably involved in the regulation of these genes.

We divided the D genes into genes associated with 1–2 mapped regions (D1 genes, 1718 genes) and genes associated with three or more distal regions (D2 genes, 437 genes) (Table S4). Interestingly, the D2 group genes are, on average, expressed at a higher level than the P and D1 genes (Figure 5a). We also found that gene groups (P-, D1-, and D2 genes) differ from each other in the content of genes of different functional categories (Figure 5b). In particular, we distinguish genes involved in epithelial–mesenchymal transition (EMT-related genes) [46], genes significantly upregulated in pancreatic cancer [45], and genes of transcription factors showing increased expression in MIA PaCa-2 cells [47].

The vast majority of genes belong to one of the three categories reviewed, though some genes may belong to more than one category (Figure S6). As shown in Figure 5b, the proportion of genes related to these three categories is higher in gene group D2 than in gene group D1 and especially than in gene group P. At the same time, the content of the “housekeeping” genes [48] is higher among P genes. These results are consistent with the data of Zabidi et al., 2015 [49] that enhancers of “housekeeping” genes are more often concentrated near the transcription start sites (TSS) of genes than enhancers of genes with complex and tissue-specific regulation.

D2-genes of the categories “EMT-related genes” (54 genes) and “Genes significantly upregulated in pancreatic cancer” (53 genes) are, on average, expressed at higher levels in MIA PaCa-2 cells than other genes in this group (Figure 5c and Table S4). Many of the identified D2 genes encode transcription factors with the expression characteristic for MIA PaCa-2 cells, such as HMGA1, MYC, ID1, NR2F1, CEBPB, STAT3, TP63, SOX9, ZEB1, and others [47]. These results indicate that the obtained library of DNA fragments contains specific enhancer-like (regulatory, promoter, and enhancer) elements that are involved in the regulation of key identity genes that determine the individual type of MIA PaCa-2 cells.

For a more detailed analysis of the functions of the genes associated with the mapped regions, we investigated their features using the Gene Ontology database. We showed that the D2 gene group is enriched in genes that are involved in processes relevant to carcinogenesis, in particular, in the regulation of cell death (GO:0010941) and apoptosis (GO:0043065) [50], cell migration and motility (GO:0030334; GO:0040012) [51], and ERBB-linked signaling pathway (GO:1901185) [52,53] (Figure S7a). The D1 group has significantly fewer genes associated with these processes, and they are virtually absent among the P genes (Figure 6). At the same time, the group of P genes associated with proximal regions is enriched in genes involved in fundamental processes typical for housekeeping genes—regulation of DNA metabolism (GO:0051052) and chromosome organization (GO:0033044), cell cycle processes (GO:0022402; GO:0010564), responses to virus transcription (GO:0019083), DNA damage (GO:0006974) and others (Figure S7b). Thus, the gene repertory defined by Gene Ontology is consistent with the result obtained by direct analysis of gene groups (see above and Figure 5).

Approximately 1 × 10⁷ logarithmically growing MIA PaCa-2 cells were fixed, pelleted, and lysed as detailed in the methods in Supplementary Materials. The lysate was sonicated with a Cole-Parmer CP750 ultrasonic processor and precleared with Dynabeads Protein G (Life Technologies, Carlsbad, CA, USA). The lysate was incubated with anti-H3K27ac antibodies (Abcam, ab4729, Cambridge, UK), and immune complexes were sorbed on Dynabeads Protein G, sedimented, and eluted from the beads. ChIP DNA and Input DNA were extracted using a PCR purification kit (Qiagen, 28104, Hilden, Germany). Aliquots were used to assess ChIP DNA enrichment by functional elements relative to Input DNA by quantitative real-time PCR (qPCR). For this, promoter and non-promoter regions from three loci of housekeeping genes PSMB5, PSMB2, and COPZ1 were amplified. Non-promoter regions were selected from the regions that do not overlap with potential cis-regulatory elements (UCSC Browser track “ENCODE Candidate Cis-Regulatory Elements (cCREs) combined from all cell types”). Primers are listed in Table S5.

293T cells were transfected with 12 μg lentiviral plasmid, 6 μg pMD.G [66], 12 μg pCMV∆R8.91 [67], and 7.5 μg pcDNA3.1 puro Nodamura B2 (Addgene, Watertown, MA, USA, cat.#17228 [34]) using Lipofectamine 2000 (Invitrogen, Waltham, MA, USA) according to the manufacturer’s recommendations. Cells were incubated, and the virus was collected from a culture medium. Virus stocks were immediately used for MIA PaCa-2 cell transduction.

The MIA PaCa-2 cells were transduced with both non-diluted virus and with dilutions 1/32, 1/64, 1/128, and 1/256. After incubation, the viral-containing medium was replaced with a fresh medium. Then, the numbers of GFP-positive cells in each well were determined by flow cytometry, and the titer was calculated using a percentage of GFP-positive cells.

The End Repair/dA-Tailing Module (NEB, E7442S, Ipswich, MA, USA) was used to convert ChIP DNA to repaired DNA having 5’ phosphorylated dA-tailed ends. Repaired DNA fragments were ligated to adapter, obtained by annealing equimolar amounts of the two oligonucleotides (oligo 1; oligo 2; Table S5). After ligation, the DNA was purified, combined with library primer (Table S5), and amplified. An aliquot was used to assess amplified library enrichment by functional elements relative to Input DNA by qPCR.

PCR-amplified ChIP DNA was digested with Xho I restriction endonuclease and ligated into XhoI- and SalI-treated pLVPGm.1-mP. The ligated vector was used for transformation of DH5-alpha E. coli, and pool of plasmid DNAs (pLVPGm.1-mP-L) with ChIP DNA library fragments was isolated. An aliquot of the library was used to assess cloned ChIP DNA enrichment by functional elements relative to Input DNA qPCR as described above. In total, 30 colonies grown on agar plates were analyzed by PCR with primers P1 and P2 (Table S5) for the presence of inserts with primers flanking inserted DNA fragments.

We analyzed the intersection of mapped regions and random genomic regions with known ChIP-seq H3K27ac peaks obtained for MIA PaCa-2 and Capan-2 cells [38], with with “DNase I Hypersensitivity peak clusters” (95 cell types) and “Transcription Factor ChIP-seq Clusters” (340 factors, 129 cell types) from ENCODE [73], with known super-enhancers from the SEA database [74], with MIA PaCa-2 super-enhancers from SEdb 2.0 [75] and CAGE peaks from FANTOM5 project [40]. ChIP-seq peaks of H3K27ac for MIA PaCa-2 and Capan-2 lines as BED files with accessions GSM1574239 and GSM1574236, respectively, were downloaded from the NCBI GEO database [76]. BED files with ENCODE data were downloaded using the Table Browser tool of the UCSC genomic browser (http://genome.ucsc.edu (accessed on 7 April 2022)). A BED file with human 93,873 super-enhancers derived from H3K27ac ChIP-seq data for various cell types was downloaded from the SEA version 3.0 database (http://sea.edbc.org (accessed on 7 April 2022)). The coordinates of the MIA PaCa-2 super-enhancers were extracted from file SE_02_0723_SE_hg38.bed (hg38), downloaded from the SEdb 2.0 database (http://www.licpathway.net/sedb (accessed on 21 November 2022)). Coordinates of the CAGE peaks were obtained using the Table Extraction Tool of the FANTOM5 project (https://fantom.gsc.riken.jp/5/ (accessed on 7 April 2022)). Random genomic regions were generated using the ShuffleBed tool [71] on the web-based Galaxy platform (https://usegalaxy.org (accessed on 7 April 2022)) [77]. The number of random regions, their size, and chromosome distribution are in perfect accordance with those for the mapped regions. Intersections of the mapped regions and random regions with the genomic intervals of the downloaded datasets were determined using the tool “Intersect the intervals of two datasets” on the Galaxy platform.

The program SeqMonk version 1.48.1 (https://www.bioinformatics.babraham.ac.uk/projects/seqmonk/ (accessed on 7 April 2022)) was used to analyze and obtain data on the relative distribution of mapped regions around the 5-ends of genes.

GREAT [43] was used to identify associations between mapped regions and genes by applying the following rules: basal+extension: 2000 bp upstream, 2000 bp downstream, 500,000 bp max extension, curated regulatory domains excluded. To perform functional annotation of gene sets with Gene Ontology (GO), we used DAVID [54] (version 6.8, https://david.ncifcrf.gov/home.jsp (accessed on 7 April 2022). Obtained GO terms with corresponding p-adjusted values (p-adj < 0.05) were used to visualize GO terms clusters with REVIGO [78] (http://revigo.irb.hr/ (accessed on 7 April 2022)) using the following options: remove obsolete GO terms, choose species Homo sapiens (9606), SimRel (default) semantic similarity measure.

The NGS data obtained in this study can be found at https://disk.yandex.ru/d/xB-nGQbDsRoT3A (accessed on 24 November 2022) or can be obtained on request from the corresponding author.