Condemned or Not to Die? Gene Polymorphisms Associated With Cell Death in Pemphigus Foliaceus.

Pemphigus foliaceus (PF) is an autoimmune blistering skin disease that occurs sporadically across the globe and is endemic in Brazil. Keratinocyte adhesion loss (acantholysis) is associated with high levels of anti-desmoglein 1 IgG autoantibodies, but the role of cell death is poorly understood in PF. Current evidence disqualifies apoptosis as the major cell death mechanism and no other process has yet been investigated. To approach the role of variation in genes responsible for cell death pathways in pemphigus susceptibility, we systematically investigated the frequencies of 1,167 polymorphisms from genes encoding products of all 12 well-established cell death cascades (intrinsic and extrinsic apoptosis, necrosis, necroptosis, ferroptosis, pyroptosis, parthanatos, entotic, NETotic, lysosome-dependent, autophagy-dependent, and immunogenic). By multivariate logistic regression, we compared allelic and genotypic frequencies of 227 PF patients and 194 controls obtained by microarray hybridization. We found 10 variants associated with PF (p < 0.005), belonging to six cell death pathways: apoptosis (TNF, TRAF2, CD36, and PAK2), immunogenic cell death (EIF2AK3, CD47, and SIRPA), necroptosis (TNF and TRAF2), necrosis (RAPGEF3), parthanatos (HK1), and pyroptosis (PRKN). Five polymorphisms were associated with susceptibility: TNF rs1800630*A (OR = 1.9, p = 0.0003), CD36 rs4112274*T (OR = 2.14, p = 0.0015), CD47 rs12695175*G (OR = 1.77, p = 0.0043), SIRPA rs6075340*A/A (OR = 2.75, p = 0.0009), and HK1 rs7072268*T (OR = 1.48, p = 0.0045). Other five variants were associated with protection: TRAF2 rs10781522*G (OR = 0.64, p = 0.0014), PAK2 rs9325377*A/A (OR = 0.48, p = 0.0023), EIF2AK3 rs10167879*T (OR = 0.48, p = 0.0007), RAPGEF3 rs10747521*A/A (OR = 0.42, p = 0.0040), and PRKN rs9355950*C (OR = 0.57, p = 0.0004). Through functional annotation, we found that all associated alleles, with the exception of PRKN rs9355950*C, were previously associated with differential gene expression levels in healthy individuals (mostly in skin and peripheral blood). Further functional validation of these genetic associations may contribute to the understanding of PF etiology and to the development of new drugs and therapeutic regimens for the disease.

Introduction

Pemphigus foliaceus (PF) is an autoimmune bullous disease of the skin, characterized by the production of autoantibodies that recognize the desmosome protein desmoglein 1 (DSG1) (1, 2). The binding of antibodies to this cell adhesion molecule is accompanied by acantholysis (keratinocyte detachment) and lesions in the superficial granular layer of the epidermis. In contrast, a different form of pemphigus called pemphigus vulgaris (PV) is characterized by anti-desmoglein 3 antibodies, which bind in deeper layers of the skin and cause blisters in the skin and mucosa (3). PF occurs sporadically around the world−0.75–5 cases/million per year (4). In Brazil, PF is endemic and commonly known as fogo selvagem (“wild fire” in Portuguese). The prevalence of PF in Limão Verde, located in the Brazilian state of Mato Grosso do Sul, is of 3.04%, one of the highest prevalences ever reported for autoimmune diseases (5). PF pathogenesis is multifactorial, resulting from poorly understood interactions between multiple environmental and genetic factors (6, 7).

Cell Death in Pemphigus: An Unsolved Issue

