s2fieldsofstudy: Biology Medicine Chemistry

text:
NF-κB transcriptional activation by TNFα requires phospholipase C, extracellular signal-regulated kinase 2 and poly(ADP-ribose) polymerase-1

Background The nuclear enzyme poly(ADP-ribose) polymerase-1 (PARP-1) is required for pro-inflammatory effects of TNFα. Our previous studies demonstrated that PARP-1 mediates TNFα-induced NF-κB activation in glia. Here, we evaluated the mechanisms by which TNFα activates PARP-1 and PARP-1 mediates NF-κB activation. Methods Primary cultures of mouse cortical astrocytes and microglia were treated with TNFα and suitable signaling pathway modulators (pharmacological and molecular). Outcome measures included calcium imaging, PARP-1 activation status, NF-κB transcriptional activity, DNA damage assesment and cytokine relesease profiling. Results TNFα induces PARP-1 activation in the absence of detectable DNA strand breaks, as measured by the PANT assay. TNFα-induced transcriptional activation of NF-κB requires PARP-1 enzymatic activity. Enzymatic activation of PARP-1 by TNFα was blocked in Ca2+-free medium, by Ca2+ chelation with BAPTA-AM, and by D609, an inhibitor of phoshatidyl choline-specific phospholipase C (PC-PLC), but not by thapsigargin or by U73112, an inhibitor of phosphatidyl inisitol-specific PLC (PI -PLC). A TNFR1 blocking antibody reduced Ca2+ influx and PARP-1 activation. TNFα-induced PARP-1 activation was also blocked by siRNA downregulation of ERK2 and by PD98059, an inhibitor of the MEK / ERK protein kinase cascade. Moreover, TNFα-induced NF-κB (p65) transcriptional activation was absent in cells expressing PARP-1 that lacked ERK2 phosphorylation sites, while basal NF-κB transcriptional activation increased in cells expressing PARP-1 with a phosphomimetic substitution at an ERK2 phophorylation site. Conclusions These results suggest that TNFα induces PARP-1 activation through a signaling pathway involving TNFR1, Ca2+ influx, activation of PC-PLC, and activation of the MEK1 / ERK2 protein kinase cascade. TNFα-induced PARP-1 activation is not associated with DNA damage, but ERK2 mediated phosphorylation of PARP-1.

