Pancreatic cancer ranks fourth for cancer-related deaths in western countries and is projected to become the second by 2030 (1, 2). It is a very deadly disease with the mortality rate closely approaching to the incidence rate. The median survival is less than 18 months, and the 5-year survival rate remains as low as 3 ~ 15% (3–5). The poor prognosis is mainly due to the late onset of symptoms, diagnosis at an advanced stage and subsequent rapid progression. A comprehensive identification of the risk factors can be instrumental to a better understanding of the disease etiology and to the development of methods for risk stratification, that in turn could facilitate early detection, which at the moment remains elusive.

Genetic factors play an important role in the etiology of pancreatic cancer (6). Genome-wide association studies (GWAS) have identified various frequent genetic variants associated with pancreatic cancer risk. The two largest pancreatic cancer GWAS done in European populations are the Pancreatic Cancer Cohort Consortium (PanScan) and the Pancreatic Cancer Case-Control Consortium (PanC4), and a total of 21 susceptibility loci associated at genome-wide significance level have been discovered, and studied individually and in combination (7–15). However, the identified SNPs explain only 4.1% of the total phenotypic variance of pancreatic cancer, which do not fully account for the overall 21.2% estimated genetic heritability (16). This can be explained by the relatively small effect sizes of the individual risk loci, and by the strict multiple testing correction required for GWAS (typically p < 5x10^-8), which is likely to result in a large number of false negatives.

Over the past decade, GWAS have achieved substantial success in discovering many common variants underlying the genetic architecture of complex diseases (17), including pancreatic cancer. Standard models for implying specific relationships between genotypes and phenotypes include additive, recessive and dominant models (18). The association of biallelic single nucleotide polymorphisms (SNPs) having alleles A/a with a given endpoint (e.g. disease risk) is typically analyzed with a logistic regression model logit(P) = α + β (X), where in an additive model X = 0, 1 or 2 depending on the genotype (homozygotes A/A, heterozygotes A/a and homozygotes a/a, respectively), thus the risk of disease is increased exp(β)-fold for subjects with genotype A/a and exp(2β)-fold for subjects with genotype a/a. A recessive model compares rare homozygotes a/a (who will have X=1) versus the rest (combining heterozygotes A/a and common homozygotes A/A, who will have X=0); a dominant model compares A/A (X=0) versus A/a + a/a (X=1). As most GWAS studies assume that allelic effects are additive, most of the associations reported in GWAS consider only the additive model of inheritance. But for variants which do not follow an intermediate model of inheritance, the recessive or the dominant genetic model can have more power to detect associations. Reanalysis of GWAS data with the recessive model of inheritance, considering homozygotes for the minor allele as the only “exposed” category could help to identify additional risk loci for non-negligible subsets of SNPs (19).

To discover novel recessively inherited pancreatic cancer risk loci, we performed a secondary analysis using genotyping data from all published pancreatic cancer GWAS conducted in subjects of European origin, and then replicated the most promising variants in cases and controls collected from the PANcreatic Disease ReseArch (PANDoRA) consortium. Better understanding the genetic background of the disease could be an invaluable tool to stratify the population by individual risk and increase our chances of early detection.

In the first phase, which was conducted at a genome-wide scale, we re-analyzed the data from the PanScan+PanC4 GWAS dataset according to a recessive model of inheritance, and we observed six SNPs that showed specific recessive associations with PDAC risk with p<10^-5 while p>10^-3 using the additive model ( Figure 1 , Supplementary Table 2 ). The correlated SNPs in these regions (r²>0.8 in LD) did not show evidence of stronger association under an additive model ( Supplementary Table 3 , Supplementary Figure 1 ). The associations of the genotyped SNPs with PDAC risk under the additive and recessive genetic models are shown in Figure 1 ( Supplementary Table 2 ). In the validation phase in PANDoRA, no statistically significant associations (p<0.05) were observed, using the recessive model, except for rs2066357. However, this SNP showed high heterogeneity, with an opposite effect compared to the discovery phase under the recessive genetic model.

In meta-analyses, none of the six SNPs reached the conventional genome-wide significance threshold (p<5×10⁻⁸). However, MIR96 rs4626538 (OR=0.93; p=4.42×10⁻⁶), RP5-991O23.1 rs7008921 (OR=1.42; p=9.68×10⁻⁶) and ARHGAP27 rs147904962 (OR=0.75; p=4.08×10⁻⁶) maintained a specific recessive effect compared to the additively inherited effects (p=1.06×10⁻³, p=0.05 and p=3.97×10⁻³, respectively), and the p-values of rs4626538 and rs147904962 in the meta-analysis were slightly lower in comparison with those observed in the first phase. The results after multiple imputation were generally consistent with those without multiple imputation ( Supplementary Table 2 ). Results did not change when we added the PANDoRA cases and controls from Brazil, who are ethnically admixed ( Supplementary Tables 4, 5 ).