Despite the pathogenic relevance of IgG autoantibodies in the acantholytic process, the exact mechanisms that lead keratinocytes to death remain unknown (8). Apoptosis has been suggested to play an important role in some dermatoses with positive Nikolsky's sign (skin detachment after slight rubbing) as in PV and PF (9–12), either prior (10, 13–15), or after the loss of cell adhesion (8, 16–18). As early as 1998, Gniadecki et al. reported many apoptotic keratinocytes in acantholytic tissue and in the cohesive epidermis just under the blisters of sporadic PF and PV lesional skin biopsies, as judged by positive TUNEL signs (terminal deoxynucleotidyl transferase dUTP nick end labeling) (9). Rodrigues et al. (19) also found TUNEL-positive keratinocytes in 12/13 biopsies of perilesional skin of endemic PF lesions. Among them, 10/13 presented keratinocytes with intense expression of proapoptotic inducible nitric oxide synthase (iNOS) and 3/13, with rather discrete-moderate expression of proapoptotic FAS protein. Anti-apoptotic Bcl-2 occurred in 4/13 biopsies only, being much more abundant in the inflammatory infiltrate, which also had discrete-moderate expression of interleukin 1, interferon gamma, and tumor necrosis factor alpha (TNF-α) proinflammatory cytokines (11/13) (19). After the passive transfer of sporadic PF-antibodies in the experimental neonatal mouse model, keratinocytes expressed the proapoptotic Bax factor, followed by activation of caspases (CASP) 3 and 6, and down-regulation of the anti-apoptotic Bcl-x(L) factor. In this model, apoptotic inhibitors abrogated the acantholytic process (14). Nevertheless, p38MAPK signaling occurred in this same model in two phases, and the first peak of activation coincided with acantholysis, prior to the second peak that induced activation of CASP-3 (18). Taking into account ultra-violet (UV) light exposure as a known risk factor for endemic PF, it is interesting that caspases-3 and -7 cleave desmoglein-1 intracellularly and that knock-down of desmoglein-1 protects cells from UV induced apoptosis in irradiated keratinocytes (20).

On the other hand, only few apoptotic cells were detected in skin biopsies of endemic PF patients, whereas p63 marked many undifferentiated cells distributed over the whole epidermis, both in injured and non-injured skin (21). Electron microscopy did not reveal any morphological signs of apoptosis—retraction of pseudopods, pyknosis, karyorrhexis, and plasma membrane blebbing—in acantholytic tissue of PV and PF patients (8, 22). There were, as well, no signs for activation of caspases (as cleaved CASP3 and CASP8, fractin, or nuclear poly (ADP-ribose) polymerase—PARP) in PV and PF biopsies, nor in PV or PF IgG–incubated healthy breast reduction skin biopsies (8). A possibility suggested by Schmidt and Waschke (16) is that caspase signaling adds in desmosome destabilization, but that apoptosis itself is not responsible nor needed for acantholysis to occur.

Thus, whereas some authors state that apoptosis causes cell death in PF (9, 19, 23), others found no clear evidence of such event (8, 17, 18, 22). The uncertainty about how cell death takes place in PF, as well as the scarcity of genetic association studies on this topic, prompted the current investigation encompassing genes of all known cell death routes. In fact, there are several pathways orchestrating cell death, classified following morphological, biochemical, and functional features. In 2018, the Nomenclature Committee on Cell Death (NCCD) proposed 12 pathways orchestrating cell death, supported by genetic, biochemical, and functional results: intrinsic apoptosis, extrinsic apoptosis, mitochondrial permeability transition (MPT)-driven necrosis, necroptosis, ferroptosis, pyroptosis, parthanatos, entotic, NETotic, lysosome-dependent, autophagy-dependent, and immunogenic pathways (24). All of them are classified as regulated cell death (RCD) routes. All routes depend on the molecular machinery (causing the activation of one or more signal transduction pathways), which can be pharmacologically and/or genetically modulated. RCD begins with excessive cellular stress, caused by non-coped perturbations of the intra- and extracellular environment (24, 25). Given the poorly understood underlying mechanisms leading to keratinocyte death in PF, we intended to identify genetic variants associated with PF analyzing variants from genes whose products are known to be directly involved in RCD routes.

Materials and Methods

Population Sample

A total of 227 PF patients and 194 unrelated controls with no diagnosis or familial history of autoimmune illnesses were analyzed in this study. Patients received clinical and/or immunohistological diagnosis. Patients and controls were recruited from 1984 to 2015 in South/Southeastern/Central-Western Brazilian hospitals: Hospital Adventista do Pênfigo (Campo Grande, Mato Grosso do Sul), Lar da Caridade—Hospital do Fogo Selvagem (Uberaba, Minas Gerais), Hospital das Clínicas—University of São Paulo (Ribeirão Preto, São Paulo), Hospital de Clínicas—Federal University of Paraná, Hospital de Dermatologia Sanitária São Roque, and Hospital Santa Casa de Misericórdia (Curitiba, Paraná). All individuals enrolled in this study were unrelated, predominantly of European ancestry and living in rural endemic areas. No history of other autoimmune diseases was reported for patients, as well as no history of any autoimmune disease for the controls. The median age was 40.9 years for patients (minimum 6, maximum 83) and 44.8 years for controls (minimum 11, maximum 86). In both groups, 52% were women. Peripheral blood was collected, from which DNA was extracted by the phenol-chloroform-isoamyl alcohol protocol (26). Patients and controls voluntarily agreed to participate in the study and gave their informed consent. This study was carried out according to the Declaration of Helsinki, with the approval from the Brazilian National Ethics Committee (CONEP) under the protocol number CAAE 02727412.4.0000.0096 and approval number 505.988.