Background
Tumor necrosis factor alpha (TNFα) has pro-inflammatory effects in the nervous system and other tissues. These effects are due in part to its ability to induce nuclear factor kappa B (NF-κB) activation (as reviewed by [1,2]). The high-affinity TNFα receptors, TNFR1 and TNFR2, are both expressed on astrocytes and microglia of the central nervous system [3,4], and TNFα binding to either of these receptors initiates signaling cascades that activate the IKK kinase complex [1,2,5]. In the canonical NF-κB signaling pathway, the IKK kinase complex phophorylates the IκB subunit of cytosolic NF-κB complex, leading to the release and nuclear translocation of NF-κB dimers. These NF-κB dimers bind to κB sites in promoter regions of target genes, including many genes encoding pro-inflammatory proteins (as reviewed by [2,5]). However, nuclear translocation and DNA binding by NF-κB is not sufficient to induce gene transcription. This requires for NF-κB to form a functional transcription complex with other proteins. Evidence suggests that poly(ADP-ribose) polymerase-1 (PARP-1) is a key component of the NF-κB transcription complex [6][7][8]. PARP-1 forms branched poly(ADP) ribose (PAR) polymers on several proteins involved in DNA repair and transcription, including PARP-1 itself (PARP-1 automodification). In addition, many nuclear proteins express binding sites for these PAR polymers [9,10]. The protein-protein interactions resulting from PAR formation and PAR binding modulate the function of several transcription factors, including NF-κB [6,11,12], thus providing a link between PARP-1 and cellular inflammatory responses. PARP-1 enzymatic inhibitors can suppress inflammatory responses in microglia [11,[13][14][15] and other cell types [12,13,[16][17][18][19][20][21]. However, it has also been reported that enzymatically inactive PARP-1 can facilitate NF-κB transcriptional activation [22][23][24], suggesting that the antiinflammatory effects of PARP-1 inhibitors might occur by a different mechanism.
The signal transduction pathway leading from TNFα receptor activation to PARP-1 activation has not been established. PARP-1 activation is classically induced by DNA damage, in which PARP-1 binding to DNA strand breaks or kinks triggers enzymatic activation and poly(ADP-ribosyl)ation of histones and other proteins associated with DNA repair [22]. PARP-1 was recently shown to also be activated under several conditions in which DNA damage is not known to occur, including DNA transcription, DNA replication, and in neurons, neurotrophin responses and memory consolidation [25][26][27][28][29]. Studies suggest that ERK2 activation and release of endoplasmic reticulum Ca 2+ stores are involved in DNA damage independent of PARP-1 activation in neurons [30,31], but little is known about this signaling pathway in inflammatory cells or TNFα signaling.
Goals of the present study were to establish the mechanism by which TNFα leads to PARP-1 activation and to determine whether inhibition of PARP-1 enzymatic activity blocks TNFα-induced transcriptional activation of NF-κB. These studies used primary cell cultures of microglia and astrocytes in order to avoid artifacts inherent in the use of immortalized cell lines. Microglia are resident immune cells of the central nervous system, and astrocytes are glial cells that contribute to the inflammatory response in brain. Results of these studies show that TNFα leads to PARP-1 activation by a pathway independent of DNA damage, involving phosphatidyl choline-specific phospholipase C (PC-PLC), and ERK. The present study also demonstrated that inhibition of this pathway and inhibition of PARP-1 enzymatic activity both prevent NF-κB-mediated gene transcription.

Cell cultures
All animal studies were approved by the San Francisco Veterans Affairs Medical Center animal studies committee and University of Manitoba animal care and use committee, and follow the NIH guidelines for humane care of animals. Cultures were prepared from PARP-1 −/− , and wild-type (wt) mice as described previously [14]. The PARP-1 −/− mice were descendants of the 129S-Adprtl tmlZqw strain, originally developed by Wang et al. [32] and obtained from the Jackson Laboratory (Bar Harbor, ME). These mice were outbred for seven generations with wt CD-1 mice, and wt CD-1 mice were used as controls for the PARP-1 −/− mice. Astrocyte cultures were prepared from 1-day old mouse pups of both sexes in 24-well plates as described previously [33]. At confluency, microglia were harvested by gently shaking and re-plated at density of 5 × 10 5 cells/well on a 24-well plate [14]. The microglia cultures were subsequently maintained in glia-conditioned medium. The confluent astrocyte cultures were treated with 22 μM cytosine β-Darabinofuranoside for 2 days to inhibit the proliferation of remaining microglia. The astrocyte cultures were subsequently maintained in Eagle's minimal essential medium (MEM) supplemented with 3 % fetal bovine serum (FBS).

Drug incubations
Studies using astrocyte cultures were initiated by replacing the culture medium with a physiologically balanced salt solution (BSS) containing 3.1 mM KCl, 134 mM NaCl, 1.2 mM CaCl 2 , 1.2 mM MgSO 4 , 0.5 mM KH 2 PO 4 , 15.7 mM NaHCO 3 , and 2 mM glucose, pre-equilibrated to pH 7.2 in a 5 % CO 2 atmosphere. Drugs were prepared as concentrated stock solutions in BSS. Experiments were performed at 37°C in a 5 % CO 2 atmosphere and were terminated by complete medium exchange and replacement with fresh BSS. Studies with microglial cultures were performed identically but used MEM rather than BSS.

Immunostaining
Fixation and immunostaining of cell cultures was performed as previously described [14]. For detection of poly(ADP-ribose), incubations were performed with a 1:1000 dilution of rabbit anti-poly(ADP-ribose) (Trevigen, Gaithersburg, MD, Cat# 4336-BPC) at 4°C for 24 h, followed by incubation with Alexa Fluor 594-conjugated anti-rabbit IgG (Molecular Probes, Eugene, OR), 1:500 dilution, for 2 h at room temperature. For evaluation of NF-κB p65 subunit translocation, incubations were performed with a 1:30 dilution of rabbit polyclonal anti-p65 (Cell Signaling Technology, Inc., Danvers, MA, Cat# 3987) at 4°C for 24 h, followed by incubation with Alexa Fluor 594-conjugated anti-rabbit IgG (Molecular Probes, Eugene, OR), 1:500 dilution for 2 h. In some studies, the cells were subsequently incubated with 2 μg/ml propidium iodide for 5 min to obtain nuclear counterstaining. Confocal photomicrographs of astrocytes were obtained with a Zeiss LSM 510 META with a dual-PMT detector and a 32-channel META detector and Zeiss AIM imaging software. Microglia images were obtained with standard epifluorescence microscopy because these cells could not be cultured in a resting state on glass coverslips. Controls prepared in the absence of primary antibody showed no staining under the conditions described (not shown). DNA strand breaks were detected by DNA-polymerase I-mediated biotin-dATP nick translation (PANT) [34]. Cultures were permeabilized in 1 % Triton X-100 in phosphate buffered saline (PBS) for 30 min, and endogenous peroxidase was quenched with 2 % hydrogen peroxide. The cultures were then incubated for 90 min with the PANT reaction mixture: 5 mM MgCl2, 10 mM 2-mercaptoethanol, 20 μg/ml bovine serum albumin, 30 μM dGTP, 30 μM dCTP, 30 μM dTTP, 1 μM dATP, 29 μM biotinylated dATP, and 40 U/ml of DNA-polymerase I in PBS, pH 7.4. The biotinylated areas of DNA strand breaks were visualized by subsequent incubation with streptavidin-horseradish peroxidase and DAB peroxidase substrate. Controls incubated with PANT reaction mixture without DNA-polymerase I showed no staining.
siRNA downregulation of ERK1/2 Astrocytes were incubated with siRNA transfection complexes at day 10 in vitro, at which time they are 95 % confluent. siRNA incubations were performed with RNAiFect reagents (Qiagen Inc., CA) according to the manufacturer's instructions and as previously described [35]. The culture medium of each well was replaced with 400 μl of transfection complex mixture containing 1 μg of ERK1 siRNA or 3 μg of ERK2 siRNA and 6 μl of RNAiFect in optiMEM. The siRNA sense sequence for ERK1 was 5′ ACAAGCGCAUCACAGUAGAtt 3′ and for ERK2 was 5′ CAAAGUUCGAGUUGCUAUCtt 3′ (Ambion Inc, TX, USA). Controls were prepared with 3 μg of a mismatch sequence lacking significant homology to any known mouse gene sequences. The siRNA complexes were removed after 6 h and replaced with culture medium. Four days later when cultures were used for experiments, the percentage of siRNA-induced reduction in ERK1 and ERK2 expressions was determined by quantitating the ERK band size in ERK siRNA vs. mismatch siRNA transfected cultures in both ctrl and TNFα conditions.

Lentivirus NF-κB reporter gene transfection
The lentivirus NF-κB reporter gene construct (pLenti-kB-dEGFP) has a destabilized enhanced green fluorescent protein (dEGFP) construct under the control of five tandem repeats of κB-enhancer elements specific for p65and a plain TATA box [36]. Transfections were performed as previously described [15]. Astrocytes were infected with the Lenti-κB-dEGFP 4-5 days before experiments. EGFP expression was evaluated 1 and 24 h after exposure to TNFα by counting EGFP-expressing cells from five random fields within each culture well.

Calcium imaging
Intracellular calcium levels were measured in astrocytes loaded with Fura-2-AM (Molecular Probes; 3 μM, 60 min). Experiments were carried out in bicarbonatefree BSS at 37°C and buffered to pH 7.2 with 10 mM PIPES. Calcium-free BSS was prepared by omitting CaCl 2 and adding 2 mM EGTA. The change in dye fluorescence (Fura 340/380 ratio) was measured with aZeiss 200-M inverted microscope outfitted with an ORCA ER II microscope system using Openlab Improvision software. Excitation and emission wavelengths were controlled by Sutter filterwheels and shutters. Baseline values were recorded for 15-20 min prior to drug additions. At the end of each experiment, the calcium ionophore A23187-Br (10 μM) was applied to calibrate the intracellular signal. Data are presented as a change in Fura2 emission (Fura 340/380 ratio), and as a mean change in Fura 340/380 ratio during 15-, 30-, and 60-min intervals from three to four independent experiments. In each independent experiment data were acquired from six to ten cells at 90-s intervals, and at 10-s intervals in neutralizing antibody experiments. The high-affinity dye, Fura-2, was used to maximize sensitivity at low intracellular calcium concentrations, and therefore, the signal is likely saturated and underestimates the differences between Ca 2+ -free and Ca 2+ -containing medium.

Western blotting
Western blots were prepared and quantified as described [14]. Membranes were incubated with 1:1000 dilutions of rabbit anti-ERK1/2 or anti phospho-ERK1/2 (Cell Signaling Technology, Beverly, MA, Cat# 9102 and 4376), or with a 1:500 dilution of mouse monoclonal anti-poly(ADP-ribose) (Trevigen, Gaithersburg, MD, Cat# 4335-MC). After washing, the membranes were incubated for 2 h with peroxidase-conjugated anti-rabbit or anti-mouse IgG (Vector Laboratories, Burlingame, CA) diluted 1:7500. The protein bands were visualized using ECL TM WB Detection kit (Amersham-Pharmacia Biotech) and X-OMAT AR film (Kodak) or ChemiDoc MP imaging system (BioRad). To quantify protein loading, the membranes were re-probed with mouse monoclonal anti β-actin at a 1:10,000 dilutions, followed by peroxidase-conjugated anti-mouse IgG (1:10,000 dilutions). Controls performed in the absence of primary or secondary antibodies showed no signal (data not shown). Band densities were quantified with the NIH Image J program.

Site-directed PARP-1 mutations
The mutant PARP-1 constructs were prepared as described [35]. In short, human PARP-1 cDNA (hPARP-1, BC03754 from NIH Mammalian Gene Collection) was corrected to match the human genomic sequence at amino acid residue 762 (from alanine to valine) using Site-Directed Mutagenesis kit (Stratagene). Single (S372E) or double mutants (S372A and T373A) of these constructs were then generated and sequence-verified. PARP-1 −/− astrocytes were transfected 2 days before experiments with the PARP-1 constructs by replacing the culture medium with 300 μl of transfection complex mixture containing 0.5 μg of DNA and 1 μl of Lipofectamine 2000 (Invitrogen) in OptiMEM. The medium was replaced with culture medium 8 h after the initiation of transfection and cells were used for experiments 40 h later. The transfection efficacy of this method has been established previously by western blots detecting PARP-1 expression with a polyclonal antibody (ALX-210-302, Alexis/Enzo Life Sciences, Plymounth Meeting, PA) [35].

Cytokine assay
Cytokines were analyzed in 50-μl aliquots of cell culture medium using a Milliplex mouse multiplex immunoassay bead system according to the manufacturer's instructions (Millipore) [37]. Samples were assayed in duplicate, and the fluorescent signal corresponding to each cytokine was measured with a BioPlex 200 system (BioRad, Hercules, CA) in parallel with known standards. Cytokine concentrations measured from culture medium were normalized to the protein content of each well as determined by the bicinchoninic assay [38].

Statistics
For cell culture studies, each "n" denotes an independent experiment, where each experiment is from a separate cell culture preparation and comprised of 3-4 parallel Immunostaining for poly(ADP-ribose) (labeled as PAR) shows accumulation in both astrocytes and microglia during incubation with TNFα (15 ng/ml). Phase contrast images shows simultaneous morphological transformation of microglia. Astrocytes do not exhibit morphological changes (not shown). Scale bar, 10 μm. Representative of n = 4 treatments per condition. Results are presented as a means ± standard error. Statistical significance was evaluated by one-way ANOVA followed by the Student-Newman-Keuls multiple comparison test. p values below 0.05 were considered significant.

Results
TNFα-induced PARP-1 activation in the absence of detectable DNA damage TNFα at physiological concentration (15 ng/ml) induced a rapid accumulation of PAR polymers, a product of PARP-1 enzymatic activity, in both astrocyte and microglial cultures (Fig. 1). PAR formation in microglia was time-linked to morphological activation, as characterized by process retraction and soma enlargement (Fig. 1). Astrocytes, unlike microglia, do not undergo morphological changes in response to TNFα (not shown). Prior studies have shown that TNFα does not induce PAR accumulation in PARP-1 −/− cells [14], indicating that PARP-1 is the major source of PAR formation [15,39,40].
PARP-1 activation is classically induced by DNA damage, particularly DNA strand breaks [22,41]. To determine if DNA damage was involved in TNFα-induced PARP-1 activation, cultures were evaluated with the DNA-polymerase I-mediated biotin-dATP nick translation (PANT) technique for single-strand DNA breaks. Cultures fixed after 1 h of TNFα incubation, a time point with robust PAR accumulation, showed no increase in PANT staining. By contrast, an extensive PANT signaling was observed in cells treated with 10 μM N-methyl-N'-nitro-N-nitrosoguanidine (MNNG), a DNA alkylating agent widely used to induce PARP-1 activation [42]. The degree of PARP-1 activation induced by this concentration of MNNG was comparable to that observed in cultures treated with TNFα (Fig. 2c).
Studies with a TNF receptor 1 (TNFR1)-neutralizing antibody were performed to assess the role of this receptor subtype in PARP-1 activation. The neutralizing antibody delayed and reduced TNFα-induced intracellular Ca 2+ elevations (Fig. 5a, c, d). The TNFR1 neutralizing antibody also reduced PAR formation (Fig. 5e), suggesting that TNFR1 has a critical role in PARP-1-related inflammatory responses.

ERK2 phosphorylation of PARP-1 induces NF-κB activation
The role of ERK2-mediated phosphorylation in PARP-1 activation was further evaluated using PARP-1 constructs containing mutations at the ERK2 phosphorylation sites, S372 and T373 [35]. The S372A and T373A mutations prevent ERK2-mediated phosphorylation, and S372E provides a constitutive phosphorylation mimic. PARP-1 −/− astrocyte cultures were transfected with either wt or mutant PARP-1 constructs, and NF-κB transcriptional activation was assessed using a lentivirus-delivery of dEGFP reporter gene that is activated by the p65 subunit of NF-κB [36]. In cells transfected with wt PARP-1, TNFα induced a robust activation of the reporter gene within 1 h and a further increase at 24 h. This TNFα-induced transcriptional activation was absent in non-transfected cells. It was also absent in transfected cells treated with the PARP inhibitor DPQ, thus confirming a requirement for PARP-1 enzymatic activity in this process. Cells transfected with the mutant PARP-1 lacking ERK2 phosphorylation sites (S372A and T373A) also showed no TNFα-induced NF-κB activity. By contrast, cells transfected with the S372E PARP-1 construct, which mimics constitutive phosphorylation by ERK2, showed increased GFP expression in the absence of TNFα stimulation (Fig. 7). Of note, neither PARP-1 deficiency nor PARP-1 inhibition blocked the nuclear translocation of p65 induced by TNFα (Fig. 8). We also analyzed TNFα-induced cytokine release and saw decrease in NF-κB-dependent pro-inflammatory cytokines . a TNFα induced increase in intracellular free Ca 2+ , as indicated by the change in Fura2 emission. This increase was blocked in Ca 2+ -free medium (0 Ca 2+ + 0.5 mM EGTA, green), and the intracellular calcium response profile was changed by neutralizing TNFR1 antibody (TNFR1 ab, red). IgG antibody (IgG ab, blue) was used as a control. Arrow denotes addition of TNFα (15 ng/ml) to the astrocyte cultures. The Ca 2+ ionophore, A23187-Br (10 μM) was added at the end of each experiment (arrow head) to provide a calibration point. b TNFα-induced intracellular free Ca 2+ (Fura2 ratio) increase was prevented in Ca 2+ -free medium but was not affected in cells pre-treated with thapsigargin (2 μM). * p < 0.05 compared to control at each time point; # p < 0.05 compared to TNFα at each designated time point, n = 5-6. c, d TNFR1 neutralizing antibody (TNFR1 ab) delayed calcium influx rise onset and reduced the influx plateau time, but the calcium influx maximum of the level of intracellular free Ca 2+ (Fura2 ratio) was not affected between TNFR1 ab and IgG ab-pre-treated TNFα-stimulated cells. * p < 0.05, n = 5. e PARP-1 activation profile detected by PAR formation was also changed when TNFR1 was blocked. TNFR1 neutralizing antibody reduced the length and level of PAR formation, mimicking effects of PARP inhibitor. * p < 0.05, n = 3 Fig. 6 TNFα-induced PARP activation is regulated by ERK2. a TNFα (15 ng/ml) induces ERK1/2 phosphorylation in astrocyte cultures. Representative of n = 3. b PAR immunoblots shows that TNFα-induced PAR formation is reduced after siRNA downregulation of ERK2, but not ERK1. PAR is seen as a discrete band at 116 kDa attributable to auto-poly(ADP)-ribosylation of PARP-1 and as a smear over other molecular weights resulting from PAR formation on other proteins. Actin band shows protein loading. ERK1 and ERK2 expression are shown below the PAR western blot. The siRNA treatments reduced their expression by 72 ± 8 % and 72 ± 5 %, respectively, with protein loading shown by the actin band. c Quantification of PAR western blots, conditions as in b; * p < 0.05, n = 3. d Immunostaining for PAR also shows that TNFα-induced PAR formation is suppressed by siRNA downregulation of ERK2 expression. Images are representative of four independent experiments. All PAR formation studies are conducted 60 min after TNFα stimuli  (Table 1), providing more functional outcome measure for PARP-1 mediated NF-κB regulation. Interestingly, PARP-1 inhibition did not reduce anti-inflammatory cytokines (IL-4, IL-10, IL-13, TGFβ) ( Table 1).

PLC does not have a role in PARP-1 activation induced by DNA damage
Taken together, these findings suggest a signaling pathway in which TNFα interaction with TNFR1 triggers prolonged Ca 2+ influx, PC-PLC activation, ERK2 phosphorylation, and activation of PARP-1. We therefore evaluated whether this pathway is involved in PARP-1 activation induced by DNA damage. Astrocytes were exposed to 10 μM MNNG to induce DNA damage sufficient to produce PARP-1 activation comparable to that induced by TNFα (Fig. 2c). DNA damage-induced PARP-1 activation, like TNFαinduced activation, was blocked by the MEK/ERK inhibitor PD9089; however, it was not blocked by the PLC inhibitor D609 or by Ca 2+ -free medium and was only partly blocked by BAPTA-AM (Fig. 9).

Discussion
The mechanism by which PARP-1 is activated by TNFα and mediates TNFα pro-inflammatory effects has not been previously established. Results of these studies suggest that TNFα induces PARP-1 activation by a pathway involving TNFR1, calcium entry, activation of PC-PLC, and activation of the MEK1/ERK2 protein kinase cascade. This signaling pathway differs in key respects from the pathway involved in classical PARP-1 activation by DNA damage (Fig. 10). These results also show that PARP-1 enzymatic activity is required for TNFα-induced NF-κB transcriptional activation and related cytokine production in glia.
NF-κB activation involves its translocation to the nucleus, transcription complex formation and binding to DNA [2]. However, these events alone do not necessarily  Data is presented as a pg of cytokine detected from the cell culture medium, normalized to mg of protein in corresponding cell culture well (mean ± stdev). PJ34 is a PARP-1 inhibitor. * p ≤ 0.05 compared to ctrl (wt), ** p ≤ 0.05 compared to TNFα (wt), n = 3-4 result in NF-κB-driven gene transcription; transcription of the target gene can be prevented by parallel binding of repressors or other contradicting processes [12,52,53]. Results presented here with an NF-κB reporter gene assay provide an example of this; PARP-1 inhibitors allowed NF-κB translocation to the nuclei, but blocked NF-κBmediated gene transcription of pro-inflammatory cytokine, while some anti-inflammatory cytokines (IL-10, IL-13), which have a κB site in their promoter [54,55] were not reduced by PARP-1 ablation. The reporter gene approach also clarifies the importance of PARP-1 enzymatic activity in NF-κB activation. Several reports indicate that PARP-1 inhibitors can block glial inflammatory responses and the associated upregulation of proteins with NF-κB promotor sites [11, 13-15, 18, 56], but other reports suggest that PARP-1 enzymatic activity is not required for NF-κB transcriptional activation [24]. The present results demonstrate that PARP-1 enzymatic activity is required for transcriptional activity induced by the p65 subunit of NF-κB in astrocytes. However, we cannot exclude a non-enzymatic interaction between PARP-1 and other NF-κB DNAbinding subunits or in non-canonical NF-κB activation [57]. Both genetic depletion and pharmacologic inhibition of PARP-1 were found to inhibit NF-κB-dependent proinflammatory cytokines, but interestingly, anti-inflammatory cytokines were not inhibited. This suggests that immune modulation ability of PARP-1 [37] is a result of complex transcriptional regulation, possibly involving other transcription factors [58]. Although PARP-1 activation is classically associated with DNA damage, in which PARP-1 binding to DNA strand breaks or kinks triggers its enzymatic activation, PARP-1 activation also occurs in settings such as gene transcription and neurotrophin signaling, in which DNA damage is not known to occur [25][26][27][28][29][30][31]. PARP-1 activation induced by TNFα would similarly be expected to occur independent of DNA damage, but it is also possible that TNFα could upregulate pro-oxidant inflammatory responses and thereby lead to DNA damage [59]. This possibility was excluded here using the PANT assay, which showed no Our results indicate some commonalities in the signal transduction pathways mediating DNA damage-induced PARP-1 activation and TNFα-induced PARP-1 activation. PARP-1 activation induced by DNA damage requires PARP-1 to be phosphorylated at an ERK2 phosphorylation site [35]. Here, PARP-1 activation induced by TNFα was found to be blocked by siRNA downregulation of ERK2 expression and by PD98059, which inhibits ERK2 activation. In addition, cells expressing mutant PARP-1 that lack the ERK2 phosphorylation at sites S372 and T373 showed no PARP-1 activation [35] and no NF-κB transcriptional activation, whereas cells expressing mutant PARP-1 containing the phosphomimetic S372E mutation showed an increase in both basal PARP-1 activity and basal NF-κB transcriptional activity. These findings suggest that ERK2 phosphorylation of PARP-1 at S372 and/or T373 is required for PARP-1 activation triggered by either DNA damage or by TNFα. We cannot exclude the possibility that other mitogen-activated protein kinases (MAPKs) have a role in TNFα-induced PARP-1 activation. TNFα certainly stimulates phosphorylation of p38 MAPK and cJun N-terminal kinase (JNK) [60], which could at least indirectly facilitate PARP-1 activation or parallel processes. Indeed, activation of JNK is a required downstream event in PARP-1-mediated (MNNG-induced) necrotic cell death [61,62], and p38 interaction with PARP-1 is a crucial step in hydrogen peroxide-induced cell death pathway [63]. Notably, phosphorylation of both of these MAPKs were identified downstream of PARP-1 activation. Our previous studies assessed the role of p38, JNK, and ERK1/2 in MNNGinduced PARP-1 activation and cell death, and only ERK1/2 inhibition reduced both [35].
Our results also indicate differences in the DNA damage-induced and TNFα-induced PARP-1 activation pathways. PARP-1 activation by MNNG, a DNA alkylating agent, was not inhibited by the PC-PLC inhibitor D609 or by Ca 2+ -free medium and only partially inhibited by the intracellular Ca 2+ chelator BAPTA-AM. By contrast, TNFα-induced PARP-1 activation was completely blocked by D609, Ca 2+ -free medium, and BAPTA-AM. The Ca 2+ i increase induced by TNFα was abrogated by Ca 2+ -free medium but not by thapsigargin pre-treatment, suggesting that the increased Ca 2+ i originates from the extracellular space rather than the ER. Fig. 10 The schematic summary of TNFα-induced PARP-1 activation pathway in glial cells. PARP-1 activation induced by TNFα does not involve DNA damage, which is a classical inducer of PARP-1 activation. A MNNG-induced DNA damage activation pathway differs from a TNFα-induced one, but they both share the requirement of ERK2-mediated PARP-1 phosphorylation. TNFα binding to its receptor TNFR1 induces calcium influx, and activates phosphotidyl choline-specific phospholipase (PC-PLC) leading to phosphotidyl choline cleavage to diacyl glycerol (DAG). This process leads to activation of mitogen-activated protein kinase pathway (MEK1/2-ERK1/2) and ERK2-mediated phosphorylation of PARP-1, which allow PARP-1 activation. Enzymatically active PARP-1 co-activates NF-κB transcriptional activation promoting pro-inflammatory secretion profile. TNFα-induced PARP-1 activation was prevented by calcium chelators (EGTA, BAPTA-AM), TNFR1-neutralizing antibody (ab), inhibitors of PC-PLC-DAG (D609), MEK1/2-ERK1/2 (PD98059), PARP (DPQ) and ERK2 siRNA, while MNNG-induced PARP activation was affected only by inhibitors of ERK2 and PARP. TNFα-induced NF-κB transcription activation was prevented by inhibitors of PC-PLC-DAG, ERK1/2 and PARP Also arguing against a role for ER calcium release was the finding that TNFα-induced PARP-1 activation was not blocked by U73112, an inhibitor of PI-PLC-induced IP 3 formation, but was inhibited by D609, which inhibits PC-PLC-induced DAG formation. Similarly, neutralizing antibody against TNFR1, extracellular receptor for TNFα, significantly diminished the TNFα-induced increase in Ca 2+ i and PAR formation. Moreover, the PLC-activating compound m-3M3FBS induced PARP-1 activation and NF-κB transcription activity, which were prevented by inhibitor targeting PLC-DAG signaling in the absence of other stimuli, suggesting that this pathway is sufficient to induce PARP-1 activity independent of DNA damage. These findings are congruent with prior work; TNFα binding to its TNFR1 receptor can stimulate PC-PLC activation and production of DAG (as reviewed by [50]. DAG precedes PKC activation, which in turn activates the MEK/ERK protein kinase cascade in astocytes [64]. In macrophage cultures, LPS-induced activation was reported to involve PC-PLC, DAG, PKCζ, and ERK1/2 [65]. While PKCs are likely mediators between DAG and MEK/ERK activation also in the signaling pathway described in our present report, further studies are required to identify whether they have a role in PARP-1 activation and which subunits are involved. The potential route for calcium influx could be through store-operated channels (SOCs), which are known to be regulated by TNFα-induced K + current upon TNFα stimulation in human microglia [66]. PLC-DAG activation has also been reported to trigger opening of SOCs [67,68].
Our data demonstrate the critical role of TNFR1 in PARP-1-related inflammatory responses, which is fitting as TNFR1 is a dominant pro-inflammatory receptor expressed by microglia and also associated with NF-κB-dependent anti-apoptotic transcriptional activity [69]. However, we cannot completely exclude the role of TNFR2, as the data relies on specificity of TNFR1 antibody, and there is no commercially available neutralizing antibody for TNFR2. It also remains to be established how the described signaling components (i.e., PC-PLC, calcium influx, and ERK2) upon TNFR1 stimulation interact with the classical TNFR adaptor proteins sucs as TRADD, TRAF-2, and RIP1.
A key role for ERK2 in DNA damage-independent PARP-1 activation has been previously attributed to a non-enzymatic interaction between ERK2 and PARP-1 [30,31]. Prior studies have also suggested a role for PLCmediated liberation of IP 3 and ER calcium from within the nucleus itself as a factor stimulating PARP-1 activity [31]. The present findings do not address this possibility directly, but the failure of thapsigargin or D609 to block PARP-1 activation, along with the suppression of PARP-1 activation in Ca 2+ -free medium, argues against ER Ca 2+ release as a major factor in the setting of TNFα stimulation. Whether these experimental differences are attributable to cell-type differences, differing cell stimuli, or other factors remains to be established. PARP-1 is by far the most abundant and active PARP species in the cell nucleus, but PARP-2 and PARP-3 are also present [40,70,71]. Recent findings have identified distinct roles for PARP-2 and PARP-3 in the astrocyte inflammatory responses induced by S. Aureus [72]. A major contribution of PARP-2 or PARP-3 to the PAR formation and NF-κB activation observed in the present studies is unlikely because no PAR formation or NF-κB activation was observed in PARP-1 deficient cells. However, the differing PARP species interact [73,74]. It is therefore possible that the different PARP species and interactions between these species may have signalspecific effects on transcription factor activities.

Conclusions
TNFα stimulation of glial cells leads to PARP-1 activation via a mechanism that is independent of DNA strand breaks. TNFα stimulates PARP-1 activation by a pathway involving TNFR1, calcium entry, activation of PC-PLC, and activation of the MEK1/ERK2 protein kinase cascade. PARP-1 enzymatic activity is required for TNFα-induced NF-κB transcriptional activation and pro-inflammatory cytokine release.


extfieldsofstudy: Biology Medicine

text:
Assessing the Pathogenicity of Insertion and Deletion Variants with the Variant Effect Scoring Tool (VEST‐Indel)

ABSTRACT Insertion/deletion variants (indels) alter protein sequence and length, yet are highly prevalent in healthy populations, presenting a challenge to bioinformatics classifiers. Commonly used features—DNA and protein sequence conservation, indel length, and occurrence in repeat regions—are useful for inference of protein damage. However, these features can cause false positives when predicting the impact of indels on disease. Existing methods for indel classification suffer from low specificities, severely limiting clinical utility. Here, we further develop our variant effect scoring tool (VEST) to include the classification of in‐frame and frameshift indels (VEST‐indel) as pathogenic or benign. We apply 24 features, including a new “PubMed” feature, to estimate a gene's importance in human disease. When compared with four existing indel classifiers, our method achieves a drastically reduced false‐positive rate, improving specificity by as much as 90%. This approach of estimating gene importance might be generally applicable to missense and other bioinformatics pathogenicity predictors, which often fail to achieve high specificity. Finally, we tested all possible meta‐predictors that can be obtained from combining the four different indel classifiers using Boolean conjunctions and disjunctions, and derived a meta‐predictor with improved performance over any individual method.

Introduction
The average human exome contains over four hundred naturally occurring microinsertion/deletion variants (referred to as indels throughout) [Consortium, 2010]. In-frame indels account for Additional Supporting Information may be found in the online version of this article. † These authors contributed equally to this work. * Correspondence to: Rachel Karchin ß50% of these variants, and result from the insertion/deletion of an integer number of codons, and ultimately amino acids. Frameshifts account for the other ß50% of indels, and result from contiguous nucleotide insertions/deletions of a length not divisible by three. This fractional change in the number of codons shifts the translational reading frame, resulting in an entirely new downstream sequence thereby shifting the position at which the first stop codon is encountered. Thus, frameshift indels translate to protein that is very distinct from the native protein, particularly if the indel occurs early in the transcript sequence. Both in-frame and frameshift indels alter protein sequence and length.
Because they drastically alter protein primary structure, but are also highly prevalent in healthy populations, indels present a unique classification challenge. Clearly, the principles governing pathogenicity are not identical to those governing changes in protein function and stability; otherwise, most indels would be pathogenic. The challenge then arises because protein sequence, structure, function, and stability are typically considered when assessing a variant of unknown impact on disease liability [Steward et al., 2003]. These protein-based criteria could lead to a high false positive rate, and therefore low specificity, because the fraction of indels that appreciably impact health might be overestimated.
The increased utilization of high-throughput genomic sequencing technologies and hopes for their clinical application, coupled with the high prevalence of indels, has led to a demand for bioinformatic predictors of indel pathogenicity. Most computational methods for assessing genetic variation initially focused on missense variants; more recently, several groups have extended these methods to handle indels [Choi et al., 2012;Hu and Ng, 2012;Hu and Ng, 2013;Zhao et al., 2013;Folkman et al., 2015]. Most of these methods utilize supervised machine learning classifiers and are trained on two classes of indel: pathogenic from disease mutation databases and benign from either population variation databases or tolerated interspecies variations derived from genomic alignments. DDIG-in is based on a support vector machine, and the authors of this method reported a sensitivity of 0.86 and specificity of 0.72 for frameshift indels [Zhao et al., 2013], and a sensitivity of 0.89 for in-frame indels [Folkman et al., 2015]; the authors did not report prediction specificity for in-frame indels. PROVEAN uses an unsupervised approach that compares the reference protein sequence with a sequence that incorporates a variant of interest [Choi et al., 2012]. The authors of PROVEAN reported high sensitivities of 0.93 and 0.96 for inframe insertions and deletions, respectively, and a specificity of 0.80 for in-frame insertions and 0.68 for in-frame deletions; PROVEAN does not assess frameshift indels. SIFT-indel, based on a J48 Decision Tree [Hall et al., 2009], achieved good balanced accuracies for in-frame (sensitivity = 0.81; specificity = 0.82) [Hu and Ng, 2013] and frameshift indels (sensitivity = 0.90; specificity = 0.78) [Hu and Ng, 2012]. However, the neutral dataset used in those studies comprised indels derived from cross-species comparisons. As the authors state, SIFT-indel was trained to predict impact on gene function, irrespective of impact on disease. Indeed, when the method was applied to variants from human variation databases, the majority of the indels were predicted to be deleterious; thus, specificities would be below 50% for predicting indel pathogenicity. The CADD classifier utilized a unique approach, in which a support vector machine was trained to discriminate fixed (or nearly fixed) derived alleles in humans from a set of simulated variants [Kircher et al., 2014]. The CADD classifier was developed to predict deleterious variants rather than variant pathogenicity or impact on protein function, but with the stated assumption that these quantities are all related. The authors of CADD reported classifier performance on missense variants and indels together, but not on indels separately [Kircher et al., 2014].
We present a novel method for assessing the impact of indels with an emphasis on discriminating between benign and disease-causing indels. This new approach expands the functionality of our Variant Effect Scoring Tool (VEST), which was initially designed to assess the impact of missense variants; we refer to the new functionality as VEST-indel, throughout. In addition to more conventional protein sequence and functional considerations, VEST-indel includes a feature based on PubMed search results for the gene of interest, which is a measure of known relevance to human health. Utilizing one dataset for cross-validation and a second, non-overlapping dataset for testing, we report high sensitivity and specificity for both inframe and frameshift indels. We provide direct comparison between VEST-indel and existing indel prediction methods by constructing an additional benchmark dataset of neutral and disease-causing indels on which none of the methods had been trained, and show that VEST-indel achieves improved performance. Next, we exhaustively test Boolean combinations of all tested classifiers and identify a meta-predictor that further improves performance relative to any of the methods individually. Finally, we used Tajima's D statistic [Tajima, 1989] to detect signatures of positive, balanced, or relaxed selection for variants that resulted in false positive predictions. The presence of these selective pressures could suggest a context dependence in which indels defined as neutral could also be considered pathogenic, in turn justifying false positive predictions.

Data Collection
A curated set of in-frame and frameshift indels (microdeletions and microinsertions) of ࣘ20 base-pairs in length, annotated as being pathogenic from publications in the biomedical literature, was downloaded from Human Gene Mutation Database [Stenson et al., 2014. Only high-confidence annotations with the "DM" designation were included. A second curated set of in-frame and frameshift indels was downloaded from the NCBI ClinVar database on August 7, 2014. Only entries annotated as "likely pathogenic" (Clinical Significance 4) or "pathogenic" (Clinical Significance 5) and not annotated as a somatic mutation were included. Any entry from ClinVar that was also present in HGMD was removed from the ClinVar set. Annotated in-frame and frameshift variants were downloaded from the Exome Variant Server using (ESP6500SI-V2-SSA137) [Fu et al., 2013] and from the 1000 Genomes Project Phase 3 (ftp://ftp-trace.ncbi.nih.gov/1000genomes/ftp/) [Clarke et al., 2012]. To increase the likelihood that variants from the Exome Variant Server and 1000 Genomes Project were benign common polymorphisms, and to retain sufficient variants for our training set, we only used variants with a minor allele frequency (MAF) ࣙ0.01 and occurring in either African individuals or those of African ancestry. In ESP600, these were identified as "African-American" and in 1000G as the AFR superpopulation comprising YRI (Yoruba in Ibadan, Nigeria), LWK (luhya in Webuye, Kenya), GWD (Gambian in Western Divisions in the Gambia), MSL (Mende in Sierra Leone), ESN (Esan in Nigeria), ASW (Americans of African ancestry in SW USA), and ACG (African Caribbeans in Barbados). The other populations represented in ESP6500 and 1000G are believed to have experienced severe bottlenecks in recent history, and hence individuals from these populations may harbor potentially pathogenic variants at higher MAF than individuals of African ancestry [Lohmueller et al., 2008;MacArthur and Tyler-Smith, 2010;Ng et al., 2008;Hu and Ng, 2012]. A curated set of putatively benign in-frame and frameshift indels, derived from pairwise genome alignments of human and cow, dog, horse, chimpanzee, rhesus macaque, and rat, was generously provided to us by Pauline Ng and Jing Hu. This set had been previously used to train their SIFT-indel classifier [Hu and Ng, 2012;Hu and Ng, 2013]. Additional background information about these data sets, including probability densities for indel length and MAF, are shown in Supp. Figure S1.

Feature Selection
The Random Forest Feature Importance Z-score [Breiman, 2001] was used to rank a set of 49 candidate features from [Wong et al., 2011] and five additional features (Supp. Table S1), using PARF software (http://code.google.com/p/parf), with 100 trees and default parameters. To avoid overfitting, an independent feature-selection set was used (500 pathogenic and 500 benign examples for each of the in-frame and the frameshift classifiers). We used a greedy algorithm to identify a good, minimum set of features. Briefly, beginning with the top-ranked feature, a Random Forest was trained using only that feature and 10-fold cross-validation was used to estimate the classifier's area under the receiving operator characteristic (ROC) curve (AUC). We successively added the next top-ranked feature until all candidate features were included. For the in-frame classifier, the maximum AUC was achieved with 23 features and for the outof-frame classifier, the maximum AUC was achieved with 16 features (Supp . Table S2). These features were used for the remainder of the work described here. The selected features include measures of gene importance, the damaging effect of the variant on protein activity, evolutionary conservation and protein local environment (Supp.

Classifier Training Protocol
Random Forest classifiers were trained to classify in-frame and frameshift variants (using PARF software with 100 trees and default parameters). For in-frame classifier training, 2,475 pathogenic and 1,877 benign examples were available, whereas 24,478 pathogenic and 1,350 benign examples were available for frameshift classifier training (Table 1). To handle class imbalance, the in-frame classifier was trained on a randomly selected set of 1,877 pathogenic examples and all 1,877 benign examples. Ten frameshift classifiers were trained on a randomly selected set of 1,350 pathogenic examples and all 1,350 benign examples (repeated 10 times, sampling without replacement). All ten classifiers were used to score frameshift variants, by computing 10 scores for each variant and averaging them.

Modeling the Analytical Null Distribution
We developed an analytical null score distribution based on Random Forest classifier scores of putative benign variants (1000 Genomes Project Phase 3, MAF ࣙ0.01, African ancestry). The scored variants (537 for in-frame insertion/deletions, 346 for frameshift insertion/deletions) did not overlap the examples used for Random Forest feature selection, training or the independent test set. An empirical cumulative distribution of scores was calculated and modeled as a Generalized Pareto Distribution (GPD) [Pickands III, 1975;Knijnenburg et al., 2009].
where k and a are the GPD shape and scale parameters, respectively (in-frame k = 0, a = -8.48, frameshift k = 0, a = -9.04) and z is the score (0 ࣘ z ࣘ 1).

Analytical Null Distribution for Missense Variants
We developed an analytical null distribution by modeling a GPD based on an empirical null distribution of VEST missense Random Forest classifier scores (Eq (1), k = 0, a = -5.26). The missense Random Forest  was trained on 47,724 pathogenic missense mutations from HGMD v2012.2 and 45,818 putatively benign missense variants from ESP6500 accessed 07/2012. The empirical null scores were generated from 28,509 variants from the 1000 Genomes Project with AF ࣙ 1%, none of which overlapped with training data.

Combined Prioritization of Indel and Missense Variants
For each In-frame, frameshift, or missense variant, a VEST score was computed using homology-restricted 10-fold cross-validation with the appropriate Random Forest (Performance Assessment). Then, a P value was calculated using the analytical null (Eq. (1)) for its respective type. To assess whether VEST could correctly prioritize a pathogenic variant over a benign variant, irrespective of whether the variant was in-frame, frameshift or missense, we ranked the combined set of variants according to p-value, and computed area under the ROC curve [Fawcett, 2004].
To compare VEST results with CADD, which also provides combined prioritization, we scored the same variants using the CADD Webserver (Materials and Methods). Variants were ranked according to their scaled C-scores and area under the ROC curve was computed.

Homology-restricted cross-validation
The in-frame indel Random Forest and each of the 10 frameshift Random Forests were assessed for sensitivity, specificity and balanced accuracy, using a rigorous ten-fold cross-validation protocol. The same protocol was applied to the missense Random Forest to assess combined prioritization of all three mutation types. To avoid overestimating our performance, we ensured that any examples from genes whose protein products had ࣙ35% sequence identity were included in the same fold. BlastP with default parameters [Altschul et al., 1997] was used for pairwise alignment of protein sequences and sequence identity calculations [Altschul et al., 1997]. Evidence suggests that homology-restricted cross-validation is important to avoid overly optimistic estimates of pathogenicity classifier performance [Capriotti and Altman, 2011].

VEST-indel independent test set
An independent set of examples, having no overlap with feature selection, training, or empirical null sets, was constructed. We removed any test set examples whose protein products had ࣙ35% sequence identity with any training examples [Altschul et al., 1997]. Pathogenic in-frame and frameshift mutations were taken from ClinVar [Landrum et al., 2013] and benign in-frame and frameshift variants were taken from the interspecies set (Data Collection). All data were "cleaned" so as to ensure that there was no overlap between these examples and the other three data sets. The VEST-indel independent test set is included in Additional File 1 (http://karchinlab.org/vest_indel_additional_files/Additional_File_ 1.xlsx).

Comparison of Insertion/Deletion Variant Pathogenicity Predictors
Four previously published methods were selected for comparison with VEST-indel. Three of the methods (SIFT-indel, DDIG-in, CADD) [Hu and Ng, 2012;Hu and Ng, 2013;Zhao et al., 2013;Folkman et al., 2015] handle both in-frame and frameshift insertion/deletions and the fourth method PROVEAN [Choi et al., 2012]   handles only in-frame variants. According to their publications, each of these methods used a customized dataset for classifier training and validation. Data sources included the UniProtKB, HGMD versions from 2010, 2011, 2012 and 2014, the ESP6500, the 1000G project, an interspecies collection of putatively benign variants developed by the SIFT-indel team, a set of variants fixed in recent evolution and a simulated set of variants developed by the CADD team [Boutet et al., 2007;Consortium, 2010;Hu and Ng, 2012;Fu et al., 2013;Hu and Ng, 2013;Stenson et al., 2014] (Table 2). To our knowledge, SIFT-indel, DDIG-in, PROVEAN and CADD are the four most widely used methods in the field. All provide an easy-to-use Web interface that can handle batch submissions: PROVEAN URL is http://provean.jcvi.org/genome_submit_2.php?species=human, SIFT-indel URL is http://sift.bii.a-star.edu.sg/www/SIFT_indels2. html, DDig-In URL is http://sparks-lab.org/ddig/, and CADD URL is http://cadd.gs.washington.edu/. The Web interfaces ensured that we applied each of these methods in the manner intended by their respective authors.

Multi-method benchmark set
To

Performance assessment
SIFT-indel and DDIG-in provide a categorical classification for each example (damaging/neutral or disease/neutral) and a confidence measure. PROVEAN, CADD, and VEST-indel provide a numerical score for each example (-40 to 12.5, 1 to 99, and 0 to 1). To compare methods, PROVEAN scores were assigned to categories of damaging (<-2.5) or neutral (ࣙ-2.5) (as recommended by the authors) and VEST-indel scores were assigned to categories of pathogenic (ࣙ0.5) or benign (<0.5), which represents a majority vote of decision trees in the Random Forest classifier. CADD scaled "C-scores" were assigned to categories of deleterious (<15) or not deleterious (ࣙ15) (as recommended on their Webserver). Sensitivity (TP/(TP + FN)), specificity (TN/(TN + FP)) and balanced accuracy ((sensitivity + specificity)/2) were calculated for each method, where TP = the number correctly classified as pathogenic (

Meta-Predictors that Combine Classifications of Multiple Methods
In these Boolean expressions, each method is represented by a variable X i , which is set to TRUE when the method classifies an example as pathogenic and FALSE when the method classifies an example as benign. For combinations of two methods, candidate metapredictors were (X 1 and X 2 ) and (X 1 or X 2 ). For combinations of three methods, candidate meta-predictors (X 1 and X 2 and X 3 ), (X 1 or X 2 or X 3 ),(X 1 or X 2 or X 3 ), ((X 1 and X 2 ) or X 3 ), ((X 1 or X 2 ) and X 3 ), ((X 1 and X 3 ) or X 2 ),((X 1 or X 3 ) and X 2 ), ((X 2 and X 3 ) or X 1 ), ((X 2 or X 3 ) and X 1 ). For combinations of four methods, there are 64 possible combinations (Supp . Table S4). We used a brute-force approach and limited the number of methods in the meta-predictor to a maximum of four to avoid a combinatorial explosion. All possible four-way combinations of the five methods were explored. HUMAN

Selective Pressures on Genes and False Positive Classifications
A standard method for identifying genes under selection is Tajima's D statistic [Tajima, 1989], and for each gene harboring, a variant in the multi-method benchmark set, we computed this statistic based on its longest annotated RefSeq transcript [Pruitt et al., 2007]. If RefSeq transcripts were not available for the gene, the longest annotated Ensembl transcript [Cunningham et al., 2015] was used. These calculations were performed using SNPs in 1000 Genomes Phase 3 AFR samples and the PopGenome package in R [Pfeifer et al., 2014]. Each gene was assessed for the presence of statistically significant positive or balancing selection (P < 0.05). The PopGenome package estimates P values by simulation using Hudson's coalescent model [Hudson, 2002].

Availability
VEST-indel is freely available for non-commercial use via a batch Web service at http://cravat.us. Choose the analysis program "VEST". Results for small-size submissions (<100) are returned on average in less than five minutes. Very large submissions are supported, and results for a submission of ß1,000,000 mutations are returned on average in less than 24 hr.

Results
In protein-coding exons, in-frame indels generally have a less severe impact than frameshifts [Ng et al., 2008]. Since the biological effect of in-frame and frameshift indels is different, we chose to develop two distinct Random Forest classifiers to handle these two distinct variant types. We assessed the performance of the classifiers in three phases: (1) we estimated their sensitivity and specificity with stringent, homology-restricted 10-fold cross-validation (pathogenic class from HGMDv2014.3 [Stenson et al., 2014], benign from ESP6500 African Ancestry [Fu et al., 2013]) ( Table 2); (2) we re-estimated sensitivity and specificity on an independent test set of variants (pathogenic class from ClinVar [Landrum et al., 2013], benign Interspecies alignments [Hu and Ng, 2012;Hu and Ng, 2013]) ( Table 2) that had not been used in classifier training and had been filtered for homology overlap with the cross-validation set; (3) we re-estimated sensitivity and specificity on a second independent test set of variants (pathogenic class from new entries in HGMDv2014.4, benign from 1000 Genomes Phase III African Ancestry [Consortium, 2010] that did not overlap with any training data used by previously published methods (multi-method benchmark set). These experiments are detailed in Materials and Methods.
In our cross-validation experiments, VEST-indel achieved a sensitivity and specificity of 0.90 for in-frame indels (Table 3). Crossvalidation performance for frameshift indels was slightly lower, with a sensitivity of 0.83, specificity of 0.88, and balanced accuracy of 0.85. During the testing phase, VEST-indel maintained good performance for classifying in-frame indels, with a balanced accuracy of 0.82, and sensitivity and specificity of 0.80 and 0.85, respectively (Table 3). Frameshift prediction improved slightly from the cross-validation to testing phase. For testing, frameshift classification resulted in a balanced accuracy of 0.87, with a sensitivity of 0.89 and specificity of 0.86. Table 4 compares VEST-indel performance with that previously reported for other methods. The VEST-indel specificity of 0.90 for classifying in-frame indels was higher than that achieved by SIFT-indel (0.82) or PROVEAN (0.80 for insertions and 0.68 for  deletions); the authors of the DDIG-in method did not provide a specificity for in-frame classification. The specificity of 0.88 achieved by VEST-indel for classifying frameshift indels was 0.10 higher than that achieved by SIFT-indel and 0.16 higher than that reported for DDIG-in (Table 3); PROVEAN does not classify frameshift indels. The higher specificity achieved by VEST-indel results from a relatively low false positive rate, indicating improved ability to discriminate neutral indels from those that are pathogenic. Four of the compared methods attained reasonably high sensitivities, demonstrating the ability to identify truly pathogenic variants (Table 4). PROVEAN had the highest sensitivity for classifying in-frame variants, which was higher that realized by VEST-indel (0.90), DDIG-in (0.89), or SIFT-indel (0.81). Conversely, the SIFTindel sensitivity of 0.90 was highest among compared methods for classifying frameshift indels, with DDIG-in at 0.86 and VEST-indel at 0.83. The authors of the CADD method do not report sensitivity or specificity statistics for either in-frame or frameshift indels [Kircher et al., 2014].
Although the above-cited comparison of previously published results provides some indication of relative performance, the five methods compared in Table 4 utilized different datasets for crossvalidation and testing (see Table 2). The datasets being different, to some extent, limits the insights that can be gleaned from comparison. In addition to using the same dataset to benchmark all methods, direct and fair comparison requires that no training examples from any of the methods be present in the benchmark dataset. We constructed a multi-method benchmark set, comprising indels not present in any method's training set, by using mutations recently added to HGMD and variants from the latest release of the 1000 Genomes project (Phase 3; Table 2). These indels were new to all of the methods, and were intended to reduce overly optimistic estimates of performance. Table 4 compares performance achieved by the five methods for classifying neutral and disease-causing indels from the multimethod benchmark set. VEST-indel shows superior specificity for classifying both in-frame (0.96) and frameshift (0.95) indels. These high specificities further validate the ability of VEST-indel to accurately reject neutral variants as disease causing. All methods had reasonably high balanced accuracies for in-frame indel classification, with VEST-indel and PROVEAN yielding the highest balanced accuracy of 0.88. Of note, VEST-indel and PROVEAN achieved nearly identical balanced accuracies with approximately equal trade-offs in sensitivity and specificity (Table 3). For frameshift variants, VESTindel outperformed the other methods, having a balanced accuracy of 0.90, compared with 0.77 for DDIG-in, 0.59 for SIFT-indel, and 0.52 for CADD. In the case of DDIG-in, VEST-indel showed substantially improved sensitivity and specificity (Table 4). The dramatic gain in performance achieved by VEST-indel, relative to SIFTindel and CADD, resulted from a marked gain in specificity (0.95 vs. 0.25 for SIFT, and 0.05 for CADD); this is consistent with previous reports for SIFT-indel, which maintains good specificity when predicting protein-damaging indels, but suffers low specificity when predicting pathogenicity [Hu and Ng, 2012;Hu and Ng, 2013].
Although they are often confounded, bioinformatics prediction of protein-damaging variation is not necessarily the same as prediction of pathogenic/disease-causing variation [Capriotti et al., 2012]. A variant that reduces protein stability, function, or even results in complete loss of protein production is certainly protein damaging, but implications for health and disease will also depend on the importance of that particular protein in complex networks of interacting molecules [Capriotti et al., 2012]. Thus, although all diseasecausing variants are likely to be protein damaging, all proteindamaging variants are not invariably disease causing. A good set of predictive features should capture both protein damage and gene importance. VEST-indel uses 23 features for in-frame classification and 16 features for frameshift classification; these features were top-ranked by a Random Forest Z-score feature selection method and a greedy algorithm that maximized Random Forest ROC AUC (see Materials and Methods). Features include measures of DNA sequence conservation, DNA natural variation, gene-level annotations, transcript-level annotations, computational predictions of protein local structure, protein local regional sequence composition and protein-level annotations from UniProtKB (Supp. Table S1). Some of these features are similar to those adopted by previously published methods, whereas others have not been previously applied to indel classification. In the first category, we include DNA sequence conservation scores, indel length, indel location within the transcript, predicted highly flexible or disordered regions in protein structure, predicted solvent accessibility in protein structure, and occurrence of the indel within a repetitive (low-information) sequence. In the second category, we include a hidden-Markov modelbased score [Karchin et al., 2005] of an alanine substitution, several properties of amino-acid residue sequence in a 15-residue window around the position where the indel begins, the density of single nucleotide polymorphisms in the coding exon where the indel begins, and the number of results returned from a PubMed search of the HUGO-approved gene name harboring the indel. A feature matrix that compares features used in VEST-indel to those used by SIFT-indel, DDIG-in, PROVEAN, and CADD is provided (Supp .  Table S4).
Among all VEST-indel features, the PubMed feature was consistently the most informative, for both frameshift and in-frame variants. For this feature, the algorithm searches the title and abstract of every publication indexed by PubMed, and returns the number of publications mentioning the HUGO gene name (no aliases) in which the mutation occurs [Schuler et al., 1996;Wheeler et al., 2007]. Thus, genes that have been subject to relatively greater attention in the biomedical literature will be scored as more relevant to human health. This result is consistent with previous reports that literature mining is a useful proxy for the importance of a given gene to human health [Schuler et al., 1996;Karchin et al., 2005;Wheeler et al., 2007;Masica and Karchin, 2011;Capriotti et al., 2012;Clarke et al., 2012]. Insertion/deletion variants in less important genes may knock-down or knock-out protein activity but may not necessarily cause disease [Ng et al., 2008].
Not all variants in genes important for human health are necessarily pathogenic, and over reliance on gene-level disease relevance might ultimately decrease performance. Therefore, we compared the five tested methods on a difficult set of variants, limited to genes that contained at least one pathogenic and one benign variant (141 pathogenic and 78 benign in-frame variants in 57 genes, and 561 pathogenic and 88 benign frameshift variants in 86 genes). Using the scores from the cross-validation experiments (Table 3), we recomputed sensitivity and specificity for VEST-indel for these variants only, and also scored them using the DDIG, SIFT, PROVEAN, and CADD Webservers (Supp . Table S5). VEST-indel had the highest balanced accuracy of the five methods on this difficult set (0.78 in-frame, 0.67 frameshift).
Although VEST-indel, SIFT-indel, DDIG-in, PROVEAN, and CADD share some similarities with respect to training sets and features, we considered that they might be different enough to provide independent information about an indel of interest. Therefore, they could be combined into a meta-predictor to yield improved performance. This approach has had some success in predicting the pathogenicity of missense variants [Gonzalez-Perez and Lopez-Bigas, 2011;Frousios et al., 2013;Martelotto et al., 2014]. Using the multi-method benchmark set, we assessed the classification performance resulting from each pair, trio, or quartet of methods combined using Boolean conjunctions and disjunctions. See Supp. Tables S6 and S7 for a complete list of the tested combinations.
For in-frame classification, the combination of ((VEST-indel AND PROVEAN) OR (CADD AND DDIG-in)) yielded a substantially improved sensitivity (0.93) while retaining good specificity (0.97), when compared to VEST-indel alone (sensitivity = 0.81, specificity = 0.96), and indeed any of the methods alone (Supp .  Table S6). This result indicates that these methods are highly complementary when combined in the described fashion. Conversely, for frameshift classification, the combination of ((VEST-indel AND (SIFT-indel OR DDIG-in)) had roughly equivalent sensitivity (0.83) and specificity (0.97) to VEST-indel alone (sensitivity = 0.85, specificity = 0.95). This results because the most specific method (VESTindel) is combined using the AND operation (i.e., sensitivity could not possibly increase, nor could specificity decrease).
The strategy of classifying a variant as pathogenic if any of the classifiers predicted it to be pathogenic (i.e., combining classifiers with a Boolean OR) did not yield good results. For the in-frame classifier, HUMAN MUTATION, Vol. 37, No. 1, 28-35, 2016 the combination (VEST-indel OR SIFT-indel OR PROVEAN OR CADD) had a sensitivity of 1 but a specificity of 0.56, with balanced accuracy of 0.78. Combining four classifiers or three classifiers with the OR operator consistently yielded good sensitivity but a substantial decrease in specificity. This result is, to some extent, expected because combining classifiers with the OR operation increases the possibility of accepting a variant as pathogenic. Conversely, requiring that all classifiers agree (i.e., combining classifiers with a Boolean AND) reduces the probability of a pathogenic classification. Indeed, all meta-predictors that used only AND operators had high specificity, but low sensitivity. For example, the (VEST-indel AND SIFTindel AND CADD AND DDIG-in) meta-predictor had a specificity of 1.00 and sensitivity of 0.46. Taken together, these results highlight the benefit of developing meta-predictors that combine Boolean conjunctions and disjunctions, rather than considering only a single type of Boolean operation.
The benign examples in the multi-method benchmark set were taken from the 1000 Genomes Phase 3 samples, limited to individuals in the AFR (African) super-population [Consortium, 2010] and having MAF ࣙ 0.1. Whereas common variants are generally considered to be non-pathogenic [Tennessen et al., 2010], datasets of common variants may be contaminated by pathogenic variants if they occur in genes that are not under purifying selection [Hu and Ng, 2012]. We assessed the possibility that the multi-method benchmark set might include common pathogenic variants. If this were the case, a false positive call from one of the methods might represent the correct identification of a truly pathogenic variant. Genes not subject to purifying selection might alternatively be under positive, balancing, or relaxed (neutral) selection [Hu and Ng, 2012]. For each of VEST-indel, SIFT-indel, DDIG-in, PROVEAN, and CADD we assessed the relationship between variants under selective pressure and those called as false positives, using Fisher's exact test (twotailed, α = 0.05). None of the benign variants were under balancing selection, defined as a statistically significant (nominal P < 0.05) positive Tajima's D statistic [Tajima, 1989]. Thirteen frameshift and 56 in-frame variants were under positive selection, defined as a statistically significant (nominal P < 0.05) negative Tajima's D statistic. With the exception of a borderline P value for DDIG-in in-frame variant classification (P = 0.051), there were no statistically significant relationships between positively selected variants and variants that were called as false positives. For frameshift variants, P = 1.0 for VEST-indel, P = 0.26 for SIFT-indel, P = 0.18 for DDIG-in, and P = 0.40 for CADD. For in-frame variants, P = 0.25 for VEST-indel, P = 0.56 for SIFT-indel, P = 0.13 for PROVEAN, and P = 0.79 for CADD.
VEST-indel P values for in-frame and frameshift indels are comparable to VEST P values for missense variants and as a result, multiple variant types can be jointly prioritized. We assessed joint prioritization performance by combining variants from the VESTindel in-frame and frameshift training sets (Table 1) and variants from the VEST missense training set  (2,475 pathogenic and 1,877 benign in-frame indels; 24,478 pathogenic and 1,350 benign frameshift indels; 38,221 pathogenic and 38,221 benign missense variants). We also assessed performance in a balanced set, in which we randomly selected 1,350 pathogenic and 1,350 benign variants of each type for the combined set. VEST P values and scaled CADD scores were used to compute ROC area under the curve (AUC) as described in Materials and Methods. For the combined set, VEST and CADD achieved a similar ROC area under the curve (AUC) of 0.90 and 0.88, respectively. For the balanced set, VEST classification resulted in an AUC of 0.91 and CADD classification resulted in an AUC of 0.74.

Discussion
In this study, we sought to develop a method for predicting indel pathogenicity. This functionality is distinct from existing classifiers that were developed to predict indel impact on protein structure or function. Although clinical utility appears to be a common goal for much of bioinformatics methods development, indel pathogenicity prediction presents the challenge of distinguishing variants that affect protein structure and function from those that adversely affect health [Consortium, 2010]. Given the enrichment for protein sequence and annotation features available for algorithmic development [Boutet et al., 2007;Pollard et al., 2010], the difficulty discriminating neutral and disease-causing indels might be unsurprising. Our newly developed classifier, VEST-indel, partially addresses previous methodological limitations, and achieves high balanced accuracy even when tasked with sorting disease-associated indels from those present in the general population. In particular, VEST-indel realized substantial gains in specificity relative to existing methods, highlighting reductions in falsely classifying neutral variants as pathogenic. To realize these performance gains, VEST-indel heavily utilized a new feature that captures the known relevance of a gene to human health. This new "PubMed" feature leverages decades of community-wide biomedical research. Thus, the algorithm uses features that ultimately estimate indel impact on protein, and the PubMed feature additionally estimates the biological context of the protein. Given that poor specificity also limits the utility of methods aimed at assessing the pathogenicity of missense variants [Chan et al., 2007;Hicks et al., 2011;Thusberg et al., 2011;Shihab et al., 2013], the approach presented here might prove beneficial for variant classification in general.
The in-frame meta-predictor ((VEST-indel OR DDIG-in) AND (PROVEAN)) achieved excellent sensitivity (0.93) and specificity (0.94) when applied to our multi-method benchmark dataset. This complementarity results because the two high-specificity classifiers are combined using the OR operation, which is then combined with the high-sensitivity classifier PROVEAN, using the AND operation (see individual classifier performance, Table 4). The Boolean OR operation increases the possibility of pathogenic classification; importantly, pathogenic classification from VEST-indel and DDIGin is complementary rather than entirely overlapping, hence the increased sensitivity relative to either method alone. As expected, however, the specificity of the (VEST-indel OR DDIG-in) classifier decreased (see Table S6). Next, even though the highly sensitive PROVEAN is slightly more prone to false positives, the specificity of the meta-predictor cannot decrease owing to the unanimity required by the AND operation; on the contrary, the complementarity of true-negative calls among these three classifiers restores a high specificity. We are deliberate in this explanation because metapredictor derivation relying on a single Boolean operation type is limiting and can result in significant trade-offs in sensitivity and specificity. As our results show, taking advantage of the complementarity that can result from combining Boolean conjunctions and disjunctions can be beneficial when maximizing balanced accuracy is desired.
Our new VEST-indel method can be used in combination with VEST scoring of missense variants to yield a jointly prioritized list of both variant types. This analysis requires a single batch submission to the CRAVAT server . To our knowledge, the only other automated method available for such joint ranking is CADD. For data sets in which the number of missense variants far exceeds the number of indels, VEST and CADD have similar performance. However, when variant types (indels and missense) 34 HUMAN MUTATION, Vol. 37, No. 1, 28-35, 2016 and classes (pathogenic vs. benign) are evenly distributed, VEST significantly outperforms CADD.

Author Contributions
C.D., D.L.M., and R.K. conceived of the study and designed the experiments. C.D. performed the experiments and wrote the software prototype. D.M.G., R. Kim, and M.R. wrote the production version of the software. P.D.S. and D.N.C. obtained and compiled the HGMD data and provided critical feedback about the manuscript. C.D., D.L.M., and R.K. wrote the manuscript.


Differences in Activation Area Within Brodmann Area 2 Caused by Pressure Stimuli on Fingers and Joints

Abstract In this study, a constant pressure stimulus was applied on the 3 joints (first [p1], second [p2], and third [p3] joints) of 4 fingers (index, middle, ring, and little fingers), and the activation areas within Brodmann area 2 (BA 2) were compared for these different fingers and joints by using functional magnetic resonance imaging. Eight healthy male college students (25.4 ± 1.32 years) participated in the study. Each session was composed of 3 blocks, and each block was composed of a Control phase (30 seconds) and a Pressure phase (30 seconds). No pressure stimulus was applied in the Control phase, during which the subjects would simply lay comfortably with their eyes closed. In the Pressure phase, a pressure stimulus was applied onto one of the joints of the selected finger. For each finger and joint, BA 2 areas activated by the pressure stimulus were extracted by the region of interest method. There was a significant difference in the activation areas for the different fingers (P = .042) as well as for the different joints (P = .050). The activation area decreased in the order of the little, index, and middle fingers, as well as in the order of p1, p3, and p2.

INTRODUCTION
M any studies have examined, using functional magnetic resonance imaging (fMRI), the neurological changes in the primary somatosensory area (S1) caused by the application of various tactile stimuli. [1][2][3][4][5][6][7] Such studies provide important insights into not only the normal organization of sensory cortices, but also into how neural systems are affected by disease or damage and how they may change over time during remodeling and with intervention. [8][9][10] Data relating to the application of various tactile stimuli may also be useful in longitudinal studies designed to monitor cerebral plasticity and reorganization, for example in sensorimotor recovery after neurosurgery, cerebral ischemia, 11 or in understanding mechanisms of supraspinal pain processing. 12 The S1 area can be divided into Brodmann areas (BA) 1, 2, and 3. BA 3 is activated by vibrations, pressure, and most other tactile stimuli; BA 2 is activated by pressure, joint position, and complex touch; and BA 1 usually responds to vibrotactile stimuli. 13 The above-mentioned studies examined changes in activation areas within the somatosensory area, especially in BA 3, in response to vibration stimuli on fingers. 1 -7,14 Differences in the activation areas within BA 3 were observed across joints 7,14 and across fingers when a vibration stimulus was applied to the first joint. [1][2][3][4][5][6] These studies point to the consistent interest in BA 3 activation areas in relation to application of vibration stimuli to fingers and joints.
Other studies have focused on how pressure stimuli lead to changes in the activated area of the somatosensory area. For example, some studies observed changes in BA 3 activation areas upon application of a pressure stimulus to the first joints of the fingers. [15][16][17] Specifically, van Westen et al 17 found that when the pressure stimulus was applied to the first joint of the little finger, the activation area was larger than when the pressure stimulus was applied on the middle finger. Overduin and Servos 15,16 found that the activation area triggered by the index finger was larger than that triggered by the little finger when the same pressure stimulus was applied to the first joints of both these fingers. According to previous studies, activation areas within BA 3 decrease in the order of the index, little, and middle fingers when a pressure stimulus is applied to the first joints of each of these fingers. Although some differences have been reported in BA 3 activation areas triggered by pressure stimuli applied to different fingers, no study seems to have addressed differences in activation areas for different joints. Moreover, previous researchers focused only on BA 3, referred to as the ''primary somatosensory cortex,'' and not on BA 2, which is known to respond well to pressure stimuli. As such, further research is needed in order to investigate activation areas in the case of pressure stimuli applied to fingers. In other words, it is important to conduct research on differences in activation areas for pressure stimuli applied to different joints as well as on those for pressure stimuli applied to different fingers and to expand the focus from BA 3 to BA 2, which has high sensitivity to pressure stimuli.
The present study was conducted focusing on the following 2 points in order to overcome the limitations of preceding studies. First, to date, preceding studies have compared activation areas in response to vibration and pressure stimuli applied to the first joints of different fingers, whereas in the present study, we have carried out the comparison between fingers as well as between joints. Second, most preceding studies focused on BA 3 of the somatosensory area when vibration and pressure stimuli were applied to fingers, whereas the present study has focused on BA 2, which has been reported as the most sensitive to pressure stimuli.
Accordingly, in this study, pressure stimuli were applied to the 3 joints (first [p1], second [p2], and third [p3] joints) of 4 fingers (index, middle, ring, and little finger) and the resulting activation areas within BA 2 were compared between these different fingers and different joints.

METHODS
To date, there have been no reports on the differences in activation patterns in response to tactile stimuli applied to different joints on the fingers. As a first step toward accurately assessing these differences, a preliminary experiment was first performed to assess the contribution of sex with a higher sensitivity because results of a preceding study showed differences in sensitivity to tactile stimuli by sex. 18 Specifically, a preliminary experiment including 5 men (23.2 AE 1.2 years) and 5 women (22.6 AE 0.7 years) was performed to determine their ability to differentiate between pressure stimuli by joint. A pressure stimulus (8.5 psi) was applied to the 3 joints of the index finger with the intention to study the contribution of sex with a higher sensitivity. We found that men more accurately differentiated pressure stimuli by joint than women did. Therefore, 8 healthy male college students (25.4 AE 1.32 years) participated in the main study as subjects. None of them reported having a history of psychiatric or neurological disorders. The overall procedure was explained to all of them, and they provided informed consent for the procedure. All experimental procedures were approved by and performed under the regulations of the Institutional Review Committee of Korea University (KU-IRB-11-46-A-1).
Pressure stimuli were applied using an MR-compatible pressure stimulator. 19 This method stimulates finger joints by injecting air into a noninvasive blood pressure cuff (M1866A, Philips, The Netherlands) ( Figure 1A). The air pressure generated from the air pump is delivered to the cuff, which measures 6.4 cm Â 2.5 cm, via a 7 m air tube. E-Prime software (Psychology Software Tools, Inc. Sharpsburg, PA) was used to control parameters such as the pressure intensity and time. The functional images were obtained while applying a constant pressure (8.5psi) on a joint (p1, p2, and p3) of each of the 4 right fingers (index, middle, ring, and little fingers).
Three pressure stimulators were attached to the 3 joints (p1, p2, and p3) of an arbitrarily selected finger (index, middle, ring, or little finger) ( Figure 1A). Figure 1B shows the paradigm for the fMRI experiments. A session consisted of 3 blocks, and each block was composed of a Control phase (30 seconds) and a Pressure phase (30 seconds). No stimulus was applied in the Control phase, during which the subjects would simply lie comfortably with their eyes closed for 30 seconds. In the Pressure phase, a 30-second pressure stimulus was applied onto one of the joints (p1, p2, or p3) of the selected finger (index, middle, ring, or little finger). Because there were 3 Pressure phases in the 1 session, all 3 joints were subjected to pressure stimuli. Each session was repeated thrice for each finger. During these 3 sessions, the 3 joints of the selected finger were subject to pressure stimuli for a total of 3 times. The remaining 3 fingers were subjected to the same procedures. The order of providing the stimuli for each finger and joint was counterbalanced, and a 5 min break was given between sessions. All subjects participated in 12 sessions (3 sessions/finger T 4 fingers ¼ 12 sessions). They wore headsets and kept their eyes closed to eliminate audiovisual influences. The subjects were asked to minimize hand and head movements.
Images were scanned with a 3-T MRI system (Magnetom TrioTim, Siemens Medical Systems, Erlangen, Germany) using a standard 32-channel head coil. Anatomical images were obtained using a T1-weighted 3D-MPRAGE sequence with a repetition time (TR) of 1900 ms, echo time (TE) of 2.48 ms, flip angle of 98, field of view (FOV) of 200 mm, and spatial resolution of 0.8 T 0.8 T 1 mm 3 . Functional images were obtained using a T2Ã-weighted gradient echo-planar imaging sequence with TR of 3000 ms, TE of 30 ms, flip angle of 908, FOV of 192 mm, slice thickness of 2 mm, and in-plane resolution of 1.5 T 1.5 mm 2 .
The fMRI data were analyzed using SPM 8 (Wellcome Department of Cognitive Neurology, London, UK). All functional images were aligned with the anatomical images from the study using affine transformation routines built into SPM 8. The realigned scans were coregistered to a subject's anatomical images obtained within each session and normalized to SPM8's template image that uses the space defined by the Montreal Neurological Institute. Motion correction was performed using sinc interpolation. Time-series data were filtered with a 240second high-pass filter to remove artifacts because of cardiorespiratory and other cyclical influences. The functional map was smoothed with a 3 mm isotropic Gaussian kernel before statistical analysis. Statistical analysis was performed both individually and as a group using the general linear model and the theory of Gaussian random fields incorporated into SPM 8. Statistical parametric maps with the t-statistic were computed. The analysis of data for individual subjects was performed at a significance threshold of P < 0.05 corrected for multiple comparisons. Furthermore, random effect analysis was used to compute group activation. The region of interest method was used to observe activation areas only within BA 2 upon application of pressure stimuli. The subtraction method (Pressure phase-Control phase) was used to extract from BA 2 the number of activated voxels for stimuli on each finger and joint by subject, and activation areas were then calculated (activation area [mm 2 ] ¼ number of activated voxels T 1.5 mm T 1.5 mm). Statistical differences in activation areas for the different fingers and different joints were determined by performing repeatedmeasures two-way ANOVA (PASW Statistics 18). Moreover, in order to examine if there was a statistically significant difference between activated BA 2 areas activated by stimuli to each joint (p1, p2, or p3) by finger and activated BA 2 areas in response to stimuli by joint in each finger (index, middle, and little fingers), repeated-measures one way ANOVA (PASW Statistics 18) was performed for each area. In addition, a pairwise analysis revealed significant differences in activation areas for the different fingers and different joints. Table 1 shows BA 2 activation areas (mean AE SE mm 2 ) following application of pressure stimuli to the different fingers and joints for each subject. In addition, the group analysis results of all subjects in response to pressure stimuli are shown in Figure 2 according to finger and joint. In the case of stimuli applied to the joints of the ring finger, activation areas in BA 2 were observed only for 3 out of 8 subjects. This result was regarded as insignificant and was thus excluded from further analysis. Repeated-measures two-way ANOVA revealed significant differences in activation areas between different fingers (P ¼ 0.042) and between different joints (P ¼ 0.050) ( Table 2).

RESULTS
With regard to the stimulated finger, the activation area in BA 2 was largest following stimulation of the little finger, followed by the index finger and finally the middle finger. Statistical analyses of activated BA2 areas by finger identified a significant difference by finger in case of p1 (P < 0.001), whereas no significant differences by finger were identified in case of p2 and p3. Pairwise analysis of p1 found significant differences between the index finger and the middle finger (P ¼ 0.032), between the index finger and the little finger (P ¼ 0.040), and between the middle finger and the little finger (P ¼ 0.001) ( Figure 3A).
With regard to the stimulated joint, the activation area in BA 2 was largest in response to stimulation of p1, followed by p3 and finally p2. Statistical analysis of activated BA 2 areas by joint showed a significant difference by joint in the index finger (P ¼ 0.049) and the little finger (P ¼ 0.001), whereas the middle finger exhibited no significant difference by joint. Pairwise analysis identified significant differences between p2 and p3 (P ¼ 0.002) in the index finger, and between p1 and p2 (P ¼ 0.004), and between p2 and p3 (P ¼ 0.028) in the little finger ( Figure 3B).

DISCUSSION
A constant pressure stimulus was applied on the 3 joints (p1, p2, and p3) of 4 fingers (index, middle, ring, and little finger), and the activation areas within BA 2 for the different fingers and different joints were compared using fMRI.
In our study, it was found that the activation area within BA 2 was largest when the stimulus was applied to the little finger, followed by the index finger and finally the middle finger. Our results are consistent with those of a study by Peters et al, 20 who determined the tactile sensitivity of the first joints of the 5 fingers through a subjective assessment and reported higher sensitivity for smaller fingers. This is because the smaller the finger, the more elaborate the distribution of the sensory receptors, which in turn leads to higher tactile sensitivity. The little finger, being the smallest of the 5 fingers, has the highest sensitivity. However, our result is inconsistent with the results of some previous work, [15][16][17] in which the activation area (BA 3) decreased in the order of the index, little, and middle fingers. The discrepancies between the results of our study and those of previous studies 15-17 on the activation area according to finger stimuli can be explained by one of the following reasons. In the previous studies, it was difficult to observe high levels of sensitivity to specific stimuli (pressure), as the focus was on BA 3, which is an area that responds to tactile stimuli such as vibrations and pressure. However, the present study observed the activation area within BA 2, which is an area that usually responds to pressure stimuli, and this is likely to have produced results with higher sensitivity. The differences could also have arisen from the choice of experimental method (the type of stimuli, strength, time, etc.). Nevertheless, further studies are required to determine the exact causes.  Because past research involved providing stimuli only to the first joint of each finger, it was impossible to compare the differences in activation areas between different joints. In contrast, the present study found a decrease in the activation area within BA 2 in the order of p1, p3, and p2. Weinstein 21 applied pressure stimuli of various strengths to the joints of 5 fingers and measured the tactile sensitivity of subjects based on a subjective assessment. The p1 was found to have the highest sensitivity, followed by the p3 and finally the p2; this result is in agreement with the results of our study. This implies that the activation area within BA 2 accurately reflects the tactile sensitivity to pressure stimuli on joints.
In a study of tactile sensitivity, Weinstein 21 found that the ring finger was less sensitive than the other 4 fingers. This may explain why we observed some inconsistencies when the pressure stimulus was applied to the ring finger in the present study.
In conclusion, the largest activation area within BA 2 was associated with the little finger in the case of pressure stimuli applied to different fingers, and with the p1 in the case of pressure stimuli applied to different joints. The results of previous work conducted with a focus on tactile sensitivity 20,21 were verified from a neurological perspective. This study showed that differences in tactile sensitivity to pressure stimuli on different fingers and joints can be reflected in changes in the activation area within BA 2. These findings are expected to serve as a basis for neurological research on the effects of pressure stimuli provided to fingers and joints on the S1 (BA 2).
The results of the present study could be referred to while carrying out a comparison between normal subjects and abnormal subjects for studies of somatosensory activation in clinical trials (eg, studies on cerebral plasticity and reorganization following neurosurgery or ischemic brain injury). In general, pressure stimuli have been used mainly to diagnose patients who have neurological disorders that affect tactile sensation. However, tactile stimuli have been applied only to the first joint of the index finger for the diagnosis of patients. According to the results of the present study, the p1 of the little finger, which showed the largest activation in response to pressure stimuli among the 3 joints of the 4 fingers, would be useful in the diagnosis of neurological disorders in the somatosensory area. However, it is necessary to conduct various further studies to verify this potential clinical application.
In the present study, we applied pressure stimuli to each of the 3 joints of the 4 fingers, compared the activation areas in BA 2 in the S1 area by finger and joint, respectively, and then presented the results. Despite not being presented in the results, data on the activation areas in BA 1 and BA 3 were also extracted by finger and by joint. However, activation areas appeared in only 2 out of 8 subjects in total. In other words, the BA 2 area that dominantly receives information of pressure tactile stimuli clearly showed activation by finger and by joint, whereas BA 1 and BA 3 did not show significant activation. It is possible that these results may have been influenced by the small number of subjects in the present study. It is necessary to increase the sample size to analyze activation patterns not only in the BA 2 area, but also in BA 1 and BA 3, and we also need to increase the sample size for a comparison with BA 2. In addition, it is necessary to study differences by sex and by age, and to perform a comparative study between normal groups and patient groups, which would create a more solid basis for clinical application.