FGFR1–4 RNA-Based Gene Alteration and Expression Analysis in Squamous Non-Small Cell Lung Cancer

While fibroblast growth factor receptors (FGFRs) are involved in several biological pathways and FGFR inhibitors may be useful in the treatment of squamous non-small cell lung cancer (Sq-NSCLC), FGFR aberrations are not well characterized in Sq-NSCLC. We comprehensively evaluated FGFR expression, fusions, and variants in 40 fresh-frozen primary Sq-NSCLC (stage IA3–IV) samples and tumor-adjacent normal tissues using real-time PCR and next-generation sequencing (NGS). Protein expression of FGFR1–3 and amplification of FGFR1 were also analyzed. FGFR1 and FGFR4 median gene expression was significantly (p < 0.001) decreased in tumors compared with normal tissue. Increased FGFR3 expression enhanced the recurrence risk (hazard ratio 4.72, p = 0.029), while high FGFR4 expression was associated with lymph node metastasis (p = 0.036). Enhanced FGFR1 gene expression was correlated with FGFR1 protein overexpression (r = 0.75, p = 0.0003), but not with FGFR1 amplification. NGS revealed known pathogenic FGFR2,3 variants, an FGFR3::TACC3 fusion, and a novel TACC1::FGFR1 fusion together with FGFR1,2 variants of uncertain significance not previously reported in Sq-NSCLC. These findings expand our knowledge of the Sq-NSCLC molecular background and show that combining different methods increases the rate of FGFR aberrations detection, which may improve patient selection for FGFRi treatment.