Selection of Candidate Genes and Genotyping

In agreement with NCCD, we selected genes encoding proteins involved in essential aspects of at least one of the following cell death cascade processes: intrinsic apoptosis, extrinsic apoptosis, mitochondrial permeability transition (MPT)-driven necrosis, necroptosis, ferroptosis, pyroptosis, parthanatos, entotic, NETotic, lysosome-dependent, autophagy-dependent, and immunogenic pathways (24). We identified the genomic positions of each gene, considering regulatory regions of one thousand base pairs upstream and downstream from the transcription start and end sites, respectively, of the longest transcript, according to the GRCh37/hg19 human genome version (available at: http://www.lgmh.ufpr.br/data/Supplementary_material_1_Bumiller-Bini_2019.xlsx).

Genotyping was performed using microarray hybridization (CoreExome-24 v1.1 Illumina®) in 194 and 227 DNA samples from controls and endemic PF patients, respectively. A total of 3,798 SNPs were extracted from DNA microarray data, filtered based on genotyping quality, and on population genetics criteria: excluding those SNPs with minor allele frequencies <1%, with genotypic distributions deviating from those expected by Hardy-Weinberg equilibrium in controls (p < 0.05) and/or with high linkage disequilibrium (r2 ≥ 0.8). After filtering, a total of 1,167 SNPs remained for subsequent analyses (Figure 1, set available at: http://www.lgmh.ufpr.br/data/Supplementary_material_2_Bumiller-Bini_2019.xlsx).

Figure 1

Representation of the distribution of all 1,167 SNPs used in this study according to their related cell death pathway(s). Traced lines were used to represent additional shared SNPs (unions) between two pathways that were not otherwise graphically united. Genes from the intrinsic and extrinsic apoptosis pathways were represented together. Unrepresented SNPs, shared among pathways: one between apoptosis, immunogenic, lysosome, and necrosis pathways; two between apoptosis, immunogenic, lysosome, and pyroptosis pathways; two between necroptosis, necrosis, pyroptosis, NETotic, and immunogenic pathways; three between pyroptosis and necroptosis; three between apoptosis, necroptosis, necrosis, and NETotic pathways; four between apoptosis, necrosis, and immunogenic pathways; and four between apoptosis, entotic, and immunogenic pathways).

Association Analysis

Association analysis was carried out by binary logistic regression, using sex and two principal components (PCA) as co-variables. The PCA eliminates spurious associations due to possible population stratification. The significance level was set to p < 0.005 (27). The analyses were performed using PLINK software version 1.1.9 (28). We applied Fisher's exact test to perform the stratified association analysis (29), in this case, with a significance level set to p < 0.05.

In silico Analysis

We performed in silico analysis to explore the genetic associations identified in this study. Linkage disequilibrium within the Iberian population (major ancestral of the Euro-Brazilian population) was evaluated with LDlink (30). The predicted regulatory impact of the genetic variants was verified using ENSEMBL (31), GTEx portal (32), UCSC (https://genome.ucsc.edu/cgi-bin/hgGateway) (33), HaploReg (https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php) (34), Innatedb (https://www.innatedb.com) (35), and Blood eQTL (expression quantitative trait loci) browser (36), which assemble public datasets, and published data. GTEx portal and the Blood eQTL browser inform whether a SNP is an expression quantitative trait loci (eQTL). Innatedb informs whether there is a physical interaction between proteins.

Results

We found 10 SNPs associated with PF (p < 0.005) located in 10 different genes participating in six RCD routes: apoptosis (CD36, PAK2, TNF, and TRAF2), immunogenic cell death (CD47, EIF2AK3, and SIRPA), necroptosis (TNF and TRAF2), necrosis (RAPGEF3), parthanatos (HK1), and pyroptosis (PRKN) (Table 1). Altogether, the less frequent alleles of five SNPs were associated with PF susceptibility, while the less frequent alleles from the other five SNPs, with protection. All associated variants were located within non-protein coding regions.

Table 1

Cell death-related gene variants associated with PF.

Gene	SNP	MAF (%)		Model	CONTR	CASE	OR	95%CI	p
		Contr	Case
APOPTOSIS AND NECROPTOSIS
TNF	rs1800630	15.7	26.2	add	61/327	117/329	1.90	[1.34–2.70]	0.0003
6p21.33	C>a			rec	7/187	14/209	1.79	[0.70–4.54]	0.2214
	Promoter			dom	54/140	103/120	2.24	[1.49–3.38]	0.0001
TRAF2	rs10781522	49.2	38	add	191/197	171/279	0.64	[0.48–0.84]	0.0014
9q34.3	g>A			rec	47/147	34/191	0.58	[0.35–0.95]	0.0305
	Intron 9			dom	144/50	137/88	0.52	[0.34–0.79]	0.0024
APOPTOSIS
PAK2	rs9325377	49.5	42	add	196/192	190/162	0.72	[0.54–0.94]	0.0173
3q29	a>G			rec	55/139	36/190	0.47	[0.29–0.77]	0.0023
	Intron 1			dom	141/53	154/72	0.81	[0.53–1.25]	0.3449
CD36	rs4112274	7.51	15.27	add	29/357	69/383	2.14	[1.33–3.44]	0.0015
7q21.11	C>t			rec	1/192	5/221	3.92	[0.45–34.2]	0.2169
	Intron 3			dom	28/165	64/162	2.23	[1.34–3.68]	0.0018
IMMUNOGENIC
EIF2AK3	rs10167879	17.6	9.9	add	67/313	43/393	0.48	[0.31–0.73]	0.0007
2p11.2	C>t			rec	5/185	4/214	0.62	[0.16–2.40]	0.4940
	Intron 14			dom	62/128	39/179	0.41	[0.26–0.66]	0.0002
SIRPA	rs6075340	33.5	42.1	add	130/258	191/263	1.50	[1.13–2.00]	0.0055
20p13	a>G			rec	17/177	47/180	2.75	[1.51–4.98]	0.0009
	Intron 2			dom	113/81	144/83	1.35	[0.90–2.02]	0.1449
CD47	rs12695175	11.4	17.9	add	44/342	81/371	1.77	[1.20–2.63]	0.0043
3q13.12	T>g			rec	4/189	10/216	2.23	[0.68–7.30]	0.1850
	Intron 6			dom	40/153	71/155	1.95	[1.23–3.07]	0.0041
NECROSIS
RAPGEF3	rs10747521	40.1	36.7	add	155/231	165/285	0.81	[0.61–1.1]	0.1676
12q13.11	a>G			rec	36/157	23/202	0.42	[0.23–0.76]	0.0040
	Intron 1			dom	155/231	165/285	1.03	[0.20–0.69]	0.8837
PARTHANATOS
HK1	rs7072268	43.0	53.1	add	167/221	241/213	1.48	[1.13–1.94]	0.0045
10q22.1	t>C			rec	37/157	70/157	1.57	[1.02–2.42]	0.0412
	Intron 5			dom	130/64	171/56	1.87	[1.18–2.95]	0.0076
PYROPTOSIS
PRKN	rs9355950	34.3	22.9	add	133/255	104/350	0.57	[0.42–0.78]	0.0004
6q26	T>c			rec	25/169	11/216	0.35	[0.16–0.73]	0.00502
	Intron 4			dom	108/86	93/134	0.55	[0.37–0.81]	0.0024

Logistic regression association tests were done with allele frequencies (“add”), frequency of homozygotes for the minor allele (“rec”- recessive model), and summed frequencies of heterozygotes and homozygotes for the minor allele (“dom”—dominant model). The minor alleles in our sample are given in lowercase; they are the reference for the associations. In bold, significant associations (p < 0.005). SNP, single nucleotide polymorphism; MAF, minor allele frequency; CONTR, controls; CASE, cases; Model, association tests; OR, odds ratio; CI, confidence interval; TNF, tumor necrosis factor alpha; TRAF2, TNF receptor associated factor 2; PAK2, p21 (RAC1) activated kinase 2; CD36, CD36 molecule; EIF2AK3, eukaryotic translation initiation factor 2 alpha kinase 3 (also known as PERK); SIRPA, signal regulatory protein alpha; CD47, CD47 molecule; RAPGEF3, Rap guanine nucleotide exchange factor 3 (also known as EPAC); HK1, hexokinase 1; PRKN, parkin RBR E3 ubiquitin protein ligase (also known as PARK2).

Functional Annotation in silico Analysis

To explore the potential effects of all 10 genetic variants associated with PF, we used functional annotation available in reference public databases. As outlined below, most of the associated variants are associated with gene expression in different tissues (Table 2). Nevertheless, the association of these variants (with disease and with gene expression) may also be explained by strong linkage disequilibrium with other causal variants (hitch-hiking effect).

Table 2

Cell death-related gene variants associated with PF.

Variant	Allele associated with PF	eQTL effect (25, 29)	Tissues
TNF_rs1800630	A risk	low expression: TNF (p = 3.71e−21) and LTA (p = 2.88e−24) (36).high expression: DDX39B (p = 2.6e−4), CSNK2B (p = 3.14e−11), HCP5 (p = 2.5e−3) and MICB (p = 2.0e−15) (36).	Peripheral blood.
TRAF2_rs10781522	G prot	low expression: C8G (p = 2.6e−5), AGPAT2 (p = 9.1e−4), and CLIC3 (p = 4.79e−17) (32, 36).high expression: TRAF2 (p = 2.73e−4), PHPT1 (p = 1.2e−6)	Skin, peripheral blood, muscle.
PAK2_rs9325377	A prot	low expression: PAK2 (p = 1.2e−6) (32).high expression: PIGX (p = 9.7e−7). PIGX maps just upstream of PAK2, on the same strand of chromosome region 3q29 (32).	Skin, spleen.
CD36_rs4112274	T risk	low expression: CD36 (p = 2.37e−29) (36).	Peripheral blood.
EIF2AK3_rs10167879	T prot	low expression: AC062029.1 (p = 3.0e−4) (32)high expression: EIF2AK3 (p = 3.5e−4), ANKRD36BP2 (p = 1.6e−6) (32).	Skin, whole blood.
SIRPA_rs6075340	A risk	low expression: SIRPA (p = 4.1e−7) (32).high expression: SIRPB1 (p = 6.6e−98) (32).	Esophagus, skin.
CD47_rs12695175	G risk	low expression: CD47 (p = 3.37e−4) (36).high expression: MYH15 (p = 1.3e−8) (32).	Peripheral blood, skin.
RAPGEF3_rs10747521	A prot	low expression: RAPGEF3 (p = 2.0 e−8) (32).high expression: PCED1B-AS1 (32).	Lung, heart.
HK1_rs7072268	T risk	high expression: HK1 (p = 2.3 e−6) (32).	Tibial nerve.

PRKN_ rs9355950 is not an eQTL.

eQTL, expression quantitative trait loci; LTA, lymphotoxin alpha; DDX39B, dead box polypeptide 39 B; CSNK2B, casein kinase II beta; HCP5, HLA complex P5; MICB, major histocompatibility complex class I-chain related gene B; PHPT1, phosphohistidine phosphatase 1; C8G, complement C8 gamma chain; AGPAT2, 1-acylglycerol-3-phosphate O-acyltransferase 2; CLIC3, chloride intracellular channel 3; PIGX, phosphatidylinositol glycan anchor biosynthesis class X; ANKRD36BP2, ankyrin repeat domain 36B pseudogene 2; SIRPB1, signal regulatory protein beta 1; MYH15, myosin heavy chain 15.

In accordance with functional annotation, TNF _rs1800630, TRAF2_rs10781522, PAK2_rs9325377, EIF2AK3_rs10167879, SIRPA_rs6075340, PRKN_rs9355950, and HK1_rs7072268 SNPs are located within predicted transcription factor binding sites. On the other hand, the location of TNF_rs1800630, TRAF2_rs10781522, PAK2_rs9325377, SIRPA_rs6075340, CD47_rs12695175, CD36_rs4112274, HK1_rs7072268, and PRKN_rs9355950 SNPs overlap with enhancers and/or promoters that are important in several tissues, including skin, T and B cells from blood (enriched in H3K27Ac) (33, 34). Despite being located downstream of the last exon of LTA and far from the transcription start site of TNF, TNF_rs1800630 is located in a DNase hypersensitive region bound by RNA polymerase II subunit A, providing strong evidence for active transcription of this region in B lymphocytes (33, 34).

Furthermore, the SIRPA protein interacts physically with CD47, suggesting a possible gene interaction/synergistic effect of the associated polymorphisms of both genes on disease susceptibility (35). In fact, carriers of SIRPA_rs6075340*A and CD47_rs12695175*G alleles were more frequent among patients than controls (OR = 2.02 [95%CI = 1.08–3.81], p = 0.0202). TNF also interacts physically with the TNF receptor TRAF2, as shown in a cervical cancer cell line (35). Thus, an additive susceptibility effect was evident in carriers of TNF_rs1800630*A and TRAF2_rs10781522*A (OR = 4.78 [95%CI = 2.18–10.94], p < 0.0001) (Supplementary Tables 1, 2).

Discussion

Although the underlying molecular mechanisms of RCDs overlap considerably, our approach allowed us to suggest the possible role of non-apoptotic RCDs and raise hypotheses to explain the genetic associations that we observed.

TNF, TRAF2, and PAK2: Apoptosis and Necroptosis

The genetic associations of TNF_rs1800630*A (also known as −863*A, OR = 1.9, p = 0.0003) and TRAF2_rs10781522*G (OR = 0.64, p = 0.0014) with PF point to a specific role of the death receptor pathway in the disease. As the major proinflammatory cytokine mediating apoptosis and necroptosis, TNF binds TNF type I receptor (TNFRI), activating NF-kB through TRADD (TNFR-associated death domain), and TRAF2 recruitment. This sequence of events culminates in cell survival and inflammation (37, 38). In the absence of TRAF2, TNF binding to TNFRI builds “death-inducing signaling complexes” that can activate either necroptosis or caspase-dependent cell death. In the latter case, this leads to the activation of CASP8 and of the apoptotic cascade (37, 39). On the other hand, TRAF2 was recently reported to positively regulate caspase-2 activation, which initiates apoptosis and is a negative regulator of necroptosis (40, 41). For necroptosis to ensue, caspase inactivity or absence must prevail [thus, in the absence of TRAF2 (38)] (42). Interestingly, although TNF_rs1800630*A carriers present reduced TNF expression, TRAF2_rs10781522*G is associated with higher TRAF2 gene expression in peripheral blood (36). In contrast with TNF, the allele from TRAF2 was associated with PF protection. Since both molecules are known to interact physically (35), certain allelic combinations of these genes present an additive effect toward susceptibility to the disease.

Moreover, many associations of TNF_rs1800630*A with diseases have been reported. Those of enhanced susceptibility to cancer seem to indicate that reduced TNF levels increase the chance of inappropriate cell proliferation, due to insufficient signaling for apoptosis/necroptosis [e.g., hepatocellular carcinoma (43), non-Hodgkin lymphoma (44, 45), gastric cancer (46), and colon cancer (47)]. On the other hand, the same allele was associated with progression to and severity of infections, as for severe malaria, including cerebral malaria (48, 49), HPV-associated oral squamous carcinoma (50), HBV chronification (51), as well as with chronic disorders, which may rely on insufficient immunological response to different kinds of aggressors [steroid-induced osteonecrosis of the femoral head (52), progression of acute pancreatitis to systemic inflammation and multi-organ dysfunction syndrome (53), chronic periodontitis (54), and cardiovascular heart disease (55)].

The TNF_rs1800630*A was associated with predisposition to autoimmune disorders affecting the skin and mucosal tissue, as vitiligo (56), systemic lupus erythematosus (SLE) (57), and Crohn's disease (58). As in these autoimmune skin disorders, the rs1800630*A was associated with susceptibility to endemic PF, in the Brazilian population (OR = 1.9, p = 0.0003). Interestingly, TNF microsatellite polymorphisms were associated with susceptibility to PF in Tunisia, where the disease is also endemic (59). The rs1800630*A causes decreased TNF transcription and lower serum TNF levels (51, 60–62). Furthermore, carriers of rs1800629*A (also known as −308*A) presented higher susceptibility to pemphigus (both PV and PF) (63). This allele is associated with TNF expression levels in whole blood, similar as observed for rs1800630*A. Thus, although both occur in different haplotypes (rs1800630_rs1800629*CA and AG, with CG representing more than 70% of the allelic combinations in the Iberian population), both present the same effect on reducing gene expression and increasing PF susceptibility.

In addition, the TRAF2 PF protective allele is associated with altered expression of at least five genes implicated in cell survival or death (PHPT1, PTGDS, LCNL1, C8G, CLIC3). It is also associated with the expression of one gene related to size and inflammatory itching reaction, after mosquito bites (AGPAT2) (64–70). The AGPAT2 association is particularly interesting, given the epidemiological association of endemic PF with massive and continued exposure to mosquito bites (71, 72).

Another PF susceptibility allele, CD36_rs4112274*T, is associated with decreased CD36 gene expression in blood. This receptor is a mediator of both endoplasmic reticulum stress and the generation of reactive oxygen species in the intrinsic apoptosis pathway (73). Its predicted down-regulation in most PF patients is an additional argument favoring a less prominent role (if any) for apoptosis in the disease. Its expression in keratinocytes correlates with early wound healing (74). Thus, lower expression is also predicted to increase the extent of PF epidermal lesions.

The PAK2 rs9325377*A/A was associated with PF protection (OR = 0.48, p = 0.0023). The PAK2 product is activated through proteolytic cleavage, by caspase-mediated apoptosis. Cleavage of PAK2 regulates morphological changes in apoptotic cells and always correlates with apoptotic cell death (75). The variant PAK2_rs9325377*A occurs in strong linkage disequilibrium with the PAK2_rs6583176*A in the Iberian population (D' = 0.93), and both alleles were associated with increased susceptibility to gastric cancer (76). Thus, it is conceivable that they are associated with apoptosis failure, through yet unknown mechanisms. The association of the PAK2_rs9325377*A with higher PIGX expression further reinforces this possibility, since PIGX knockdown may inhibit breast cancer cell growth (77).

Taken together, decreased TNF gene expression (suggested by TNF_rs1800630*A association) increases PF susceptibility, whereas higher TRAF2 gene expression (suggested by TRAF2_rs10781522*G association) seems to protect against the disease, as well as PAK2_rs9325377*A/A. Moreover, TRAF2_rs10781522*G is associated with lower expression of C8G in the skin, expected to decrease the production of a complement component protecting keratinocytes against apoptosis and necroptosis (32). Higher TNF and TRAF2 levels are thus predicted to be protective against the disease, both preferentially leading to cell survival and inflammation. The fact that we did not find a clear association between any of these gene polymorphisms and apoptosis or necroptosis is in agreement with previous findings using electron microscopy (8, 22).

EIF2AK3, CD47, and SIRPA: Immunogenic Cell Death

We identified genetic associations of EIF2AK3_rs10167879*T with PF protection, and of SIRPA_rs6075340*A/A and CD47_rs12695175*G with susceptibility to PF. These genes encode products that participate in the immunogenic activation by CALR (calreticulin) (78), highlighting CALR as a molecule likely associated with damage-associated molecular patterns (DAMP) that may follow or precede the immunogenic pathway in PF. This cascade can be activated by a relatively limited set of stimuli, that may also include environmental factors associated with PF susceptibility, as UVB (79), thiol and other calcium-sequestering components (80), and components of fly saliva (71, 72) (Figure 2).

Figure 2

Proposed role of the immunogenic cell death in pemphigus foliaceus. 1: Keratinocytes exposed to exogenous PF risk factors, such as mosquito saliva (71, 72), ultraviolet irradiation (79), and thiol molecules (80), present several stress-derived peptides through HLA class I molecules. These peptides are recognized by cytotoxic T CD8+ cells and natural killer (NK) cells, which launch the immunogenic pathway. 2: This pathway induces damage-associated molecular patterns, including membrane exposure of CALR (81). CALR is transported from the endoplasmic reticulum to the outer cell membrane (81), a process mediated by EIF2AK3 (24). Calreticulin then interacts with SIRPA on the surface of dendritic cells, giving it a signal known as “eat me,” responsible for stimulating phagocytosis (81). Relaying the opposite signal is CD47, an outer surface membrane protein that also interacts with SIRPA and antagonizes the activity of CALR, inhibiting phagocytosis (81). We hypothesize that in PF the “eat me” signals prevail over the “don't eat me” signals, increasing the phagocytosis of keratinocyte debris by dendritic cells. 3: Dendritic cells then present keratinocyte peptides, as those derived from desmoglein 1, to T helper cells. 4: T lymphocytes stimulate auto-antibody production by B lymphocytes. 5: The immunogenic cell death process initiated by exogenous stimuli also activates an adaptive immune response, which includes recruitment, and activation of both cytotoxic T lymphocytes and natural killer cells.

Besides its association with PF protection, EIF2AK3_rs10167879*T is also associated with higher EIF2AK3 and lower EIF2AK3-DT (its antisense lncRNA) expression (32). This is the first documented association between this allele and any disease. The association is supported by the fact that deregulation of EIF2AK3 (also known as PERK) has been reported as one of the earliest pathogenic events in PV independently of IgG. Also, EIF2AK3 is activated in keratinocytes exposed to PV serum, as demonstrated by an increase in its phosphorylated levels and in phosphorylation of its target, eIF2α. Decreased EIF2AK3 expression, mediated by small interfering RNA, reduced the effects of PV serum on cell cycle and keratinocyte viability, two PV hallmarks (82). In agreement with this, Cipolla et al. (83) formerly proposed that PF IgG and/or non-IgG extracellular factors may lead to endoplasmic reticulum (ER) stress, resulting in C/EBP-homologous protein (CHOP) induction via EIF2AK3 (PERK), and activation of transcription factor 6 (ATF6). The association with EIF2AK3_rs10167879*T is puzzling, since it would imply that high EIF2AK3 levels are protecting against PF. Nevertheless, it may be explained in the light of its pleiotropy and a possible effect on differential transcription of its two poorly characterized isoforms (31). The SIRPA_rs6075340*A and CD47_rs12695175*G alleles are associated with lower expression of their respective genes (36), in addition to their synergistic association with PF susceptibility. The interaction between CD47 and SIRPA is capable of antagonizing the activity of surface-exposed CALR, responsible for emitting “eat me” signals for phagocytosis (81). The lower expression of those two physically interacting molecules is expected to lead to excessive internalization of cell debris and antigen presentation, increasing PF autoantibody production (Figure 2). This hypothesis agrees with our prior results on complement receptor polymorphisms, which would protect against PF development by modulating scavenging efficiency of acantholytic cell debris (7). A further argument in favor of the importance of the immunogenic cell death pathway in PF is the up-regulated expression of immunogenic deadly granzyme GZMA and GZMB genes in T lymphocytes of untreated patients with generalized lesions (84).

PRKN, HK1, and RAPGEF3: Pyroptosis, Parthanatos and Necrosis

PRKN_rs9355950*C was associated with PF protection (OR = 0.57, p = 0.0004). This gene encodes the mitophagy-regulating Parkin protein, which prevents cell death by causing ubiquitination of mitochondrial proteins presented by damaged mitochondria. Pyroptosis amplifies mitochondrial damage through caspase 1-driven cleavage and inactivation of PRKN (85). Since the associated effects of this variant on gene expression are unknown, the importance of pyroptosis on PF susceptibility remains elusive.

HK1_rs7072268*A was associated with increased susceptibility to PF (OR = 1.48, p = 0.0045) and also with higher gene expression levels in nervous tissue (32). Furthermore, this variant is associated with glycated hemoglobin levels (86). It is conceivable that it may be actually associated with enhanced poly(ADP-ribosyl)ation of hexokinase 1, occurring through the activation of PARP [poly(ADP-ribose) polymerase], with consequent glycolysis inhibition and induction of parthanatos-induced cell death (87).

Finally, the RAPGEF3_rs10747521*A/A genotype of the necrosis pathway, associated with protection against the disease (OR = 0.42, p = 0.004), is also associated with lower RAPGEF3 gene expression, as well as of its neighboring gene, SLC18A1. This protein has a role in inhibiting the p38MAPK pathway (88), activated in PF (83). Higher PCED1B-AS1 lncRNA expression is also associated with the same variant, predicted to enhance monocyte apoptosis and reduce autophagy (89). The reasons subjacent to this association await the results of future functional studies.

Conclusion

For the first time, SNPs located within genes involved in all known cell death cascades were systematically investigated in a single disease. The genetic association profile with TNF, TRAF2, CD36, and PAK2 variants favors cell survival and inflammation, instead of apoptosis/necroptosis, to explain resistance against the disease. On the other hand, susceptibility is conferred by variants of CD47 and SIRPA of the immunogenic cell death pathway, proposed to lead to excessive internalization of cell debris, and antigen presentation, which may increase PF autoantibody production. The importance of other pathways as pyroptosis, parthanatos, and necrosis, represented by one association each in our setting, cannot be disclosed and shall be further investigated. Functional validation of these associations, especially of genes encoding common isoforms, as EIF2AK3, will provide a better understanding of PF pathogenesis and contribute to the development of new drugs and to therapeutic improvement for the disease.

Data Availability Statement

The datasets analyzed for this study can be found in the LGMH repository: http://www.lgmh.ufpr.br/data/Supplementary_material_1_Bumiller-Bini_2019.xlsx; http://www.lgmh.ufpr.br/data/Supplementary_material_2_Bumiller-Bini_2019.xlsx.

Ethics Statement

This study was carried out in accordance with the guidelines of the Conselho Nacional de Ética em Pesquisa (CONEP) with written informed consent from all subjects. All subjects gave written informed consent in accordance with the Declaration of Helsinki. The protocol was approved by CONEP (No. 505.988).

Author Contributions

MP-E, AB, VB-B, GC, MS, DA, and MB contributed to conception of the work. AB and VB-B designed the study. AB, GC, and MB supervised the study. MP-E provided the samples. MP-E and AB provided resources. DA performed microarray hybridization. VB-B and MS did the statistical analysis and drafted the manuscript. AB, MP-E, GC, DA, and MB edited it. All authors revised the work critically for intellectual content and approved the final version of the work.

Conflict of Interest

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.