MiR‐3130‐5p is an intermediate modulator of 2q33 and influences the invasiveness of lung adenocarcinoma by targeting NDUFS1

Abstract Genome‐wide association studies (GWAS) have reported a handful of loci associated with lung cancer risk, of which the pathogenic pathways are largely unknown. We performed cis‐expression quantitative trait loci (eQTL) mapping for 376 lung cancer related GWAS loci in 227 TCGA lung adenocarcinoma (LUAD) and reported two risk loci as eQTL of miRNA. Among the miRNAs in association with lung cancer risk, we further predicted and validated miR‐3130‐5p as an intermediate modulator of risk loci 2q33 and the tumor suppressor NDUFS1. We assessed the phenotypic impacts of the interaction between miR‐3130‐5p and NDUFS1 in both lung cancer cell lines and mice xenograft models. As a result, miR‐3130‐5p directly regulates the expression of NDUFS1 and the corresponding tumor invasiveness, migration and epithelial‐mesenchymal transition (EMT). Our findings provide important clues for the pathogenic mechanism of 2q33 in lung carcinogenesis which informs clinical diagnosis and prognosis of LUAD. We performed a cis‐eQTL analysis for 376 lung cancer risk loci based on the expression profiles of 251 miRNAs in a cohort of 227 TCGA lung adenocarcinoma. We report a novel pathogenic pathway of 2q33 via miR‐3130‐5p and NDUFS1.


| INTRODUCTION
Genetic variations cause multiple human complex diseases, including cancers. GWAS have identified a vast amount of SNPs associated with diverse human traits by comparing large cohort of cases and healthy controls. 1,2 Previous GWAS reported approximate 680 risk loci for lung cancer and lung-related diseases. However, most of the trait-associated loci (TAL) are located in the noncoding genomic regions 3 and do not change the amino acid sequences of the protein, suggesting that these GWAS risk SNPs act as regulatory elements via influencing gene expression. Therefore, understanding the expression regulatory activities of TAL in non-coding regions becomes a demanding question of medical genetics and a milestone of the clinical implications. eQTL serves as a bridge between non-coding germline variants and variations of phenotypes, and thereby reveals the molecular mechanisms of TAL. 4 The rationale underlies eQTL is that germline variations alter the transcription abundances of mRNA or other transcripts to cause phenotypic changes in the cell. Many studies demonstrated that eQTL can further affect a wider spectrum of transcripts by targeting at genes with trans-acting regulatory activities. 5,6 Several large data consortia have been established by eQTL mapping based on large parallel sequencing projects. For example, the "Genotype-Tissue Expression" (GTEx) Project identified eQTLs in normal human tissues; and Gong et al. established "PancanQTL," "ncRNA-eQTL," and "GPSno" databases for cis-and trans-eQTLs of non-coding RNAs, protein coding RNAs and small nucleolar RNAs respectively in 33 cancer types from TCGA. [7][8][9][10] In particular, many recent studies successfully elucidated the susceptibility pathways of cancer risk loci by eQTL analyses in cancer cohorts, such as 15q25.1 and 8q24. [11][12][13] The landscape of eQTLs is highly specific to cell lineage, differentiation statuses as well as tissue types. 14 Therefore, eQTL mapping for cancer expression is highly informative of cancer biology. Several prior studies identified the regulatory axes for eQTLs in cancers and thereby revealed the biological processes that drive carcinogenesis. The discovery of eQTLs in cancers and the corresponding regulatory axes can inform cancer studies in two ways: first, eQTLs can be used to elucidate the functions of the non-coding cancer risk loci and thereby fine-map the causal pathogenic alleles; second, the effectors of the cancer eQTL are highly potential oncogenes or tumor-suppressors and serve as biomarkers in cancer diagnosis and therapy.
MiRNA is a class of endogenous small non-coding RNAs which functions as transcription regulator to inhibit gene expression by directly binding to the miRNA response element (MRE) in the 3'UTR of target genes. Accurate computational identification of miRNA targets is particularly important to investigate biological mechanisms. Many software and web servers have been developed to predict miRNA targets mainly based on evolutionary conservation of binding sites. [15][16][17][18] Dysregulation of miRNAs has been associated with diverse diseases, including lung cancer. 19 Several miR-NAs have been used as novel biomarkers and therapeutic targets of cancers. 20,21 Here in our study, we performed a cis-eQTL mapping for a set of lung-related TAL from GWAS and reported an association between 2q33 and a cancer-associated miRNA, miR-3130-5p. We further investigated the function impacts of the association in lung adenocarcinoma (LUAD). As a result, we reported NDUFS1 as a target gene of 2q33-miR-3130-5p association in LUAD and validated the corresponding phenotypic impacts on invasion, migration and EMT based on in vitro and in vivo analyses.

