The genotype of rs9390123 surrounding 1-Mb segment (chr6: 143443314-144443314 relative to human genome build 37) was downloaded from the 1000 Genomes Project for three representative populations in the world, namely CEU (Utah Residents with European Ancestry), CHB (Han Chinese in Beijing) and YRI (Yoruba in Ibadan), by using an online program, VCF (Variant Call Format) to PED (pedigree) converter (http://grch37.ensembl.org/Homo_sapiens/Tools/VcftoPed). The obtained genotype data file (PED format) and locus information file (INFO format) were loaded into Haploview 4.2 (Broad Institute) (7) and r² was calculated with default parameter. When r²≥0.8, the SNPs were supposed to be in strong LD.

The human lung epithelial cell line Beas-2B (cat. no. KCB200922YJ), lung cancer cell line A549 (cat# KCB200434YJ) and 293TN cells (cat. no. KCB2013068YJ) were maintained in high-glucose DMEM (HyClone; Cytiva) supplemented with 10% fetal bovine serum (FBS, Biological Industries; Sartorius AG) in 5% CO₂ at 37°C. All cell lines were purchased from Kunming Cell Bank of Type Culture Collection of Chinese Academy of Sciences (http://www.kmcellbank.com/). The cell bank performed short tandem repeat genotyping, karyotype and isoenzyme analyses to verify the cell identity in order to avoid potential contamination.

To include all the potential cis-regulatory elements and amplify the segment efficiently, rs9390123 and its nearby region (~1.5 kb; chr6: 143942517-143944088) were amplified via PCR with the forward primer 5′-GGCAAATCCTTCCCATAGTTCC-3′ and the reverse one 5′-TAGGGCCAAATTTCACAGGTCTTA −3′ containing MluI and BglII restriction sites, respectively. The DNA isolated from Beas-2B cells was utilized as a template. The thermocycling conditions are as follows: 98°C for 30 sec; 35 cycles of 98°C for 10 sec, 65°C for 30 sec, 72°C for 30 sec; 72°C for 2 min. To avoid potential PCR errors, amplification was performed by utilizing Q5 High-Fidelity DNA Polymerase (New England BioLabs, Inc.). Upon digestion with the aforementioned restriction enzymes (New England BioLabs, Inc.), the PCR product was inserted into the cloning site of the pGL3-promoter plasmid (Promega Corp.). The plasmid containing another allele was produced for each SNP by using a Q5 Site-Directed Mutagenesis kit (New England BioLabs, Inc.) and the primers listed in Table SI.

After seeding (~10⁵ cells) into a 24-well plate and culture for 24 h as described above, 475 ng plasmid with the aforementioned inserted segment was transfected into the Beas-2B cells by using Lipofectamine^® 2000 (Thermo Fisher Scientific, Inc.). After additional 36 h of culture and cell harvest, the luciferase expression was determined by using a Dual Luciferase Reporter Assay System (Promega Corp.). A total of 25 ng pRL-TK plasmid (Promega Corp.) was transfected simultaneously as an internal control to normalize the transfection efficiency, and the luciferase ratio between firefly and Renilla was utilized to represent the activity of the enhancer containing rs9390123 and rs9399451. In total, six replicates were performed for each plasmid transfection.

The spatial conformation of the enhancer and the related gene promoter was examined by using 3C as previously described (8). Briefly, after subjecting Beas-2B cells (~10⁸) to crosslink by using formaldehyde and subsequently lysing the cells, the chromatin was directly digested with BsrGI (New England BioLabs, Inc.) and further ligated with T4 DNA ligase (New England BioLabs, Inc.). The ligation product was purified by using the standard phenol-chloroform method (9). Simultaneously, a BAC (bacterial artificial chromosome) RP11-1012I24 (BACPAC Genomics; http://bacpacresources.org/) harboring the 6q24.2 segment was grown in LB (Luria-Bertani) medium and purified by Large-Construct Kit (Qiagen Corp.) according to the manufacturer's protocol. After digestion, ligation and purification as aforementioned, the product was utilized as a control to normalize the primer efficiency. The relative enrichment of the 3C product was assessed by quantitative PCR (qPCR) with the unidirectional primers shown in Table SII and iTaq Universal SYBR Green Supermix (Biorad Corp.). The thermocycling conditions are as follows: 95°C for 10 min; 40 cycles of 95°C for 5 sec, 60°C for 30 sec. The experiment was repeated in triplicate. All 3C PCR products were verified using sequencing with constant primer (Table SII) and BigDye^® Terminator v3.1 (Thermo Fisher Scientific, Inc.) according to the manufacturer's recommendation.

Potential transcription factors (TFs) were predicted by an online program Match (http://www.gene-regulation.com/cgi-bin/pub/programs/match/bin/match.cgi). ChIP was performed with EZ-ChIP kit (MilliporeSigma). In brief, cells (~10⁷) were crosslinked by formaldehyde, lysed by lysis buffer (MilliporeSigma) and fragmented into segments of 200–800 bp by using sonication. After diluting by dilution buffer and preclearing by Protein A beads (MilliporeSigma), the chromatin/protein complex was captured overnight at 4°C by addition of 2 µg anti-mouse POU class 2 homeobox 1 (POU2F1) antibody (Santa Cruz Biotechnology, Inc.; cat. no. sc-53830) or the same amount of normal mouse IgG (Santa Cruz Biotechnology, Inc.; cat. no. sc-2025) and further immunoprecipitated by using Protein A Beads (MilliporeSigma). After washing, resuspending, crosslink reversing and protein digestion with related buffer (MilliporeSigma), the DNA was recovered by using a column (MilliporeSigma), and quantified via qPCR as described above with the primers shown in Table SIII. The experiment was repeated in triplicate, and the PCR product was sequenced for validation by using forward PCR primer in Table SIII and BigDye^® Terminator v3.1 as described above.

The POU2F1 coding region was obtained by using nested PCR with Q5 High-Fidelity DNA Polymerase. The PCR primer sequences, annealing temperature and target region are displayed in Table SIV. After dual enzymatic cleavage with KpnI and HindIII (New England BioLabs, Inc.), the POU2F1 coding region was inserted into the pEGFP-N1 overexpression vector (Clontech Laboratories; Takara Bio USA, Inc.). A total of 500 ng POU2F1 overexpression plasmid was transfected into Beas-2B cells with Lipofectamine^® 2000 as aforementioned. pEGFP-N1 was utilized as a negative control. After 48 h of culture, the relative mRNA levels of POU2F1 and target genes were evaluated by qPCR with primers in Table SV. In total, three independent repeats were carried out.

The probes for both alleles of rs9390123 and rs9399451 (Table SVI) were labeled with biotin. Nuclear extracts were isolated from Beas-2B cells by using a Nuclear and Cytoplasmic Protein Extraction kit (Beyotime Institute of Biotechnology) and incubated with biotin-labeled probes (Sangon Biotech; 10 fmol). The probe-protein complexes were separated via electrophoresis in 4.9% non-denatured polyacrylamide gel and transferred to nylon membranes with positive charge (Beyotime Institute of Biotechnology). For each SNP allele, electrophoresis of only biotin-labeled probes and probe-protein complexes incubated with non-labeled probes (competitor oligonucleotides) was performed as a control. After incubation with 7.5 µg streptavidin-HRP (horseradish peroxidase) conjugate (Beyotime Institute of Biotechnology; cat. no. A0303) at room temperature for 15 min with 150 rpm rotation, the membrane was visualized with enhanced chemiluminescence.

RNA in the nucleus and cytoplasm of Beas-2B cells was separated and isolated by using an RNA subcellular isolation kit (Active Motif, Inc.). Reverse transcription (RT) was performed with SuperScript III First-Strand Synthesis System (Thermo Fisher Scientific, Inc.). PHACTR2-AS1 expression was determined by using qPCR with primer pairs from the literature (10–12) (Table SV). U6 and GAPDH (glyceraldehyde-3-phosphate dehydrogenase) are known to be expressed mainly in the nucleus (13) and cytoplasm (14), respectively. Therefore, their expression was also determined by qPCR with primer pairs shown in Table SV, and it was used as positive controls to verify the separation of nucleus and cytoplasm.

shRNA (5′-GCCCTGCATACTGTGGATTCA −3′) for PHACTR2-AS1 was designed with the online software BLOCK-iT RNAi Designer (https://rnaidesigner.thermofisher.com/rnaiexpress/; Thermo Fisher Scientific, Inc.). The restrictive sites for EcoRI and BamHI were added into shRNA sequence and both forward and reverse strands were synthesized (Sangon Biotech). After annealing, the two strands formed double-strand DNA with sticky end. After EcoRI and BamHI (New England BioLabs, Inc.) digestion, the pGreenPuro vector (System Biosciences, LCC) was ligated with the DNA segment containing shRNA by T4 DNA ligase. The lentiviral pGreenPuro vector (3 µg) was transfected into 293TN cells by Lipofectamine^® 2000 along with packaging (psPAX2; addgene# 12260; 2.5 µg) and envelope (pMD2.G; addgene# 12259; 2.5 µg) plasmids to produce viruses. 293TN is a genetically modified cell line to produce high titer virus. After 48 h of culture, the viruses were collected via centrifugation. The empty pGreenPuro vector (3 µg) was also transfected into 293TN cells along with packaging and envelope plasmids to produce viruses, which were further used as a negative control (sh-NC). Beas-2B and A549 cells were transfected with the aforementioned lentivirus and screened with puromycin (2 µg/ml final concentration) for 2 weeks in order to obtain a stable cell line.

Total RNA was purified from Beas-2B and A549 cells using TRIzol^® (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions, quantified by using Nanodrop 2000 (Thermo Fisher Scientific, Inc.) and RT was carried out as aforementioned. PHACTR2-AS1 expression after knockdown was evaluated via PCR or qPCR with the primer pair listed in Table SV. GAPDH was also amplified as a positive control by using the primer pair listed in Table SV.

Beas-2B and A549 cells with lentiviral transfection were cultured for 4 days in 96-well plates (6,000 cells per well). Every 24 h, 20 µl MTT (5 mg/ml; Beijing Solarbio Science & Technology Co., Ltd.) was added to the cells and incubated for 4 h at 37°C. Four time points were included in this assay. After removing the supernatant and adding 150 µl DMSO to each well, the plate was covered with foil and placed on an orbital shaker for 15 min at room temperature. Next, the optical density (OD) value was examined at a wavelength of 490 nm. A total of three independent repeats were carried out for each time point.

Beas-2B and A549 cells were adjusted to a density 500 cells per well and cultured for 14 days in a 6-well plate. After washing with PBS, the cells were fixed with paraformaldehyde and stained with 1% crystal violet (Beijing Solarbio Science & Technology Co., Ltd.) for 20 min at room temperature. The colony number was then determined visually using a DMi8 Automated Inverted Phase Contrast Fluorescence Microscope (Leica Microsystems). Three independent repeats were performed.

Beas-2B and A549 cells (10⁴) were resuspended into DMEM without FBS and placed in the upper chamber of a Transwell cell culture plate (MilliporeSigma), while DMEM with 10% FBS was placed in the lower chamber. After 2 days of culture at 37°C, paraformaldehyde fixing and crystal violet staining of Beas-2B and A549 cells were performed as aforementioned, and the number of cells was examined using a DMi8 Automated Inverted Phase Contrast Fluorescence Microscope as described above. Three independent repeats were carried out.

Beas-2B and A549 cells were adjusted to a density of 7×10⁵ cells per well and maintained to near full (95%) confluence in DMEM with 10% FBS. After scratching the cell monolayer with a pipette tip, the cells were washed with PBS for 2 min at room temperature (twice), and fresh DMEM without FBS was added. The wound areas were imaged at 0 and 48 h using an inverted microscope, and were analyzed using ImageJ software v1.51 (https://imagej.net/). Three independent repeats were performed.

The luciferase activity for each plasmid, ChIP enrichment, gene expression after POU2F1 overexpression, and cell proliferation, colony formation, cell migration and wound healing abilities were presented as mean ± standard deviation (SD). Independent (unpaired) Student's t-test was utilized to evaluate difference in the above data between two groups. All statistical analyses were performed with SPSS 20.0 (IBM Corp.). P<0.05 was considered to indicate a statistically significant difference.

To identify the potential SNP(s) associated with DRC, the genotype data surrounding rs9390123 in its nearby region were obtained from the 1000 Genomes Project for three representative populations in the world. For the 1M segment surrounding rs9390123, 3,927, 4,285 and 6,800 SNPs were found to exist in the CEU, CHB and YRI populations, respectively. However, only one SNP, rs9399451, presented high LD with rs9390123 in CEU (r²=0.979) and YRI (r²=0.851; Fig. S1), thus suggesting that this SNP may also be associated with DRC in these two populations. In the CHB population, the LD between these two SNPs was found to be moderate (r²=0.645; Fig. S1). Except rs9399451, all other SNPs within this region showed relatively low LD with rs9390123 (all r²≤0.668, 0.633 and 0.675 in CEU, CHB and YRI, respectively; Fig. S1). There was a distance of only 560 bp between these two SNPs, suggesting that they both could form part of the same functional element. For the entire human population, T is the minor allele for both rs9390123 and rs9399451, and their frequencies were found to be ~38% and 30% for CEU and YRI populations, respectively. In contrast, T frequencies in CHB were ~58 and ~68% for rs9390123 and rs9399451, respectively.

Since neither rs9390123 nor rs9399451 are within protein-coding regions, it was hypothesized that these two SNPs could influence target gene expression. To evaluate the function of rs9390123 and rs9399451 on gene expression, the segment surrounding rs9390123 and rs9399451 was inserted into the cloning site of the pGL3-promoter vector. The sequencing indicated that the nucleotide is T for both rs9390123 and rs9399451 in the plasmid construct (Fig. 1). The plasmids with corresponding allele, i.e., C at rs9390123 and G at rs9399451, were generated via mutagenesis, and subsequently all plasmids were transfected. Since most lung cancer events derive from normal epithelial cells, the lung epithelial cell line Beas-2B was used in following functional genomics work.

As shown in Fig. 1, the C allele of rs9390123 showed ~24.8% lower relative luciferase activity than the T allele (P=0.0041), while the G allele of rs9399451 presented ~56.7% lower luciferase activity than the T one (P<10⁻⁵; Fig. 1), which indicated that both SNPs are functional in lung cells, and rs9399451 may have a more important effect than rs9390123 in regulating target gene expression.

Since rs9390123 and rs9399451 are not located at the promoter of any known genes, it was proposed that these two SNPs may be located within an enhancer region and could regulate the expression of target genes. This suggestion could be supported by the histone modification signal. Searching the data from the Encyclopedia of DNA Elements project (15) in the UCSC genome browser (http://genome.ucsc.edu) revealed that there are clear H3K27Ac (acetylation at the 27th lysine residue of the histone H3 protein) and H3K4me1 (mono-methylation at the 4th lysine residue of the histone H3 protein) peaks, which are two common histone modifications around an active enhancer (16), near this region in the human lung cancer cell line A549 (Fig. S2). However, the target gene remained unclear. To identify the potential target gene for this enhancer, 3C was used. The rationale of 3C is that the enhancer is far away from its target gene(s) in sequence but close in space. Thus, after 3C library construction, the target gene promoter can be ligated in higher efficiency with enhancer than other random genome regions, which can be disclosed by qPCR. To normalize the primer efficiency, the BAC RP11-1012I24 including nearby region was used as a control. There are only three protein-coding genes [PEX3, FUCA2 (α-L-fucosidase 2) and PHACTR2] and one lncRNA PHACTR2-AS1 in this BAC. Therefore, in our 3C assay, unidirectional primers were designed to anchor the promoter of these four genes, the aforementioned enhancer and several random regions. As shown in Fig. 2, no increase was observed in ligation frequency at the PHACTR2 promoter (corresponding to the 15th point in the x-axis, ~17.9 kb away from the enhancer), suggesting that the regulatory target of this enhancer is not PHACTR2. The FUCA2 promoter region (7th point in the x-axis) exhibited a similar result. By contrast, a higher ligation efficiency was obtained in the PEX3 promoter region (2nd point in the x-axis, ~173.6 kb away from the enhancer). In addition, the promoter of the lncRNA PHACTR2-AS1 (12th point in the x-axis, ~50.5 kb away from the enhancer), also presented a high ligation efficiency, indicating that PEX3 and PHACTR2-AS1 are the regulatory targets of this enhancer.

Since rs9390123 and rs9399451 are located within a non-coding region of the genome, it could be hypothesized that these two SNPs were located at a TF binding region, and could influence the interaction between TF and DNA. Since TFs usually bind to DNA through recognizing specific motif, the potential TF can be predicted by bioinformatics approach. If one TF can interact with specific DNA segment, the TF antibody can immunoprecipitate more DNA-TF complex than IgG, since IgG presents random interaction with TF. The surrounding region containing T of rs9390123 and rs9399451 are matching with the motif of POU2F1 (see https://jaspar.genereg.net/ for detail). In contrast, both C of rs9390123 and G of rs9399451 can disrupt the recognizing site for POU2F1, thus resulting much lower scores for POU2F1 interaction in prediction (data not shown). Therefore, it was hypothesized that rs9390123 and rs9399451 effect through influencing POU2F1 binding. To confirm this suggestion, ChIP with POU2F1 antibodies was utilized to identify the surrounding region of the potential binding TF. As shown in Fig. 3, the anti-POU2F1 antibody significantly pulled down more chromatin for the region surrounding rs9390123 and rs9399451 compared with that pulled down by IgG (P=0.001 and 0.004, respectively), thus supporting that POU2F1 has the ability to bind these two regions in lung cells.

To further examine the function of POU2F1 on gene transcription, a POU2F1 overexpression vector was constructed and transfected into Beas-2B cells and mRNA levels were evaluated by RT-qPCR. After transfection of the overexpression plasmid, POU2F1 expression was increased ~1.4×10⁷ fold (P<10⁻⁷; Fig. 4A), which verified that the transfection efficiency was very high. This could further induce ~95.0 and ~110% increase in mRNA expression of PEX3 (P=0.0002; Fig. 4B) and PHACTR2-AS1 (P=0.012; Fig. 4C) in lung cells, respectively, suggesting that POU2F1 can influence the expression of target genes.

Our prediction, ChIP and overexpression assay verify the effect of POU2F1 in the expression of these two genes. To further investigate the binding efficiency between two alleles at these two SNPs, EMSA was performed based on biotin-labeled probes. As shown in Fig. 5, the two alleles for rs9390123 and rs9399451 showed an apparent different affinity to nuclear proteins from Beas-2B cells. Furthermore, these patterns could be abolished by adding competitor oligonucleotides (Fig. 5), which verified the hypothesis that these two SNPs could alter TF binding affinity.

Since PHACTR2-AS1 is the regulatory target of the enhancer containing rs9390123 and rs9399451, the function and transcript location of this gene in lung cells were investigated. To determine the location of the lncRNA PHACTR2-AS1, nuclear and cytoplasmic RNA was isolated and quantified separately. As shown in Fig. S3, this lncRNA is predominantly expressed in the nucleus. To investigate the function of PHACTR2-AS1 in lung cells, a stable knockdown cell line was constructed using Beas-2B cells by transfecting a lentivirus containing shRNA, and the knockdown efficiency is displayed in Fig. 6A. In Fig. 6A, the left panel is the PCR result for GAPDH and the right panel is for PHACTR2-AS1. For each part, the four lanes, from left to right, are ladder (L), Mock (MC), negative control (NC) and knockdown (KD) cells, respectively. The arrow in Fig. 6A points out the position of the PHACTR2-AS1 PCR product. As shown in Fig. 6A, PHACTR2-AS1 expression is persistent in the NC but decreased to undetectable amount in the KD cells, which indicates that the transfection efficiency was high and the knockdown was successful. Cell proliferation was evaluated with an MTT assay. As shown in Fig. 6B, cell proliferation was significantly inhibited in the KD group compared with that of the control group (P=0.0017 and <10⁻⁵ for days 3 and 4, respectively). The number of formed clones was also significantly decreased after PHACTR2-AS1 knockdown (P=0.00085; Fig. 6C). Migration and wound healing abilities were also significantly attenuated after knockdown (P=0.00048 and 0.00024, respectively; Fig. 6D and E).

To substantiate the function of PHACTR2-AS1 in tumorigenic lung cells, we further transfected lung cancer cell line A549 with the lentivirus. As shown in Fig. 7A, PH2ACTR2-AS1 expression was similar in the mock and negative control (NC) (P>0.05). After RNA interference, PH2ACTR2-AS1 expression was decreased ~40% (Fig. 7A) compared with the mock and NC, which indicates that the transfection efficiency was high. Consequently, the cell count was reduced significantly at day 2, 3 and 4 (all P=0.0012, <10⁻⁵ and 0.0011, respectively; Fig. 7B). Meanwhile, the clone formation (P=0.00012; Fig. 7C), migration (P<10⁻⁶; Fig. 7D) and wound healing (P=0.00062; Fig. 7E) abilities were also significantly inhibited. Overall, these results indicated that PHACTR2-AS1 can promote lung cell proliferation and migration in vitro, and suggested that it could be a novel oncogene for lung cancer.