We used data from the GTEx consortium to investigate associations between genetic variants and RNA expression. We observed that the rs147904962-A allele was associated with increased LRRC37A4P RNA expression in adipose tissue (p=8.1×10⁻⁶). An expanded list of linked SNPs (in LD with our six candidate SNPs, r²>0.6) was also considered for the GTEx analysis; we found that the T allele of rs590097 (in LD with the A allele of rs36018702, r^{2 =} 0.74, D’=1) was associated with higher expression of BCL2L11 in pancreas (p=5.64×10⁻⁶). No expression quantitative trait loci (eQTL) associations in pancreas were found for the other SNPs. Haploreg and RegulomeDB did not show evidence for functional effect for these variants.

Using a threshold of >2 for distance-normalized chromatin interaction frequency, 3DIV predicted C1orf21 and APOBEC4 to be interaction genes for rs1339571, BUB1 for rs36018702, MIR4423 for rs2066357, SPPL2C, SLC4A1, RUNDC3A, LOC100133991, TEX34, ITGA2B, and C17orf57 for rs147904962, respectively.

Additional analyses with SNPnexus and OpenTargets Genetics did not suggest any clear functional link between our candidate SNPs and pancreatic physiology or pathology. Likewise, when we reanalyzed with FUMA the results of the GWAS analysis according to the recessive model, we did not observe any noteworthy signal in the regions of the six candidate SNPs.

The gene-based analysis using MAGMA based on the p-values of the recessive model revealed that 14 genes were associated with PDAC risk at p < 0.001 ( Supplementary Table 6 ). Two of these genes showed evidence for association at p < 0.001 under the recessive model (CTSG 14q12, p = 2.53x10^-4; LEPROTL1 8p12, p = 4.34x10^-4), but not with the additive one (p = 0.20 and p = 0.10, respectively). Then we verified the expression level of the two genes in pancreatic cancer patients using GEPIA2. We found that LEPROTL1 has increased expression in pancreatic cancer tissues compared to adjacent normal pancreatic tissues of the same patients ( Figure 2 ).

GWAS data are usually analyzed according to an additive genetic model, which is generally considered to be a good surrogate for other genetic models, except for the recessive one (19). Researchers have reported risk variants that showed specifically stronger evidence under a recessive model than an additive model, for type 2 diabetes (36), schizophrenia (37), high triglycerides (38), and other traits (39, 40), but not for PDAC yet. To identify recessive susceptibility loci for PDAC risk, we performed a secondary analysis with the largest currently available pancreatic cancer GWAS datasets (PanScan and PanC4) of European ancestry and attempted the replication of the six most promising variants in additional samples collected from the PANDoRA consortium, with a combined sample size of 11,981 PDAC cases and 10,525 controls. In this study, none of our results reached genome-wide statistical significance (p<5x10⁻⁸) in either phase, or in the meta-analyses, therefore our results are not conclusive. However, for five of the six selected SNPs the results of the meta-analysis do not exclude the possible recessive association with pancreatic cancer risk. In particular, rs4626538 (7q32.2), rs7008921 (8p23.2) and rs147904962 (17q21.31) maintained a large difference in significance between recessive effects compared with the additively inherited effects.

None of the previous studies indicated a link between these loci and pancreatic cancer risk. No variants in high LD (r²>0.8) have been previously associated with any trait or disease in GWAS, although variants in low to moderate LD (r^{2 =} 0.14~0.60, D’=0.88~1 in Europeans) with rs147904962 have been reported to be associated with waist-to-hip ratio and with risk of developing allergic diseases. The minor G allele of rs7214661 (r^{2 =} 0.19, D’=0.98) was associated with higher risk of allergic disease (41) while the corresponding A allele of rs147904962 was associated with lower risk of pancreatic cancer in our study. It is consistent with the protective effect of allergy for pancreatic cancer in epidemiologic studies (42).