| Data and materials
We obtained 376 risk loci associated with traits of lung cancer from the GWAS catalog on 10 October 2017 (Table S1). 22 Then we retrieved the genotype data from 490 TCGA lung adenocarcinoma. To control the population structure, we only included Caucasian descendants. We inferred the ancestry of the cohort using 2504 individuals from 1000 Genomes project with known ancestry as a reference and followed the method described in previous study. 5 For each locus, we imputed the genotypes for variants within 1 Mb range of the tag SNP and with linkage disequilibrium (R 2 > 0.5) to the 1 KG variants with 1000 Genomes Phase 3 23 using IMPUTE2 software. 24,25 In order to control for statistical power, we filtered the imputed variants according to the following criteria: (a) info score >0.7; (b) minor allele frequency (MAF) >5%; (c) posterior probability >0.7; (d) Hardy-Weinberg Equilibrium test p value >1 × 10 −6 .
We obtained the miRNA expression profiles of 227 TCGA-LUAD, with matched somatic copy number alteration (SCNA) data and CpG methylation data. 26 The miRNA expression levels were transformed based on log 2 (RPM+1). To control the effect of outliers, we excluded individuals of which the miRNA expression levels were greater than three times of the standard deviation. Then we removed miRNA which were absent in over 20 individuals. For the simplicity of the computation, the copy number variation of a given miRNA the SCNA was categorized into gain (+1), neutral (0), loss (−1). For the CpG methylation measure, we used categorized levels of 0-3 based on the beta values of probes located within 200 and 1500 bp of the transcription start sites (TSS200 and TSS1500).

| Cis-eQTL analysis for TAL in relation with LUAD
We performed cis-eQTL mapping for 3291 SNPs in LUAD using R package, "MatrixEQTL". 27 For each SNP, we assessed the associations between the genotypes and the expression levels of miRNA located within 1 Mb either side of the SNP. To control confounding effects, we adjusted the expression levels of miRNA using somatic copy number variation, CpG methylation, sex and age of patients via the following multivariate regression model: Here, i ∼ N 0, 2 is a Gaussian error term; G i is the genotype of the ith sample; Age i is the age of the ith sample; Sex i is the sex of the ithsample; CNV i is the copy number variation of the ith sample; CpG i is the methylation status of the ith sample; 0 is the intercept, other are the regression coefficients. We called significant cis-associations between eSNP and the target miRNA (eMiR) based on a false-discovery rate (FDR) of 0.1.
To identify differently expressed miRNA in LUAD, the expression levels of miRNAs which were absent in more than 10 samples in normal or tumor samples were excluded. The rest were compared between tumors and normal lungs using t test. MiRNAs with fold change >2 and FDR <0.01 were defined as differently expressed miRNAs.