Introduction
Lung cancer is the most common cause of cancer-related death worldwide [1]. The squamous non-small cell lung cancer histotype (Sq-NSCLC) accounts for 20-30% of nonsmall cell lung cancer. The only approved novel first line systemic therapy for Sq-NSCLC is immune checkpoint inhibitors [2,3]. Therefore, it is important to identify effective targeted therapies and reliable predictive molecular biomarkers for Sq-NSCLC patients.
Recent studies have indicated the potential of fibroblast growth factor receptor inhibitors (FGFRis) as treatment due to the high rate of fibroblast growth factor receptor (FGFR) aberrations found in targetable oncogenic pathways [4] (reviewed in [5][6][7]). However, most early-phase clinical trials with FGFRis have shown only a partial response, which may be a result of the poor predictive power of FGFR1 amplification, which was initially the only predictive biomarker of response to FGFRis in Sq-NSCLC [7][8][9][10][11]. Subsequent studies revealed that FGFR mutations, fusions, or expression might provide more precise information on potential responders. In this context, FGFR gene variants and expression have been evaluated in several preclinical and clinical studies, albeit with conflicting results. For instance, FGFR mRNA or protein levels were reported as significant predictors of sensitivity to AZD4547 [12], BGJ398 [13], ponatinib [12,14], and rogaratinib [15], while others showed inconsistent results [10,16,17].
Members of the fibroblast growth factor (FGF) family, the FGFR1-4 genes, encode four highly conserved tyrosine kinase receptors (FGFR1-4). Each FGFR consists of an extracellular region composed of three immunoglobulin-like domains (Ig I-Ig III), a single hydrophobic transmembrane domain, and a cytoplasmic tyrosine kinase domain. The extracellular domains interact with FGFs, which leads to dimerization of the FGFR followed by activation by sequential autophosphorylation of tyrosine residues [7]. FGFR signaling activates the phosphoinositide-3-kinase (PI3K)/AKT, signal transducer and activator of transcription (STAT), and mitogen activated protein kinase (MAPK) pathways [18]. Deregulated FGF/FGFR signaling through FGFR gene amplification, mRNA overexpression, mutation, or gene fusion is associated with ligand-independent dimerization of FGFRs and subsequent activation of cancer-related signaling pathways (PI3K/AKT, STAT, and MAPK), affecting cell proliferation, survival, metabolism, migration, and the cell cycle. While FGFR1-4 play roles in several biological pathways, simultaneous manifestation of their aberrations and expression in Sq-NSCLC has not been well characterized. Additionally, little is known about the clinical importance of FGFR1-4 expression. Finally, the difference between tumors and the surrounding lung tissue has not yet been comprehensively explored, especially in squamous lung cancer.
Therefore, to expand our knowledge of the role of FGFRs in Sq-NSCLC as well as their use as potential biomarkers for FGFRi treatment, a comprehensive evaluation of FGFR aberrations and their clinical importance in Sq-NSCLC was conducted using realtime PCR (RT-PCR) and next-generation sequencing (NGS) to assess gene expression, fusions, and variants, as well as FGFR1-3 protein expression and FGFR1 amplification. Moreover, the difference between FGFR1-4 expression levels in primary tumors and tumoradjacent normal tissues, as well as the predictive value of FGFR1-4 mRNA expression, was investigated.
Additionally, since the anchored multiplex PCR followed by NGS can be used in gene expression studies, the FGFR1-3 expression levels obtained from the RT-PCR and NGS were compared. The analysis revealed high correlation within FGFR1, FGFR2, and FGFR3 expression assessed with both methods (r = 0.77, p = 0.000001; r = 0.86, p = 0.000001; and r = 0.95, p = 0.000001, respectively; Figure 2 and Supplementary Figure S2

Clinical Significance of FGFR1-4 mRNA Expression
Analysis of FGFR1-4 expression and disease-free survival (DFS) revealed that the increased level of FGFR3 mRNA was correlated with an increased risk of recurrence (RT-PCR: hazard ratio (HR) 4.72, p = 0.029; NGS: HR 7.9, p = 0.0049). Kaplan-Meier survival curves also showed a trend toward poorer prognosis for patients with high FGFR3 expression compared to those with low expression (Figure 3). The mean DFS time of patients with high and low FGFR3 expression was 620 days (<2 years) and 1037 days (~3 years), respectively. FGFR1 (RT-PCR: p = 0.44, NGS: p = 0.56), FGFR2 (RT-PCR: p = 0.22, NGS: p = 0.19), and FGFR4 (p = 0.97) expression levels did not affect the risk of recurrence.

Association between the FGFR1-4 Expression and Clinicopathological Features
Interestingly, analysis of FGFR1-4 expression association with clinicopathological characteristics (Table 1) revealed that the FGFR4 mRNA expression was significantly higher in patients with lymph node metastasis (median: 0.02 vs. 0.01, p = 0.036; Figure 4), while FGFR1-3 expression levels were not significantly related to any of the clinicopathological features examined.

Squamous Non-Small Lung Tumors
RNA analysis was performed for all 40 tumor and five tumor-adjacent normal samples. In the tumor group, two (5%) different ("strong" and "low confidence") FGFR gene fusions were detected ( Table 2). The first identified fusion was previously known as the FGFR3::TACC3 fusion (ex17::ex11) (M-20; Figure S3). Its exact breakpoint is indicated in the Archer Quiver database, and despite the low percent structural variation reads for GSP2 (5.49%) was categorized as "strong". The second discovered fusion event was the TACC1::FGFR1 fusion (ex2::ex2) (M-61; Figure S4), which was categorized as "low confidence" because of the low percent of structural variation reads for GSP2 (4.   FGFR variant analysis (based on the ClinVar database) revealed two pathogenic variants (5%), three variants of uncertain significance (7.5%), and one variant with no clinically significant information (2.5%) ( Table 2) Correlations of described FGFRs genetic changes with FGFRs expression level in Sq-NSCLC revealed statistically significant but weak association of FGFR2 c.2398dup (p.(Ser800PhefsTer22)) variant with the increased FGFR2 mRNA expression level (p = 0.02, r = 0.36).

Discussion
Our findings revealed significantly decreased FGFR1 and FGFR4 mRNA expression in tumor tissue compared with tumor-adjacent normal tissue. Meanwhile, expressions of both genes were enhanced in selected tumor samples, reaching values similar to surrounding normal tissue. FGFR2 and FGFR3 mRNA expression varied between tumor and tumor-adjacent normal tissues and were enhanced 2-10 times in several tumor samples. Accordingly, FGFR1 mRNA expression was lower in Sq-NCLC compared with normal tissues as per The Cancer Genome Atlas and Genotype-Tissue Expression data [19]. FGFR1 expression has recently been linked to the highly negative correlation between FGFR1 mRNA and methylation levels (average of three CpG sites: cg10823844, cg15791248, and cg27646230) in tumor samples [19]. In contrast to our results, Ren et al. [13] demonstrated that nearly 50% of tumors had an increase in FGFR1 mRNA expression compared with tumor-adjacent normal tissue; however, in their study, FGFR1 expression was normalized to GAPDH. Research has indicated that the choice of stable reference gene is crucial for accurate results in gene expression studies [20][21][22]. Our previous results assessing reference genes in Sq-NSCLC showed that GAPDH gene expression was diverse and unstable in tumor and tumor-adjacent normal tissue samples, whereas POLR2A and ACTB expression levels were the most stable among analyzed samples ( Figure S5) [22].
To our knowledge, this study is the first to show significantly lower FGFR4 mRNA expression in Sq-NSCLC compared with normal lung tissue. Huang et al. [23] as the only one compared FGFR4 mRNA expression levels in lung tumors and adjacent normal tissues with results inconsistent with our results, i.e., the FGFR4 mRNA level was significantly higher in tumor tissues. Importantly, in that study [23], the detailed tumor histotype was not mentioned and GAPDH was used for normalization. Additionally, lower FGFR4 expression at the mRNA and protein levels was reported in lung tissues obtained from patients with idiopathic pulmonary fibrosis (IPF) compared with control patients. The authors indicated that the lower FGFR expression level was related to FGFR4 downregulation by pro-fibrotic factors (TGFβ, CTGF, and ET-1) [24]. Simultaneously, our work revealed enhanced FGFR4 mRNA levels in several samples. In vitro analysis of squamous lung cancer cell lines showed that FGFR4 overexpression leads to FGFR4 auto-activation and increased cell growth, clonogenicity, soft agar colony formation [25], and enhanced cell proliferation, whereas knockdown of FGFR4 can reduce proliferation [23]. Recent studies have indicated that FGFR4 deficiency might regulate the tumor immune microenvironment by activating the antigen presentation process and cellular immunity to the change in sensitivity to immune checkpoint inhibitor treatment in NSCLC [26].
In the present study, we revealed a significant association between increased FGFR4 expression and lymph node metastasis in Sq-NSCLC, but there was no significant association between FGFR1-4 mRNA expression and other clinicopathological features. A relationship between FGFR4 and lymph node metastasis has been implied. High FGFR4 protein overexpression has been correlated with lymph node metastasis in triple-negative breast cancers [27] and gastric cancer [28,29]. Earlier studies demonstrated that FGFR4 overexpression may be a result of gene amplification, especially in breast cancer tumors with high lymph node metastases, as well as in estrogen receptor-and progesterone receptor-positive tumors [30]. Additionally, the FGFR4 p.(Gly388Arg) variant, located in the transmembrane domain, has been correlated with poorer overall and progression-free survival in Sq-NSCLC patients with lymph node involvement [31]; this can be linked to increased FGFR4 stability and sustained activation, as has been shown in prostate cancer [32].
Herein, we present the novel finding that increased FGFR3 mRNA expression might be a negative prognostic marker in terms of the risk of recurrence of squamous cell lung cancer. Due to the relatively small number of patients analyzed, our results should be interpreted with caution. Nevertheless, mRNA levels of FGFR3 mRNA have been associated with worse DFS in oropharyngeal squamous cell carcinoma (p = 0.005) [33] and in non-muscleinvasive bladder cancer (HR 3.78, p < 0.001) [34]. The FGFR3 mRNA level is also a negative prognostic factor for lung adenocarcinoma [35] and squamous cell laryngeal cancer [36], where high FGFR3 expression was significantly correlated with shorter overall survival (OS). By contrast, high FGFR3 mRNA expression was associated with better progressionfree survival in patients with primary pT1 bladder cancer (log-rank, p < 0.001) [37]. Our study did not reveal any clinical significance of FGFR1, 2, and 4 mRNA levels in terms of tumor relapse; however, the impact of FGFR1 mRNA level on squamous cell lung cancer patient survival remains controversial. For instance, elevated FGFR1 expression was reported as a negative factor that negatively impacted patient survival [38], while Wynes et al. [14] revealed no prognostic association with OS or DFS. OS is a widely used endpoint for assessing the prognostic value of studied variables. However, censoring patients at death or the date of last follow-up (for OS), especially in groups with a small number of observations, can lead to an overestimation of OS [39]. In the present study, only DFS (i.e., the time to cancer recurrence or death from any cause) was analyzed.
Additionally, the presented results show the usefulness of anchored multiplex PCR followed by NGS for FGFR1-3 gene expression analysis. The high correlation (77-95%) of FGFR expression obtained with the use of RT-PCR and an NGS panel was observed, despite the different sets of reference genes. In RT-PCR, POLR2A and ACTB were used to normalize expression, while the NGS panel contained CHMP2A, GPI, RAB7A, and VCP. Thus far, use of the relative gene expression data obtained with the applied NGS panel for genetic profiling has been described for acute myeloid leukemia [40], while it was previously limited to determining gene fusion occurrence [41,42].
Considering that FGFR mutations and fusions are detected in Sq-NSCLC (reviewed in [7]), we performed targeted FGFR sequencing with the use of anchored multiplex PCR technology followed by NGS, which was previously shown to efficiently discover genomic aberrations, including novel fusions. We found one FGFR3::TACC3 fusion (2.5%), which is the most frequent somatic translocation in Sq-NSCLC (range: 0.6-5.3%) [44,[46][47][48][49][50][51][52][53]. Interestingly, mRNA expression of FGFR3 was slightly increased (fold change = 1.5) compared with tumor-adjacent normal tissue in this sample. Nonetheless, previous reports have shown that this fusion event is not correlated with FGFR3 mRNA or protein overexpression in Sq-NSCLC [44], in contrast with glioblastomas [54]. Parker et al. [55] suggested that loss of the 3 region of FGFR3 might abolish downregulation by miR-99a and lead to overexpression of the fusion gene. We also detected a new somatic fusion event, a combination of TACC1 and FGFR1. Truncated TACC1 was fused before the FGFR1 extracellular immunoglobulin-like domain (Ig). Because this fusion did not cause truncation of FGFR1, its effect remains unknown. Additionally, FGFR1 amplification together with mRNA and protein overexpression were observed in this tumor. Two pathogenic, missense variants of FGFR2 and FGFR3, c.870G>T (p.(Trp290Cys)) and c.746C>G (p.(Ser249Cys)), were also detected at a frequency of 5%, consistent with previously published results [8,10,[56][57][58][59][60]. Both variants are located in the FGFR protein extracellular domain (Ig III) and induce constitutive dimerization and receptor activation via modest dimer stabilization in the absence of ligand [61]. Among variants of uncertain significance (based on the ClinVar database), we identified the FGFR1 (c.899T>C, p.(Ile300Thr)) and FGFR2 (c.2419G>A, p.(Glu807Lys)) missense variants and one FGFR2 frameshift (1 bp duplication) variant (c.2398dup, p.(Ser800PhefsTer22)). To our knowledge, this is the first detection of these variants in Sq-NSCLC. The FGFR1 (c.899T>C) missense variant was previously reported in craniosynostosis by Wilkie et al. [62] and nonsyndromic trigonocephaly by Kress et al. [63]. This variant results in amino acid substitution p.(Ile300Thr) located in the extracellular domain (Ig III). However, based on in silico analysis (UniProt) available via the VarSome database, it is unclear whether this variant causes a pathogenic or benign effect on the protein. Clinical significance scoring (based on the American College of Medical Genetics (ACMG) guidelines for the interpretation of sequence variants [64]) showed that it is likely a benign variant. The FGFR2 (c.2419G>A) variant results in an amino acid substitution at the end of the cytoplasmic domain of FGFR2 (p.(Glu807Lys)), aside from the kinase domain. Localization of this variant indicates a weak clinical relevance, confirmed by VarSome database scoring, suggesting a benign clinical significance. Until now, occurrence of this variant has been associated with Apert syndrome (acrocephalo-syndactyly type 1) in the ClinVar database. No further information or publications are available. The last detected variant of uncertain significance, FGFR2 (c.2398dup), results in a premature translational stop codon in the FGFR2 gene (p. (Ser800Phefs*22)). CMG classification in the VarSome database indicates a pathogenic clinical significance. Our results showed slightly association with increased FGFR2 mRNA expression level, but without protein overexpression what may arise from the issue that the frameshift is located at the end of the cytoplasmic domain of the FGFR2 protein, next to the kinase domain, and results in disruption of the last 23 amino acids of the FGFR2 protein. Therefore, its clinical relevance may be weak. This variant was also detected in tumor-adjacent normal tissue, indicating germline origin. It has not been reported previously in the literature (with the exception of the ClinVar database where it was linked with craniosynostosis), which might be a result of its low allelic frequency (~2.5%). We also report a new, single-nucleotide variant of the FGFR2 gene (c.2211G>T), which causes a missense amino acid change p.(Met737Ile) at the end of the kinase domain. ACMG classification in the VarSome database indicates a likely pathogenic effect; there is no available information on its functional and clinical significance.

Patient and Tumor Selection
This study included patients with a lung cancer diagnosis (n = 63) (Sq-NSCLC (n = 56) and adenosquamous carcinoma (n = 7)) who underwent a surgical procedure at the National Institute of Tuberculosis and Lung Diseases and for whom cancer-tissue samples (formalinfixed, paraffin-embedded (FFPE)) were obtained during the intraoperative procedure performed for diagnostic purposes. Tumor tissue samples, together with corresponding tumor-adjacent normal tissue (taken from the surgical resection margin at least 5 cm from the tumor), were snap-frozen in liquid nitrogen followed by storage at −80 • C until further analysis ("fresh-frozen" tissue samples). Cryostat sections were stained with hematoxylin and eosin and evaluated by a pathologist for cancer cell content, stromal cell contamination, and necrosis. Only tumor samples from 46 patients containing greater than 50% cancer cells and tumor-adjacent normal tissue samples without any cancer cells were selected.
Finally, the study group consisted of 40 cancer patients for whom a sufficient amount of RNA was available. All tumors were uniformly reviewed and classified histologically according to the World Health Organization guidelines [68]. Clinical stage was determined with the use of the TNM Classification of Malignant Tumors (8th edition) [69]. Clinicopathological characteristics are presented in Table 1.

RNA Extraction and cDNA Synthesis
Total RNA was extracted from fresh-frozen tissues with the use of RNeasy Plus Mini Kit (Qiagen, Hilden, Germany) according to the manufacturer's instructions. gDNA Eliminator columns allowed for the elimination of genomic DNA while avoiding RNA damage during DNase digestion. RNA quantity (260/280 ratio) was measured with the use of a NanoDrop UV spectrophotometer (ThermoFisher, Waltham, MA, USA) and Quantus fluorometer (Promega, Madison, WI, USA). Additionally, RNA quality was assessed by RNA electrophoresis with the 2100 Bioanalyzer System (Agilent, Santa Clara, CA, USA). The RNA integrity number ranged from 4.4 to 9.8 with a median value of 8.8. One microgram of total RNA was transcribed to cDNA using the High-Capacity cDNA Reverse Transcription Kit (ThermoFisher) with random primers according to the manufacturer's instructions.

Real-Time Polymerase Chain Reactions (RT-PCR)
RT-PCR was run in triplicate using the TaqMan Universal PCR Master Mix with Am-pErase™ Uracil N-Glycosylase (ThermoFisher) and approximately 10 ng of total RNA from 40 tumors and 20 corresponding tumor-adjacent normal tissue samples reverse transcribed to cDNA. Quantitative RT-PCR was run on the 7500 Fast Real-Time PCR System (Applied Biosystems, Waltham, MA, USA) with the use of the FAM-and VIC-labeled TaqMan Gene Expression Assays (ThermoFisher) for FGFR1 (Hs00241111_m1), FGFR2 (Hs01552918_m1), FGFR3 (Hs00179829_m1), FGFR4 (Hs00242558_m1), POLR2A (Hs00172187_m1), and ACTB (Hs99999903_m1). The RT-PCR results were averaged and FGFR1-4 gene expression levels were normalized to the reference genes ACTB and POLR2A ( Figure S5) [22]. Gene expression was analyzed with the use of a relative quantification method. For the 20 patients with both tumor and corresponding tumor-adjacent normal tissue samples available, data were expressed as a fold-change in expression between tumor and normal samples (−2 −∆∆Ct method).

Next-Generation Sequencing (NGS)
RNA from tumor tissue and tumor-adjacent normal tissue was sequenced and screened for gene fusions, variants, and expression of 14 genes of interest (ALK, BRAF, EGFR, FGFR1-3, KRAS, MET, NRG1, NTRK1-3, RET, and ROS1) using the FusionPlex Lung kit (Archer Dx, Boulder, CO, USA). Briefly, RNA (68-250 ng) was transcribed to cDNA using random priming. Next, cDNA quality was checked with the PreSeq RNA QC Assay (Archer Dx). Only cDNA with PreSeq result CP < 28 was used for DNA library construction according to the manufacturer's instructions. Subsequently, concentration and quality of obtained libraries were determined using the KAPA Universal Library Quantification Kit (Roche Diagnostics, Basel, Switzerland). Next, libraries were normalized, multiplexed, and sequenced using the MiSeq Reagent Kit, v3 (600 cycles) (Illumina, San Diego, CA, USA) on the MiSeq platform. Samples were sequenced in four runs with average quality parameters: QC30 of 85.4% and cluster density of 1324 k/mm 2 . The Illumina MiSeq sequencer generated paired-end sequence reads with an average of 1.15 million (range: 0.38-3.5) reads per sample (detailed sample statistics are shown in Table S2). NGS results were analyzed in Archer Analysis software v6.2 (Archer Dx), which aligns sequencing data against the human genome (version hg19) with the use of BWA (Burrows-Wheeler Alignment Tool) and Bowtie 2 (an ultrafast, memory-efficient short read aligner) for mapping. For the quality check (QC metric), four control genes (CHMP2A (charged multivesicular body protein 2A), GPI (glucose-6-phosphate isomerase), RAB7A (RAB7A, member RAS oncogene family), and VCP (valosin containing protein)) served as reliable indicators of overall RNA quality and content in the sample. A QC metric of at least 10 was required to support the targets of the assay. The median QC metric of analyzed samples was 503.88 (range: 242.5-607.4).
Gene fusions were called with the following detection limits: total number of supportive reads spanning the fusion junction ≥ 5; number of unique start sites for the fusion sequence specific primer ≥ 3; and percent of supporting reads at breakpoint supporting fusion ≥ 10%. Analysis of a control sample (described below) revealed that the percent of reads at breakpoint supporting fusion could be lowered to 3.6%. Results of the fusion annotation were split into two categories: "strong confidence" fusion and oncogenic isoform candidates and "low confidence" fusion candidates.
Gene variants were called and listed when the altered allele frequency was ≥5%, altered allele count was ≥10%, and read depth was ≥100 reads. Additionally, detected variants with an allele frequency above 2% with a minimum sequencing depth of ≥100 reads and a minimum variant depth of 10% were kept and listed if they were found in ClinVar or COSMIC databases, or had a deleterious impact on the protein. All detected variants were reviewed manually with the use of Archer Analysis software.
The relative gene expression level was assessed based on the ratio of averaged unique RNA reads originating from all GSP2 primers across the targeted and housekeeping genes (CHMP2A, GPI, RAB7A, and VCP) with Archer Analysis software.

Fluorescence In Situ Hybridization (FISH)
FGFR1 gene amplification using FISH was determined in 18 FFPE tumor samples using in a ZytoLight SPEC FGFR1/CEN 8 Dual Color Probe (containing probes specific for the 8p11 locus and the chromosome 8 centromere (CEN8)) (ZytoVision, Bremerhaven, Germany) and ZytoLight FISH-Tissue Implementation Kit (ZytoVision) according to the manufacturer's instructions. Briefly, after pre-treatment, the slides were denatured in the presence of 10 µL of probe for 10 min at 76 • C and hybridized at 37 • C overnight. Sixty tumor cell nuclei were assessed by two independent observers. The criteria of FGFR1 amplification were as follows: FGFR1/CEN8 ≥ 2.0 or the average number of FGFR1 signals per cell ≥6 or ≥10% of tumor cells containing ≥ 15 FGFR1 signals or large clusters [75].

Statistical Analysis
Comparison of the FGFR1, FGFR2, FGFR3, and FGFR4 gene expression between tumor and tumor-adjacent tissues was performed with the Mann-Whitney U test. Associations between FGFR1, FGFR2, FGFR3, and FGFR4 gene expression and clinicopathological data (Table 1), FGFR1-3 protein expression, and FGFR1 amplification were analyzed with the Kruskal-Wallis test. To estimate the association between FGFR mRNA expression level and clinical endpoint (DFS), univariate Cox proportional hazards model, the Kaplan-Meier method and log-rank test were used. Statistical calculation of correlation and strength of the relationship between mRNA expression levels from RT-PCR and NGS and FGFR protein expression were performed using Spearman's rank correlation coefficient. mRNA expression was analyzed as a continuous variable and as a categorical variable for Kaplan-Meier analysis (the median value of expression for the entire group was used as a cut-off point).
A p-value < 0.05 was considered significant. All calculations were performed using Statistica software (StatSoft, Tulsa, OK, USA).

Conclusions
FGFR1 and FGFR4 mRNA levels are significantly decreased in Sq-NSCLC tissue compared with tumor-adjacent normal tissue. Furthermore, our study shows that the increased tumor mRNA expression of FGFR3 is an unfavorable prognostic factor in terms of the risk of recurrence for Sq-NSCLC patients and the increased FGFR4 mRNA level is correlated with lymph node metastasis occurrence. We also confirm the association of increased FGFR1 mRNA with protein overexpression but not with FGFR1 amplification. Moreover, NGS revealed new and well-known FGFR variants and fusions.
FGFR mRNA and protein expression analysis in tumor and tumor-adjacent normal tissues, along with the identification of fusions and variants and investigation of amplification status, have increased our knowledge of the molecular background of Sq-NSCLC. Our data also show that the use of different methods increases the detection of FGFR aberrations, which may aid in the selection of patients most likely to respond to treatment with FGFRis.

Informed Consent Statement:
The institutional review board approved this study. No additional patient informed consent that was specific to this study was required given its retrospective nature.

Data Availability Statement: Not applicable.
Acknowledgments: The authors would like to thank Krystyna Sierota for excellent technical assistance.

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.