Additionally, GTEx showed that rs590097 regulates BCL2L11 expression in pancreas tissue. BCL2L11 is a member of the BCL2 family and plays a role in neuronal and lymphocyte apoptosis. There is evidence shown that BCL2L11 is one of the major genes contributing to apoptosis, known to be important for pancreatic biology (http://www.genome.jp/kegg/pathway.html) (43). Moreover, the observed association that rs36018702-A (correlated with rs590097-T) showed increased risk of PDAC is consistent with the higher expression of BCL2L11 in pancreatic cancer tissues than in normal pancreas tissues found through GEPIA2.

BUB1 is the interaction gene of rs36018702 predicted by 3DIV. There is evidence that BUB1 is overexpressed in PDAC tissues, suggesting a role of BUB1 in PDAC progression, and therefore corroborating the association of rs36018702 and PDAC risk (44).

The lowest p-value we observed in the meta-analysis is 4.08x10^-6 for the association of rs147904962 (17q21.31) with the risk of PDAC. rs147904962 is situated 17kb at the 3’ end of Rho GTPase Activating Protein 27 (ARHGAP27). This gene encodes a member of a large family of proteins that activate Rho-type guanosine triphosphate (GTP) metabolizing enzymes and are involved in cancer through the dysregulation of this mechanism. As ARHGAP27 mRNA is expressed in pancreatic cancer, we speculate that rs147904962 mediates regulation of cancer-associated ARHGAP27, promoting carcinogenesis through dysregulation of Rho/Rac/Cdc42-like GTPases (45). However, it has to be acknowledged that this SNP is not known to be located in a regulatory region of ARHGAP27.

Gene-based analyses based on the PanScan and PanC4 datasets (we were not able to replicate these analyses in PANDoRA, which does not have GWAS data) showed that SNPs in LEPROTL1 and CTSG were associated with PDAC risk according to the recessive, but not to the additive model. The bioinformatic analysis identified that LEPROTL1 was highly expressed in pancreatic cancer compared to matched normal pancreatic tissue of the same patients, suggesting a potential involvement in the etiopathology of PDAC. The leptin receptor overlapping transcript-like 1 gene (LEPROTL1) encodes a membrane protein, and may play a role in liver resistance by suppressing the growth hormone activity (46, 47), while the pancreatic cancer-related functions of LEPROTL1 remain unknown. The cathepsin G gene (CTSG) encodes a neutrophil serine protease of the chymotrypsin family, which was shown to affect neutrophil infiltration into the pancreas in a mouse model of pancreatitis (48). Based on this circumstantial evidence it is tempting to speculate a role for this gene and its polymorphisms in modulation of inflammation in the pancreas, which plays a role in the etiology of PDAC. However, to the best of our knowledge, a role for CTSG in pancreatic cancer has not been reported in the literature.

The lack of direct functional evidence for the SNPs of interest from bioinformatic analyses may at least in part reflect the fact that also bioinformatic tools/databases have not been designed to address effects of real recessive alleles. Ad hoc tools are needed to better understand the genetic architecture of complex genetic diseases.

It is hard to reach sufficient statistical power to detect variants with recessive effects, unless they are very frequent or have very large effects. Given the effective combined sample size of 11,981 PDAC cases and 10,525 controls, disease prevalence of 1.6%, and a significance cut-off of p<5×10⁻⁸, we had at least 80% power to detect a association with ORs equal to those observed in the discovery phase for the rare homozygote genotype for SNPs rs7008921 and rs147904962, whereas for the other SNPs power ranged between 54% to 69%. Thus, our study, in spite of the large sample size, lacked statistical power to confirm the risk with recessive model for some of the SNPs. It is worth noting that between PanScan, PanC4 and PANDoRA we have used the largest available resources for genetics of pancreatic cancer in populations of European origin. Our hypothesis that some variants may be associated with pancreatic cancer risk with a recessive model of inheritance was not disproved, but to prove it convincingly will require even larger datasets that will become available as more GWAS on pancreatic cancer risk are performed.

Identifying high-risk groups could contribute to focus surveillance and invasive screening measures, thereby improving the chance of early detection. Polygenic risk scores (PRS) approaches which could combine modest effect from each risk SNPs have demonstrated accuracies between 59% and 63% for predicting the risk of PDAC when including both non-genetic and genetic factors (14, 49–51). The accuracy of the existing PRS is not ready yet to be used in the clinical practice. It is necessary to expand the PRS with additional risk factors to improve its predictive power. For example, PRS including more SNPs that are not genome-wide significant but having noteworthy effects such as the ones we highlighted in this work may provide an additive contribution to the overall performance.

In conclusion, we propose some candidate SNPs as recessively inherited genetic variants for pancreatic cancer risk in European populations, which should be further confirmed by better powered investigations and/or meta-analysis of our results with those of other studies. Although none of the SNPs reached the genome-wide statistical significance, it is still worth to include these relevant SNPs into the PRS approach for risk stratification. A risk stratification approach with high predictive power could be used to identify subgroups at particularly increased risk of pancreatic cancer, either in the general population or in groups that are already known to have an elevated risk, such as diabetics.

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Materials and Methods, Study Populations section.

The studies involving human participants were reviewed and approved by the Ethics Commission of the Medical Faculty of the University of Heidelberg. The patients/participants provided their written informed consent to participate in this study.

DC and FC conceived and designed the study. YL performed the lab work. YL performed data curation and analysis. YL drafted the manuscript. DC, MGe, and FC reviewed and edited the manuscript. All other authors provided samples and data. All authors critically read, commented and approved the final manuscript.

This work was supported by intramural funding of German Cancer Research Center (DKFZ); and by Fondazione Tizzi (www.fondazionetizzi.it); Fondazione Arpa (www.fondazionearpa.it); the Economic Development and Innovation Operative Programme Grant (Grant number GINOP-2.3.2-15-2016-00048); the Human Resources Development Operational Programme Grant (Grant number EFOP 3.6.2‐16‐2017‐0006); Associazione Italiana Ricerca Cancro (Grant number 5x1000, 12182, and IG 17177); Fondazione Italiana Malattie Pancreas – Ministero Salute (Grant number FIMPCUP_J38D19000690001); Fondazione Cariverona: Oncology Biobank Project “Antonio Schiavi” (Grant number 203885/2017); the Ministry of Health of the Czech Republic (Grant number NV19-03-00097, FNOL-00098892); student grant from Palacky University (Grant number IGA_LF_2020_005). The ESTHER study was funded by the Baden-Württemberg State Ministry of Science, Research and Arts (Stuttgart, Germany).

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.