| Target gene prediction
In order to predict the downstream target mRNAs of eMiR, we used miRNA Data Integration Portal (mirDIP 4.1) 18 to predict target mRNAs with corresponding miRNA binding motifs in the 3'UTR regions. Then we confirmed the endogenous effect of the eSNP on the mRNA using Mendelian randomization (MR) with the eSNP as instrumental variable (Equation 2). We selected significant candidate mRNAs based on the following criteria: (a) Target genes ("Mirtarget gene") Integrated Score ≥0.1; (b) filtering these genes ("Present calls") with no less than 100 present calls in the cohort; (c) genes to which the TAL is significant genetic instrument based on Mendelian Randomization ("MR gene", weak instrument test FDR <0.05); (d) Target genes are required to be significantly differently expressed ("DE gene") between tumor tissue and the matched normal (FDR<0.01 and log 2 fold change>0.5); (e) significantly predictive ("Predictive gene") of the overall survival rate using KM-Plotter software. 28 Here, i ∼ N 0, 2 is a Gaussian error term; G i is the genotype of the ith sample; 0 is the intercept,miRNA i is miRNA expression level of the ith sample, 1 is the regression coefficients.

| Dual-luciferase reporter assays
The 3′-UTR fragments of NDUFS1 and RGS11 predicted to interact with miR-3130-5p were amplified and cloned into pmirGLO Dual-Luciferase miRNA Target Expression Vector. A549 cells transduced by those plasmids were seeded in a 24-well plate and transfected with miR-3130-5p agomir, miR-3130-5p antagomir and controls for 48 h. The cells were lysed and the luciferase activities of firefly and renilla were examined with Dual-Luciferase Reporter assay kit following the manufacturer's instructions. All experiments were performed three times. (1)

| Patients and tissue samples
Primary LUAD tissues and paracancerous tissues were obtained from 43 patients who received lung cancer resections at the Department of Thoracic Surgery of Zhongshan Hospital affiliated to Xiamen University. Written informed consents were acquired by all patients participated in the study. Our protocols were approved by Medical Ethic Council for Researchers of Zhongshan hospital and executed according to the Declaration of Helsinki. The diagnoses of LUAD were made by two independent pathologists. The resected tissue samples were quickly frozen in liquid nitrogen and stored at −80°C.

| RNA extraction and qRT-PCR
Total RNAs were extracted from cells and human tissue samples using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and purified using the PureLink RNA Mini kit (Invitrogen) according to the manufacturer's protocols. Approximately 3.5 µg of RNA was reverse-transcribed to DNA templates for qRT-PCR using the Reverse-Aid™ First Strand cDNA Synthesis kit (Fermentas, Waltham, MA, USA).

| Cell viability assay
Cells (4 × 10 3 ) were placed in 96-well plates and assessed every 24 h. CCK-8 reagent (Sangon biotech, China) was added and incubated with cells for 1 h at 37°C before the absorbance at 450 nm was measured by a microplate reader.

| Transwell assay
Cells (5 × 10 4 ) suspended in 200-μl medium without serum were seeded into the upper chamber, and lower chamber was full of 20% FBS to induce LUAD cells to migrate through the membrane. Matrigel (1:6 dilution, BD Biosciences, USA) was added on the upper chamber for the invasion assay. Twenty-four hours later, cells were fixed in 4% paraformaldehyde, stained with 0.1% crystal violet solution and counted under optical microscope.

| Wound-healing assay
Cells (2 × 10 6 ) were seeded and incubated in fresh sixwell plates. For the wound-healing assay, a confluent monolayer was scratched with a pipette tip, washed with PBS and incubated in culture medium supplemented with 1% FBS. The cultures were imaged by microscope after 0 and 48 h.

| Colony formation assay
Cells (200) were seeded and incubated in fresh 6-well plates for 14 days to allow colony formation. Colonies were fixed in 4% paraformaldehyde, stained with 0.1% crystal violet solution and counted.

| Protein extraction and western blotting
Cells or tissue samples were lysed in ice-cold RIPA lysis buffer containing protease inhibitor. An equal amount of total protein from each sample (20 µg) was loaded on 10% SDS-PAGE gels and then transferred to PVDF membranes. The membranes were blocked with 5% BSA for 2 h at room temperature and then incubated with commercially available primary antibodies for detection of NDUFS1 (Cell Signaling Technology, USA), followed by the appropriate secondary antibodies (Cell Signaling Technology, USA). GAPDH was taken as internal control and protein bands were analyzed by Image J software (NIH, USA).

| Mice xenograft experiment
A549 cells (2.5 × 10 6 ) were injected subcutaneously into the flanks of Balb/c nude mice aged 5 weeks. After formation of tumor (4 ×4 mm), mice were treated with either control or miR-3130-5p agomir (2 nmol twice a week) via tail vein for totally 4 weeks. The mice's condition, body weights and sizes of tumors were observed and recorded every three days. The tumor volumes were recorded in (width 2 × length)/2. After 2 weeks from last treatment, the animals were euthanized and tumors were removed, weighed, and frozen or fixed for biochemical or histological analyses.

| Immunohistochemistry
The paraffin-embedded tissue slides (4 μM) were deparaffinized and hydrated. Tumor sections were immunostained subsequently for NDUFS1, CD34, Collagen-IV. MaxVisionTM kit was applied to each slide (MaixinBiol, China) according to the manufacturer's instructions. Color was developed by DAB kit (MaixinBiol, China) and slides were counterstained with hematoxylin. Nonimmune mouse IgG was used to replace primary antibodies for negative controls. The immunostained slides were observed under microscope and scored by two independent pathologists using Image J.

| Clinical data collection and statistical analysis
The medical records of patients enrolled in the study were reviewed and approved by Medical Ethic Council for Researchers of Zhongshan Hospital Affiliated to Xiamen University with patients' written consents. The information of demographic characteristics, pathological characteristics of primary tumors, diagnoses and therapies were collected. These patients were followed up every 3 months until deaths or end of the investigation.
All the experiments were conducted in triplicates with consistent results. All data were listed as mean with SD and analyzed by SPSS 20.0 (IBM, USA) and Graphpad Prism 7.0 (GraphPad Software Inc.). Student's t test was adopted to confirm statistical significance. p value less than 0.05 was regarded as statistically significant. We used logistic regression to evaluate the association between NDUFS1 or miR-3130-5p and corresponding clinical features and performed survival analysis based on cox regression using "survival" and "survminer" R packages (version: R-4.0.2).

| Cis-eQTL analysis for GWAS traitassociated loci in TCGA LUAD Cohort
We retrieved 376 loci which were annotated to associate with disease traits in human lung, each of which is represented by a tag SNP (Table S1). Based on each tag SNP, we obtained 3291 SNPs in linkage disequilibrium (R 2 > 0.5). The missing genotypes of the SNPs were imputed from the genotyping data of TCGA cohort.
To define the study population, we stratified 490 TCGA LUAD samples based on the top two principal components of common variants (see Section 2). After filtering for samples with missing or extreme values, we obtained a cohort of 227 Caucasian descendants with LUAD.
We then performed a cis-eQTL analysis for 3291 SNPs and 251 miRNAs located within 1 Mb region on either side of the SNPs. Based on an FDR of 0.1, we identified 100 significant cis-associations which correspond to 93 eSNPs and 3 target miRNAs (eMiR) ( Figure 1A-C and Table S2). The eSNPs belong to two risk loci (2q33 and 8p22) and are both associated with lung cancer. 29,30 We evaluated the expression levels of three eMiRs in LUAD and matched normal samples. As a result, miR-3130-5p (log 2 fold change = 1.04, FDR =8.12 × 10 −7 ; Figure 1D) and miR-548v (log 2 fold change = 0.26, FDR = 0.01; Figure S1) were significantly differently expressed. Moreover, miR-3130-5p expression level is significantly predictive of the overall survival rate in a cohort of 323 TCGA LUAD patients (p = 0.024; Figure S2A). Therefore, we focused on miR-3130-5p for further analysis.
To reveal the downstream impacts of the TAL which act as cis-eQTL of miR-3130-5p, we retrieved 1732 mRNAs (Integrated Score ≥ 0.1; Table S3) predicted by mirDIP 4.1 database ( Figure 1E). These genes carry the binding motif of miR-3130-5p in the 3′UTR. After filtering these genes for present calls in the cohort, we yielded 1330 target genes for further analysis ( Figure 1E). We performed instrumental variable (IV) regression taking the eSNPs as genetic instruments, miRNA expression as independent variable and mRNA expression as dependent variable. As a result, we yielded 11,756 regulatory axes, including 23 eSNPs of miR-3130-5p and 596 mRNAs in trans (Table S4 and Figure 1E). To identify the mRNAs regulated by miR-3130-5p that are associated with LUAD patients' survivals, 256 significantly differently expressed genes in lung cancer were chosen for further survival analyses ( Figure 1E). As results, 189 genes significantly predictive of overall survival (FDR <0.05; Figure 1E) were selected, among which NDUFS1 is the most significantly predictive (HR: 0.46, CI 95%: 0.38-0.55, FDR = 4.25 × 10 −5 ; Figure S2B and Table S5). Other candidate genes with known functional impacts in cancer include RGS11, 31 RPS6KA2, 32 and CYB5A. 33 KEGG enrichment analysis suggested the 189 candidate genes overrepresent in two cancer-related pathways (Wnt signaling pathway and Cell cycle) at a 5% FDR level (Table S6). Based on literature research, 31,34 we chose two genes, NDUFS1 and RGS11, for further analyses.

| NDUSF1 is a target for miR-3130-5p
To verify the target gene of miR-3130-5p, we performed luciferase reporter assays for two predicted targets, NDUFS1 and RGS11 (Figure 2A and Figure S3B). In A549 cells, the reporter with wild-type NDUFS1 3′UTR exhibited significantly attenuated luciferase activities in response to the overexpression of miR-3130-5p (agomir transfection; | 3705 ZHAN et Al. p < 0.05; Figure 2B) but showed no response to the knockdown (antagomir transfection; p > 0.05; Figure S3A). In addition, the changes in the luciferase activities caused by the overexpression and knock-down of miR-3130-5p were further validated by the changes of NDUFS1 mRNA levels ( Figure 2C) and protein levels ( Figure 2D). For RGS11, the luciferase activities were not significantly affected by miR-3130-5p ( Figure S3C). In summary, our data suggested that NDUFS1 is the direct target mRNA of miR-3130-5p in LUAD.

| miR-3130-5p promotes invasion, migration, and EMT in vitro
To further verify that miR-3130-5p functions as an on-comiRNA in LUAD, we first assessed the expression levels of miR-3130-5p in a panel of cell lines derived from LUAD. Our data showed miR-3130-5p is stably expressed in A549 and H1650, which we chose for further analysis. We then conducted a series of in vitro experiments to assess the effects of miR-3130-5p on the malignant phenotypes in the two cell lines. As a result, miR-3130-5p overexpression (agomir) significantly promoted the invasion and migration of A549 (p < 0.01) and H1650 (p < 0.001 and p < 0.01) cells after transfection for 24 h (Figure 3A and B). Consistently, the knock-down of miR-3130-5p by antagomir substantially suppressed the two phenotypes. Additionally, the wound healing assays also ascertained that miR-3130-5p agomir greatly intensified the migration abilities of A549 (p < 0.01) and H1650 (p < 0.0001) cells, which were significantly reduced by miR-3130-5p antagomir (A549: p < 0.0001 and H1650: p < 0.01) ( Figure 3C). Furthermore, the western blotting assays showed suppressed epithelial markers (E-Cadherin) and activated mesenchymal markers (N-Cadherin, Vimentin, Snail) in A549 and H1650 cells driven by increased miR-3130-5p expression (p < 0.01; Figure 3D and S5), indicated miR-3130-5p may promote EMT by regulating the expression of Snail.
On the other hand, the CCK8 and colony formation assays suggested that miR-3130-5p does not affect cell proliferation and colony formation of LUAD cells ( Figure S4A and B). And miR-3130-5p had no effect on the apoptosis process of LUAD cells according to the TUNEL assay ( Figure S4C). In summary, our findings suggest that miR-3130-5p specifically promotes invasion, migration and EMT of LUAD cells.
To further verify miR-3130-5p exerts its function through modulating the activities of NDUFS1, we performed rescue experiments in A549 and H1650 cells. As a result, miR-3130-5p overexpression significantly counteracted the inhibition of invasion, migration abilities by NDUFS1 in vitro (p < 0.01; Figure 4A-C). In A549 cells, miR-3130-5p restored the repressed level of Snail and reduced the increased level of E-cadherin caused by overexpression of NDUFS1. While in H1650 cells, miR-3130-5p regained the expression of Snail which suppressed by overexpression of NDUFS1 ( Figure 4D and Figure S6). In summary, our data suggests that miR-3130-5p is a strong modulator of NDUFS1 activities in vitro and thereby influences the invasion, migration and EMT of LUAD cells.

| miR-3130-5p promotes tumor invasion in vivo
The in vitro experiments suggest that miR-3130-5p functions as a modulator of NDUFS1 with phenotypic impacts in lung cancer. Therefore, we further verified the functions of miR-3130-5p in vivo. We established xenograft tumor animal models by subcutaneously injection of A549 cells in BALb/c mice. We intravenously injected miR-3130-5p agomir and control to assess the pathological changes. As a result, there were no obvious differences in the tumor volume and body weight of mice between miR-3130-5p agomir group and control group (p > 0.05; Figure 5A and B), which is consistent with our previous observation in vitro. On the other hand, the immunohistological results showed markedly increased expression of invasion related markers, such as CD34 and collagen IV, in tumors formed in miR-3130-5p overexpressing mice compared to the control group (p < 0.01; Figure 5C), which suggests that miR-3130-5p promotes tumor invasiveness in vivo.

| DISCUSSION
Numerous studies suggest TAL act as eQTL of transcripts with regulatory activities, such as miRNA, lncRNA, or circRNA and thereby influence gene expression in cancer cells. Such regulatory axes are used to explain the biological process underlying non-coding TAL and to inform new cancer-related transcripts. In the current study, we identify miR-3130-5p as a cis-target of TAL 2q33. Further analysis suggests miR-3130-5p is a potential oncogenic micro-RNA, which intermediates the effects of 2q33 risk locus on NDUFS1 and influence the invasion phenotype of LUAD.
2q33 is reported as a pleiotropic TAL which is associated with nicotine dependence, urinal mercapturic acids concentration in smokers as well as risk of lung cancer and esophageal cancer in Northern Chinese. 30,35,36 Therefore, the eQTL activity of 2q33 helps understand the pathogenic mechanism of lung cancer. MiR-3130-5p is located 70,115 bp from the tag SNP, rs75358501. Prior studies showed that miR-3130-5p was associated with diverse diseases. [37][38][39][40][41][42] It was reported in genomic researches about melanoma and breast cancer. 39,40 Different expression of miR-3130-5p may be involved in breast cancer neoadjuvant chemotherapy resistance. 41 However, the function of miR-3130-5p is less known in other cancer types. To our knowledge, the current study is the first report of miR-3130-5p as a cis-target of GWAS loci and associated with lung carcinogenesis with a functional role.
Instrumental variable (IV) regression provides a powerful analytical tool for identifying causal relations among germline variants, intermediate transcripts and downstream effectors in trans. Many prior studies successfully used IV regression to elucidate the functional impacts of germline variants. Our study uses IV regression to suggest NDUFS1 as an effector in trans of 2q33-miR-3130-5p association, which avoids the high false discovery rates (FDR) in T A B L E 1 Analysis of the correlation between expression level of NDUFS1 in LUAD tissues and clinicopathological characteristics F I G U R E 5 MiR-3130-5p promotes tumor invasion in vivo. (A) Mice xenograft models were constructed by injecting A549 cells subcutaneously into nude mice. After the tumors developed to 4 ×4 mm, miR-3130-5p agomir and control were injected into mice regularly until mice were sacrificed 4 weeks later. Tumor tissues were resected and examined. (B) There were no significant differences in the body weight and tumor volume between miR-3130-5p agomir group and control group. (C) The protein levels of NDUFS1, CD34, and collagen IV were detected by IHC and compared between miR-3130-5p agomir group and control. *p < 0.05, **p < 0.01 direct trans-eQTL analysis and the low statistical power due to limited sample size. The method we described can be used to explain other non-coding cancer risk variants. In addition, IV regression can also provide a framework to infer causal variants from a set of variants in linkage disequilibrium.
As the downstream effector of 2q33-miR-3130-5p, NDUFS1 is the largest subunit of mitochondria complex I, which plays important role in the process of oxidative phosphorylation and production of ATP. It is well-known that dysfunction of mitochondria complex I can cause Warburg effects in cancer cells, which means accelerated aerobic glycolysis and cancer progression. 43,44 NDUSF1 is suggested to exert anti-cancer effect in ovarian cancer, 45 clear-cell renal-cell carcinoma 46 and NSCLC. 34 Consistent with previous studies, our study demonstrated NDUFS1 acts as a tumor-suppressor in the invasion, migration and EMT of LUAD cells. Moreover, low expression of NDUFS1 in cancerous tissue is significantly associated with severe lymphatic and vascular invasion pathologically. Thus, NDUFS1 may be a potential biomarker for the clinicopathological diagnosis of LUAD. However, we had not found NDUFS1 level is associated with the OS and PFS of LUAD patients. The possible reason is the heterogeneity in the cohort, that is, patients received different treatments after resection. With limited sample size it is difficult to assess the predictive power of NDUFS1. However, in a larger cohort such as TCGA, we show NDUFS1 expression is significantly associated with clinical outcome of LUAD. Therefore, NDUFS1 is a potential tumor-suppressor targeted by 2q33-miR-3130-5p axis which influences the invasiveness of LUAD. Our data suggests that NDUFS1 serves as a biomarker for diagnosis and prediction of LUAD.
The current study provides a thorough analyses to confirm the functional impacts of 2q33 and miR-3130-5p association through altering NDUFS1 transcription. However, the effect of the regulatory axis needs to be compared with known LUAD biomarkers, such as driver mutations in EGFR, ALK and ROS1. 47 Although we validated the interaction between miR-3130-5p and NDUFS1, there are multiple mechanisms through which miRNA alters the transcript abundances of mRNA. For example, miRNA affects the decay of mRNA, and induces transcriptional gene silencing in mammalian cell (RITS) by directly targeting the promoter. The current analysis does not inform the mechanism through which miR-3130-5p inhibit NDUFS1 transcription. In addition, miR-3130-5p may interact with other transcripts in LUAD and thereby influences the phenotype of cancer. For example, LncRNA and circRNA competitively bind with miRNA and regulate the expression of downstream target mRNA. 48,49 Therefore, it is possible that the pathogenicity of miR-3130-5p is through other complex pathways in LUAD.
In conclusion, we describe an eQTL analysis in LUAD to reveal the biological mechanism of TAL 2q33. We demonstrate that miR-3130-5p is an oncomiRNA which is significantly associated with 2q33 and a modulator of NDUFS1 expression in LUAD. In addition, miR-3130-5p suppresses the activity of NDUFS1 and thereby promotes tumor invasiveness and possibly migration and EMT hence a potential biomarker of LUAD. Our analyses identified and verified a novel pathogenic pathway for 2q33 in LUAD through miR-3130-5p-NDUFS1 axis. In addition, we foresee a couple of steps in the future study. First, we can further verify the function of miR-3130-5p and NDUFS1 in knockout cell lines and mice. Then, we can establish a large, wellcontrolled LUAD cohort to evaluate the clinical significance of 2q33, miR-3130-5p, and NDUFS1. We can also verify the invasiveness phenotypes using methods such as in vivo clonogenic assay. These findings will extend our knowledge of lung carcinogenesis and inform the clinical management of LUAD.