text:
A DNA break inducer activates the anticodon nuclease RloC and the adaptive immunity in Acinetobacter baylyi ADP1

Double-stranded DNA breaks (DSB) cause bacteria to augment expression of DNA repair and various stress response proteins. A puzzling exception educes the anticodon nuclease (ACNase) RloC, which resembles the DSB responder Rad50 and the antiviral, translation-disabling ACNase PrrC. While PrrC's ACNase is regulated by a DNA restriction-modification (R-M) protein and a phage anti-DNA restriction peptide, RloC has an internal ACNase switch comprising a putative DSB sensor and coupled ATPase. Further exploration of RloC's controls revealed, first, that its ACNase is stabilized by the activating DNA and hydrolysed nucleotide. Second, DSB inducers activated RloC's ACNase in heterologous contexts as well as in a natural host, even when R-M deficient. Third, the DSB-induced activation of the indigenous RloC led to partial and temporary disruption of tRNAGlu and tRNAGln. Lastly, accumulation of CRISPR-derived RNA that occurred in parallel raises the possibility that the adaptive immunity and RloC provide the genotoxicated host with complementary protection from impending infections.

PrrC's best characterized ortholog was encountered in a rare Escherichia coli strain restrictive to phage T4 mutants lacking 3 0 -phosphatase/5 0 -polynucleotide kinase (Pnk) and/or RNA ligase 1 (Rnl1) (12). PrrC's ACNase is normally silenced in this host by an associated type Ic DNA restriction-modification (R-M) protein (13,14) and turned on by a T4-encoded anti-DNA restriction peptide (15). The resulting incision of the tRNA Lys anticodon loop could disable T4 late protein synthesis and contain the infection but T4's tRNA repair proteins Pnk and Rnl1 reverse the damage (16,17). This host/phage survival cascade is probably shared by most PrrC-encoding bacteria. Namely, prrC is invariably linked to a type Ic R-M (hsd) locus. Moreover, PrrC orthologs looked at exhibited ACNase activity unless linked to a cryptic R-M system (1,7,18,19). The notion that Pnk and Rnl1 evolved as ACNase antidotes (20,21) is supported by their ubiquity among T4-like phage of PrrC-encoding bacterial clades but absence from T4-like cyanophages that are not expected to encounter an ACNase (2). Pnk also figures in mRNA decay during phage T4 development (22).
PrrC's ABC-ATPase N-domain mediates the activation of the C-proximal ACNase by hydrolysing GTP and stabilizes the active ACNase by avid binding of dTTP, which accumulates in the T4-infected cell (7,23). This dual nucleotide specificity sets PrrC apart from other ABC-ATPase containing proteins, including RloC. PrrC's C-domain harbours residues implicated in tRNA Lys recognition (24)(25)(26)(27) and a putative catalytic ACNase triad (7) shared by RloC (1). PrrC is thought to act as a dimer of dimers whose C-domains dimerize in parallel and the N-domains head-to-tail as do typical ABC-ATPases (7,27). Another view is that PrrC is a dimer whose ACNase domains do not interface (18).
RloC shares PrrC's overall organization into ATPase and ACNase domains but differs in several key respects. First, rloC is only rarely linked to an R-M system although in trans interactions of non-linked orthologs have not been excluded (1). Second, isolated PrrC has overt ACNase activity indifferent to nucleotide hydrolysis but highly unstable without dTTP (7,28). In contrast, purified RloC encoded by the thermophile Geobacillus kaustophilus (GkaRloC) exhibits ACNase activity only when turned on by its ATPase in the presence of linear DNA. Once activated, GkaRloC's ACNase can be stabilized by the non-hydrolysable ATP analogue AMPPNP but not by dTTP (2). It has been proposed that the advent of RloC's internal ACNase switch rendered unnecessary the reliance on an external R-M silencing partner. This, in turn, could facilitate RloC's lateral gene transfer and account for its broader distribution among bacteria, compared with PrrC (2). Third, while PrrC only incises its tRNA substrate, GkaRloC excises the wobble nucleotide, a lesion expected to encumber phageinduced tRNA repair and likewise contribute to RloC's broader distribution (1). Fourth, a coiled-coil/zinc-hook insert (CC/ZH) likens RloC's ATPase domain to the universal DSB sensor/effector Rad50 (5). The CC/ZH stretch of Rad50 folds back into an antiparallel coiled-coil bundle and this structure protrudes from the ATPase headdomain with the ZH motif Cys-X-X-Cys at its apex. Two ZH apices join by coordinating Zn + + to the four Cys residues of the dimerization interface. Such joints may form within a free Rad50 dimer or between two DSB-bridging dimers. In the latter state, Rad50 directs the associated Mre11 DNase to initiate DNA end resection, a key step in DSB repair (29)(30)(31)(32). A regulatory role of GkaRloC's CC/ZH is suggested by ZH mutations that uncouple its ACNase and DNA-dependent ATPase and the activation of this ACNase by linear but not circular DNA (1,2).
Further exploration of RloC's regulation reported here revealed the following. First, once turned on, RloC's ACNase is stabilized by its activating DNA and the nucleotide it hydrolysed. Second, exposure to a DSB inducer activates RloC's ACNase both in heterologous contexts and within a natural host, even when R-M mutated. Third, tRNA Glu and tRNA Gln are natural RloC substrates. Lastly, the DSB-induced activation of the indigenous ACNase was accompanied by accumulation of CRISPR-derived RNA (crRNA) (33), hinting that RloC can complement the adaptive immunity in defending the genotoxicated host from impending infections.

Bacterial strains and plasmids
Acinetobacter baylyi ADP1 was obtained from Eliora Z. Ron, Tel Aviv University, and the DrloC, DhsdR, DhsdM and DhsdS alleles of this strain were obtained from Veronique de Berardinis, Genoscope. Expression plasmids encoding wild-type GkaRloC pGkaRloC-L-His 6 and its ACNase-null derivative E696A were previously described (1).

Plasmid construction and expression of RloC proteins
Acinetobacter baylyi ADP1 orf 0152 encoding the 826aa RloC protein (NCBI accession YP_044948) was amplified by polymerase chain reaction (PCR) from genomic DNA using Pfu DNA polymerase (Stratagene). The PCR primers used introduced an NdeI restriction site at the start codon and an AgeI site at the C-end to fuse the orf via a flexible linker to a His 6 tag, as in plasmid pGkaRloC-L-His6 (1). To this end, the NdeI-and AgeIdigested PCR product replaced the GkaRloC orf of pGkaRloC-L-His6 to yield pAbaRloC-L-His6. Point mutations were introduced by Quick Change (34) and the GkaRloC Á152-478 deletion created by SLIM (35). The RloC constructs were transformed into E. coli XL-Blue 1 (Stratagene) and then into E. coli RosettaTM (DE3)pLys (Novagen). The latter transformants were grown at 37 C to 0.6 OD 600 in Luria-Bertani (LB) medium containing 100 mg ml À1 ampicilin and 34 mg ml À1 chloramphenicol. ACNase protein expression was induced by 0.05-0.1 mM isopropyl-b-D-thiogalacto-pyranoside (IPTG). Wild-type and Á152-478 GkaRloC proteins were isolated as described (2).

ACNase assays
In vivo ACNase activity was monitored as follows. Cellular RNA was phenol extracted and separated by denaturing gel electrophoresis, as such or after any 5 0 -OH termini had been phosphorylated by T4 Pnk and [g-32 P]ATP (24). ACNase products were visualized then by respective ethidium bromide staining and autoradiography. The in vitro assays of GkaRloC used as a standard substrate a synthetic ASL oligoribonucleotide that matched in sequence G. kaustophilus tRNA Glu and was 2 0 -O methylated 5 0 to the wobble base (7-2 0 Om-Glu-ASL, Figure 1A). Other synthetic in vitro substrates were the parental unmodified Glu-ASL and a derivative 2 0 -O methylated 3 0 to the wobble base (8-2 0 Om-Glu-ASL). These oligonucleotides were either [5 0 -32 P] labelled or 3 0 -end labelled using [5 0 -32 P]pCp and T4 RNA ligase 1 (36). The ACNase activation mixture (10 ml) contained 10-100 ng of the indicated GkaRloC allele, 5 ng/ml SmaIlinearized pUC19 DNA, 0.5 mM ATP, 10 mM MgCl 2 , 5 mM dithiothreitol, 70 mM Tris-HCl buffer, pH 7.5. Following incubation for 20 min at 25 C, 5-10 nM of the ASL substrate was added and the incubation at 25 C continued. The reaction was stopped with 2 volumes of 10 M urea, 0.025% each of xylene cyanoll and bromphenol blue. The products were separated by denaturing gel electrophoresis and autoradiographed.

Identification of NAL-induced RNA fragments
Acinetobacter baylyi ADP1 culture in LB medium containing 10 mg/ml ampicilin and 3.4 mg/ml chloramphenicol was grown by shaking at 25-30 C. At a density of 0.6 OD, the culture was supplemented with NAL at 30 mg/ml and the incubation continued at 25 C for 4 h. The cells were harvested then, their RNA phenol extracted and any 5 0 -OH cleavage termini in it labelled at 3000 Ci/mmol using T4 Pnk and [g-32 P] ATP. Following separation by denaturing gel electrophoresis, the labelled products were extracted and partially digested by mild alkaline hydrolysis or by the respective G or pyrimidine-specific RNases T1 and BC. The digests were separated by denaturing gel electrophoresis and autoradiographed, essentially as described (37). Labelled 5 0 -end groups released by nuclease P1 digestion were identified by 2D thin layer chromatography (38). The assumed identities of the RNA products of interest were ascertained by splint ligation essentially as described (39). Briefly, scaled up RNA preparations were weakly 5 0 -end labelled (1 Ci/mmol) to identify and isolate the species of interest. In the case of the RloC-independent product X suspected to be crRNA, the Pnk used for labelling also removed 3 0 cyclic phosphate (40) that would prevent ligation to the strongly labelled (3000 Ci/mmol) DNA oligonucleotide probe. Bridging oligonucleotide templates used to juxtapose the ligation partners matched the predicted 3 0 -proximal sequence of the RNA and the sequence of the probe. After T4 DNA ligase-mediated ligation any free 5 0 -P groups were removed by calf intestinal alkaline phosphatase. The bridging templates and ligation oligonucleotide probe used are listed in Supplementary  Table S1. Northern analysis of RNA extracted at different NAL exposure times entailed prior phosphorylation of 5 0 -OH to enable EDC-mediated covalent linkage to the nylon membrane (41). The blot was probed with a 5 0 -32 P labelled oligonucleotide matching A. baylyi ADP1 tRNA Glu residues 35-72 (Supplementary Table S1).

GkaRloC's ACNase is stabilized by its activating DNA and hydrolysed nucleotide
Thermal stability has rendered GkaRloC a convenient in vitro paradigm, superior to mesophilic orthologs. Previous data revealed that the activation of its ACNase by ATP hydrolysis requires the presence of linear DNA (2). To determine if the requisite DNA also stabilizes the activated ACNase it was degraded after the activation by DNase I. This inhibited the ACNase almost as severely as adding DNase I ahead of the activation ( Figure 1B, lanes 2-5, 10-13, [14][15][16][17]. To determine if this inhibition was due to the lack of the DNA or formation of the DNA degradation products we added to the activated ACNase a DNase I-degraded activating DNA dose. This addition had no effect (lanes 6-9), suggesting that continued ACNase activity required that GkaRloC remain tethered to the activating DNA. The ACNase substrate used 7-2 0 -Om-Glu-ASL ( Figure 7A) is a synthetic anticodonstem-loop matching in sequence G. kaustophilus tRNA Glu and containing a 2 0 -O methyl group intended to prevent the excising cleavage 5 0 to the wobble base. Using it simplified the determination of the incision rate. Surprisingly, the parental unmodified Glu-ASL was incised mainly 5 0 to the wobble base, yielding a product similar to that formed exclusively with a derivative whose 3 0 linkage was protected (8-2 0 -Om-Glu-ASL). The preferential 5 0 incision occurred whether Glu-ASL was labelled at the 5 0 -or 3 0 -end (Supplementary Figure S1A  Mutating critical ABC-ATPase residues abolishes GkaRloC's ACNase activity. Yet, the anticipated DNA- GkaRloC's ACNase was activated for 20 min in the presence of ATP and DNA followed by further pre-incubation with DNase buffer (lanes 1-5), or DNase I buffer containing purified oligonucleotides formed by a DNase I digestion of an activating DNA dose (lanes [6][7][8][9] or DNase I (lanes 10-13). In lanes 14-17 DNase I was included in the activation mixture. ACNase activity was subsequently assayed as detailed in 'Materials and Methods' section. DDR, DNase-I digest residue; After or During, DNase I added after the activation or during the activation, respectively; ASL, [5 0 -32 P]7-2'-Om-Glu-ASL; 8mer, labelled cleavage product.
dependent ATPase activity of wild-type GkaRloC could not be detected over the mutant background. This discrepancy and ability of AMPPNP to both inhibit the activation and stabilize the activated ACNase have suggested that the activation entails a single cycle of ATP binding and hydrolysis (2). To determine if AMPPNP emulates in its protective effect the nucleotide hydrolysed during the activation, ADP was included in the ACNase activation mixture or added after the activation. This protected the ACNase in either case (Figure 2A and B), but in the first, the protection sharply declined above $200 mM ADP, whereas in the second, it kept increasing ( Figure 2C). This difference was ascribed to an inhibitory ADP/ATP competition during the activation and ability of ADP to stabilize the activated ACNase, similar to AMPPNP. ADP's protective effect became apparent only after sufficient dilution of the ACNase protein ( Figure 2D). This signified perhaps synergistic effects of GkaRloC's oligomerization and the tethered DNA on the tendency of the hydrolysed nucleotide to dissociate from the protein.

DSB inducers activate RloC ACNases in a heterologous context
Because circular DNA fails to activate GkaRloC's ACNase, DSB have been implicated in the in vivo activation (2). To examine this premise, we determined if the DSB-inducing inter-strand crosslinker MMC (42) will activate plasmid-borne GkaRloC expressed in E. coli from the T7-Lac promoter. As controls served the ACNase-null mutant GkaRloC E696A (1) and the deregulated CCtruncated mutant GkaRloC Á152À478 that featured in vitro a constitutive-ACNase phenotype (Supplementary Figure  S2). In vivo ACNase activity was monitored by extracting RNA from the untreated or MMC-treated cells before or after inducing the expression of the ACNase protein with IPTG. The extracted RNA was separated by denaturing gel electrophoresis and stained with ethidium bromide. ACNase activity was indicated by the formation of the typical tRNA fragments of $42 and $33 nt and a $52 nt fragment thought to originate from a tRNA  9-12). The basal ACNase activity of wildtype GkaRloC is ascribed to partially degraded forms that could behave like the constitutive mutant and/or spontaneous DSB.
Weaker basal ACNase activity and a more pronounced response to MMC were observed when GkaRloC was replaced by the mesophilic A. baylyi ADP1 RloC (AbaRloC) ( Figure 4A). In this case ACNase, cleavage products were not detected in the absence of both IPTG and MMC (lane 1). They appeared with IPTG alone (lane 2) but at a level $4-fold lower than with MMC alone (lane 3) and $30-fold lower than with both IPTG and MMC (lane 4). Another DSB-inducing agent, the DNA gyrase inhibitor NAL (43) activated AbaRloC's ACNase similarly ( Figure. 4B, lane 2 versus 3 and 8 versus 9). In contrast, the overt ACNase of the ZH mutant AbaRloC C314G was indifferent to either DSB inducer (lanes [10][11][12]. This behavior could be demonstrated only when AbaRloC C314G 's expression was induced, possibly due to its weaker ACNase activity and/or stability.

DSB-induced activation of indigenous RloC
We chose to investigate RloC's natural activation using A. baylyi ADP1 since this soil bacterium is easy to grow and genetically manipulate (44). Moreover, a library of single gene deletion mutants exists for all its dispensable genes (45), including rloC (ACIAD0152) and the nonlinked hsdR, hsdM or hsdS genes of the respective R-M restriction, modification or specificity subunit (ACIAD343-2). However, an infectious phage is not known yet. Prolonged exposure of this host to a high MMC dose (1 mg/ml) failed to activate the indigenous RloC ACNase, consistent with its resistance to MMC (46). Unaware of this caveat we have suggested that DSB only set the indigenous AbaRloC for activation by an added, possibly phage-induced signal (2). However, exposing A. baylyi ADP1 to the toxic DSB inducer NAL sufficed to activate the indigenous RloC ACNase. This was inferred from the formation of ACNase-like cleavage products in the wild-type ( Figure 5A, lanes 1 and 2) but not DrloC allele (lanes 3 and 4). Interestingly, all three R-M mutants featured the wild-type phenotype (shown in Figure 5B for hsdR and hsdM), indicating that the R-M protein neither silences AbRloC's ACNase nor mediates its DSB-induced activation.
The level of the NAL-induced cleavage products formed by AbaRloC in its natural host relative to the total RNA was far lower than observed in the heterologous E. coli context, even in the absence of IPTG. This often impeded the detection of the ACNase cleavage products over the background of non-specific RNA fragments. Monitoring the indigenous ACNase activity was facilitated by 5 0 -end labelling, within the total RNA preparation, of the 3 0 tRNA fragments the ACNase generated. Separating the labelled RNA on a long gel resolved the RloC-specific fragments into two bands designated 'a' and 'b' ( Figure 5C). NAL induced in A. baylyi ADP1 also RloC-independent small RNA species. The most prominent of them named X was also 5 0 -end labelled, indicating that it originated by transesterifying cleavage of a larger  Figure 3. (B) Similar to (A) except that the cells encoded the indicated AbaRloC alleles and were exposed also to 30 mg/ml NAL. species. Its identification as crRNA (33) is described in a later section.

AbaRloC targets tRNA Glu and tRNA Gln
To identify the RNA species targeted by the NALactivated, indigenous AbaRloC, the 5 0 -32 P-labelled fragments 'a' and 'b' ( Figure 5C) were partially digested by base specific RNases. The product pattern obtained suggested that fragment 'a' originated from tRNA Glu and fragment 'b' from tRNA Gln (Supplementary Figure S3A and B). Complete digestion of either with nuclease P1 released mainly labelled pU (Supplementary Figure  S3C), commensurate with RloC-mediated incision of the original tRNA species 3 0 to the wobble base. The tRNA Glu and tRNA Gln assignments were ascertained then by splint ligation (39) using weakly labelled (1 Ci/ mmol) fragments 'a' or 'b' and a strongly labelled ligation partner (3000 Ci/mmol). Ligation of fragment 'a' to this probe was enabled by a bridging oligonucleotides matching a 3 0 proximal portion of tRNA Glu (Figure 6, lane 3) but not by counterparts matching Asp or Gln specific tRNAs. The weak signal obtained with the  Table S1). Following phosphatase treatment the ligation products were separated by denaturing gel electrophoresis and monitored by autoradiography. The splint ligation set-ups shown below identified the tRNA Glu fragment 'a' (left) and the tRNA Gln fragment 'b' (right). The non-specific NAL-induced product X and Xylene cyanoll were used as respective 60 and 30 nt markers.  Figure 3. X indicates the RloC-independent product accumulating in the presence of NAL. (B) The indicated A. baylyi ADP1 alleles were examined for NAL-induced ACNase activity essentially as in panel A. (C) 5 0 -OH termini in RNA extracted from the indicated A. baylyi ADP1 alleles exposed to NAL were radiolabelled, the RNA separated then on a long gel in which the 3 0 fragments generated by RloC were resolved in two major bands designated 'a' and 'b'. The broken lines between panels B and C match the stained fragments with their 5 0 -end labelled counterparts.
tRNA Arg specific probe could reflect some cross hybridization due to shared sequences but leaves open the possibility that this species was also cleaved. Fragment 'b' was ligated to the probe only with the tRNA Gln specific bridge (lane 6).

Dynamics of RNA cleavage products induced by NAL in A. baylyi ADP1
The dynamics of 5 0 -OH containing RNA fragments formed during A. baylyi ADP1's exposure to NAL were followed by monitoring the appearance and eventual decay of their in vitro labelled derivatives. The RloC specific fragments were detected in a logarithmic culture grown at 25 C only $2 h after the exposure to 30 mg/ml NAL. They reached peak intensity after $4 h and declined during longer incubation. In contrast, the RloC-independent product X accumulated ( Figure 7A). To determine if the formation of any of these products required new protein synthesis the culture was exposed to both NAL and the seryl-tRNA synthase inhibitor and stringent response inducer SH (47). SH did not prevent the NALinduced activation of AbaRloC's ACNase, causing rather the accumulation of its cleavage products ( Figure 7B). This accumulation could be ascribed to enhancement of AbaRloC's activation due to exacerbated DNA damage in the absence of newly synthesized DNA repair proteins. Alternatively, or in addition the accumulation of the ACNase products could be due to the lack of an RNase responsible for their further decay. On the other hand, the accumulation of X was severely inhibited by SH (Figure7B). Thus, NAL must have activated pre-existing AbaRloC molecules, while the enhanced production of X and perhaps also the decay RloC's cleavage products depended on new protein synthesis. It is also noteworthy that the NAL dose used in this experiment was 60 mg/ml, 2-fold higher than the standard. This increase hastened AbaRloC's activation and decay of its cleavage products as well as the accumulation of X ( Figure 7B), possibly owing to exacerbation of the DNA damage. However, SH exerted the same effects at either NAL level (not shown).
Northern analysis of RNA samples extracted during an extended exposure to NAL suggested that the activated AbaRloC cleaved only a minor fraction of its tRNA targets. As shown, at 4 h of exposure to 30 mg/ml NAL the peak signal of the 3 0 fragment of tRNA Glu amounted only to 40% of the intact tRNA Glu signal ( Figure 7C, lane  8). Moreover, during prolonged incubation the fragment signal sharply declined while that of the intact tRNA hardly changed (lanes 9 and 10). The true proportion of the cleavage product could be yet smaller because it was expected to hybridize to the probe more efficiently than the self-folding intact tRNA. Thus, the NAL-activated indigenous ACNase either partly disrupted its targets population-wide or depleted them in only a minor fraction of the cell population. The dynamics of RNA cleavage products induced by 60 mg/ml NAL were followed essentially as in A in the absence or presence of 100 mg/ml SH. (C) The NAL-activated indigenous AbaRloC does not deplete its tRNA target. RNA extracted from A. baylyi ADP1 at the indicated times of exposure to 30 mg/ml NAL was separated by denaturing gel electrophoresis, stained with EtBr (left panel) and subjected to northern blot analysis using a 5 0 -labelled DNA probe complementary to tRNA Glu residues 35-72 (right panel). Product X and 3 0 tRNA fragments generated by AbaRloC were the respective 60 and 42 nt markers.

Identification of X
Partial cleavage by RNase T1 suggested that X comprises a heterogeneous population sharing 3 0 and possibly 5 0 proximal sequences but highly variable internally (Supplementary Figure S3D). This pattern, an apparently uniform size of $60 nt and a shared 5 0 -OH containing terminal U (Supplementary Figure S3C) suggested that X comprises crRNA molecules (40,48). This assumption was confirmed by splint ligation where weakly labelled X devoid of any 3 0 cyclic-P (40) was probed with the strongly labelled ligation partner juxtaposed by a bridging template matching also the 3 0 repeat portion of A. baylyi ADP1 crRNA ( Figure 8B, Set-up I). Following T4 DNA ligasemediated ligation and phosphatase treatment to remove remaining 5 0 -P groups, Set-up I yielded a ligation product migrating as would X extended by the 13-nt ligation partner ( Figure 8C, lane 1). Thus, X could be derived from a CRISPR transcript by canonical transesterifying cleavages 3 0 to U 8 of the 28-nt repeats. X was not ligated to the probe in the control Set-ups II and III (lanes 2 and 3) where the bridging templates matched 'mock' crRNAs resulting from Cascade-mediated cleavage 1 or 2 positions 3 0 to the canonical site (lanes 2 and 3). The small ligation products seen near the 30-nt marker in lanes 1 and 2 could be derived from self-folded forms of the respective bridging templates.

DISCUSSION
The resemblance of RloC to Rad50 and PrrC and characterized properties of one of RloC's orthologs have portrayed this ACNase as a potential DSB-responsive translation disabler (1,2). The foregoing data concur with this idea: heterologously expressed RloC ACNases were activated by DSB inducers to an extent that limited their own expression (Figures 3 and 4), possibly through depletion of targeted tRNA species and consequent translation arrest. A DSB inducer activated RloC's ACNase also in a natural host ( Figure 5). However, in this case, only a minor fraction of the targeted tRNAs was disrupted (Figure 7). Whether this outcome reflected the partial cleavage of these tRNAs throughout the population or their depletion in only a minor fraction of the cells is not known. This uncertainty leads us to consider alternative ways in which RloC could benefit the genotoxicated host. Basically we consider two possibilities that need not be mutually exclusive. First, RloC benefits its host as a persistence factor augmenting tolerance to DSB inducers. Second, RloC acts as an antiviral contingency when type I DNA restriction is alleviated following DNA damage (11) and/or activated by the phage-mediated host DNA degradation. Relevant to the second task seems the possibility that RloC and the co-activated adaptive immunity ( Figures 5, 7 and 8) provide the genotoxicated host with complementary indicates that ligation is expected. In the control Set-ups II and III the bridging oligonucleotide contains one or two extra A residues expected to preclude ligation (indicated by X). (C) Ligation products were separated by denaturing gel electrophoresis and monitored by autoradiography. crRNA (X) and xylene cyanoll were respective 60 and 30 nt markers. CRL-crRNA ligated to the probe, SFL-self-folded bridges ligated to the probe. defences against impending infections. These and other notions invoked by the data are discussed below.

RloC as a DSB sensor/effector
Linear DNA was required not only for activating GkaRloC ACNase but also for sustaining its activity (Figure 1). Given that the linear DNA mimics DSB, the need for its continued presence portrays RloC as a sensor/ effector responsive only during the lifetime of its triggering DSB signal. In other words, once activated, RloC ACNase keeps functioning as long as it remains tethered to DSB termini. Its removal with the repair of the DNA lesion prompts its inactivation, thus precluding unnecessary loss of needed tRNAs. The continued activity of RloC's ACNase seems to depend also on its ability to retain the nucleotide hydrolysed during its activation (Figure 2), a process likely to entail a single cycle of ATP binding and hydrolysis (2). Stabilization of protein on/off states by persistent binding of a hydrolysed nucleotide characterizes key cellular processes, including DNA replication initiation (49) and G protein signalling (50). In the case of RloC, the bound hydrolysed nucleotide may safeguard continued ACNase activity as long as the DSB signal persists but also elicit rapid inactivation of the ACNase when its function is no longer needed, e.g. by dissociating from the protein as the latter departs from the tethered DNA.

RloC's natural tRNA targets
PrrC's natural specificity for tRNA Lys is compromised when ectopically expressed (24). Therefore, it was not possible to tell if the preferred cleavage of tRNA Glu by the ectopically expressed GkaRloC reflects the true substrate specificity of this ACNase (1). However, the identification of tRNA Glu as a target of indigenous AbaRloC ( Figure 6) and poor sequence identity of the two orthologs (33%) suggest tRNA Glu as a common RloC target. AbaRloC's dual tRNA Glu /tRNA Gln specificity suggests their shared s 2 mnm 5 U wobble base and U at the second anticodon position as RloC substrate recognition cues. However, tRNA Lys that also shares these features was not detectably cleaved, pointing at the third anticodon base as a possible added cue. Specifically, C and G that occupy the third position in the targeted species have similar Watson-Crick donor-acceptor edges. This similarity may be exploited by the ACNase to counter-select the different edge of the matching U of tRNA Lys .
Can RloC team with an R-M protein in trans?
As mentioned, a genetically linked type Ic R-M protein silences PrrC's ACNase in the uninfected cell while a phage-encoded anti-DNA restriction factor alleviates this silencing. In contrast, RloC, which is only rarely linked to an R-M protein, has an internal ACNase switch based on its DSB-responsive ATPase. This coincidence underlies a proposal that the advent of the internal ACNase switch rendered unnecessary the reliance of a PrrC-like progenitor of RloC on the external R-M regulator (2). The DSB-induced activation of AbaRloC in the R-M deficient mutants of its natural host ( Figure 5B) concurs with this assumption. Namely, it indicates that the R-M protein neither silences RloC's ACNase nor mediates its DSB-induced activation. Nonetheless, we cannot rule out that AbaRloC and other RloC orthologs not linked to an R-M protein can be still regulated by the latter in different situations such as phage infections. In favour of this possibility argues RloC's ability to excise the wobble nucleotide, as this property could have evolved to defy phage-induced tRNA repair that counteracts the milder lesion inflicted by PrrC. As T4 family members are the only phage known to encode tRNA repair (http://phage.ggc.edu/), RloC may benefit its hosts as a means to counteract this class. We consider therefore the possibility that RloC ACNases may interface an R-M protein even if not genetically linked to it. As such, they may be activated by phage anti-DNA restriction factors. However, RloC's activation during infection by a T4 family phage may be induced rather by the expected massive degradation of the host DNA, as with the T4 prototype (51). Such vigorous DNA disruption could prompt fast activation of the ACNase, sufficient to abort the infection. The hastened activation of AbaRloC's ACNase at the higher level of the DSB inducer or when the synthesis of DNA repair proteins was blocked by SH ( Figure 7B) supports this possibility.

How can RloC benefit the genotoxicated host?
Preliminary toxicity tests indicated that wild-type A. baylyi ADP1 was as sensitive to NAL as its DrloC, DhsdR, DhsdM DhsdS mutants. Although the DrloC mutant was less viable than wild type after prolonged incubation (regardless of NAL), the R-M mutants showed similar or even lower viability under these conditions (Supplementary Figure S4 and data not shown). We assume therefore that presence of the tdk-kan R cassette used to create all the deletion mutants is responsible for their lower viability (45). Thus, the DSB-induced activation of the indigenous ACNase may have provided the genotoxicated host with subtle or even no advantage. The possibility that RloC confers a subtle advantage such as an increasing tolerance to DSB inducers will have to be examined by suitable persistence tests. On the other hand, lack of advantage could reflect RloC's 'imperfection', being accidentally activated throughout the cell population, cleaving only a minor fraction of the tRNA targets in all or the majority of the cells. Such accidental activation could reflect a potential for a more vigorous activation in direr circumstances that call for cytostatic or suicidal cessation of translation, e.g. the phageinduced host DNA degradation already mentioned. RloC seems suited to ward off such phage also owing to the immediate shut-off of host gene expression their infection entails (51). The lack of host gene expression need not affect the activation of pre-existing RloC molecules but could preclude the accumulation of crRNA, which depended in A. baylyi ADP1 on new protein synthesis ( Figure 7B).
Alternatively, the NAL-activated, indigenous AbaRloC could provide advantage to a minor fraction of the cell population in which the targeted tRNA species were depleted, leading to cytostatic or even suicidal translation arrest. A cytostatic response could rescue such outliers from DSB-induced lytic prophage development, while a suicidal response could prevent superinfection of sibling cells, i.e. if the prophage lacked immunity. Relevant to such scenarios may be perhaps the P2-family prophage spanning A. baylyi ADP1 genes 2137-2200, none of which seems to encode a homolog of a known phage immunity protein. As a lytic prophage thwarter, RloC could also make up for shortcomings of the co-activated adaptive immunity. Namely, although an artificial CRISPR construct could prevent lytic phage induction (52), the adaptive immunity may not be designed to do so in nature, judged by the degenerate CRISPR repeats flanking autoimmune spacers (53).
As mentioned, DNA insults cause bacteria to shut-off a key antiviral device, the type I R-M system (9). This measure prevents degradation of fully unmodified DNA synthesized during the recovery from DNA damage but at the cost of increased susceptibility to external infections (10,11). The coincident DSB-induced activations of RloC and the adaptive immunity may therefore signify their ability to provide the genotoxicated host with compensatory complementary protection from impending infections. A case of complementing antiviral defences provided by the adaptive immunity and a type II R-M system has been recently reported (54).


extfieldsofstudy:
Biology
Computer Science
Medicine

text:
Probabilistic alignment leads to improved accuracy and read coverage for bisulfite sequencing data

Background DNA methylation has been linked to many important biological phenomena. Researchers have recently begun to sequence bisulfite treated DNA to determine its pattern of methylation. However, sequencing reads from bisulfite-converted DNA can vary significantly from the reference genome because of incomplete bisulfite conversion, genome variation, sequencing errors, and poor quality bases. Therefore, it is often difficult to align reads to the correct locations in the reference genome. Furthermore, bisulfite sequencing experiments have the additional complexity of having to estimate the DNA methylation levels within the sample. Results Here, we present a highly accurate probabilistic algorithm, which is an extension of the Genomic Next-generation Universal MAPper to accommodate bisulfite sequencing data (GNUMAP-bs), that addresses the computational problems associated with aligning bisulfite sequencing data to a reference genome. GNUMAP-bs integrates uncertainty from read and mapping qualities to help resolve the difference between poor quality bases and the ambiguity inherent in bisulfite conversion. We tested GNUMAP-bs and other commonly-used bisulfite alignment methods using both simulated and real bisulfite reads and found that GNUMAP-bs and other dynamic programming methods were more accurate than the more heuristic methods. Conclusions The GNUMAP-bs aligner is a highly accurate alignment approach for processing the data from bisulfite sequencing experiments. The GNUMAP-bs algorithm is freely available for download at: http://dna.cs.byu.edu/gnumap. The software runs on multiple threads and multiple processors to increase the alignment speed.

Background
DNA methylation is a chemical process in which a methyl group is added to the carbon-5 position of a DNA cytosine. In most vertebrates, DNA methylation typically occurs on the cytosine of a CpG dinucleotide [1,2], although some specific examples of other types of methylation have been shown to play roles in specific tissues [3][4][5][6]. Since its discovery over 60 years ago, DNA methylation has been linked to many important biological phenomena such as the suppression of gene expression http://www.biomedcentral.com/1471-2105/14/337 location for a bisulfite-treated read (BSR) is a complicated and difficult process. Although cytosine to thymine changes are allowed for when mapping the BSRs to the genomic sequence, methylated and unmethylated CpG locations are often identified and, in some cases, it is impossible to distinguish between a bisulfite (BS)-treated thymine that originated from an unmethylated cytosine, and a true thymine from a different genomic location or an individual genomic variation at that location [13,14]. This ambiguity is often magnified by the presence of sequencing errors or low-quality bases. As a result, sophisticated computational strategies are required for aligning reads from BS-seq experiments.
The first whole-genome methylation profiles were performed on Arabidopsis thaliana. To map the resultant BSRs, alignment algorithms based on a probabilistic formulation and a suffix tree [15] as well as reference genome conversion [16] were used. However, it is not computationally feasible to apply this approach to larger genomes such as the human genome, or to experiments with the current deeper sequencing depths. Later algorithms, such as BSMAP [14], constructed seed tables of locations from both the original reference sequence and the BS variants, and then extends the seeds to form a possible mapping location. This seed extension process can be somewhat unreliable because the seeds must be exact and do not take into account BS-treated variations. BS alignment methods such as BS Seeker [17], Bismark [18], and BRAT-BW [19] have been used to map BSRs. These methods employ a Burrows-Wheeler transformation [20] for fast (in)exact string matching and then combine the results with either a pre-processing or post-processing script to handle BSRs with three letters after converting all Cs to T. The strength of these methods lie with reads with fewer mismatches, but they offer very limited support for aligning reads with insertions or deletions (indels). These methods also have difficulty in aligning Ts in the reads to Cs in the genome without also (incorrectly) aligning for cysteines in the reads to thymidines in the genome. More recently, the LAST alignment program has been adapted to align BSRs [21]. LAST uses a seed extension approach similar to the one used by NCBI BLAST [22], but the speed and mapping accuracy are increased by using variable-length seeds and base quality information [23].
In this study, we present a highly accurate probabilistic mapping algorithm, Genomic Next-generation Universal MAPper for Bisulfite Sequencing (GNUMAPbs), for BSR alignment. GNUMAP-bs is an extension of the Genomic Next-Generation Universal MAPper (GNUMAP) [24] that can accommodate the alignment of BSRs to a reference genome. GNUMAP-bs was developed to achieve higher accuracy than other BS-seq approaches by including base quality scores in the alignment pro-cess. In addition, the GNUMAP-bs probabilistic mapping approach allows for the unbiased estimation of DNA methylation, especially when reads are aligned to multiple genomic locations.

Results and discussion
We compared GNUMAP-bs with several commonly-used BS alignment methods, namely BSMAP (v2.3) [14], Bismark(v0.7.9) [18], LAST [21], and BRAT-BW [19]. We also compared GNUMAP-bs with the unpublished proprietary probabilistic aligner Novoalign (v2.07.09, http://www. novocraft.com) because it uses a probabilistic approach that is similar to the approach used in GNUMAP-bs. For a fair comparison, we considered only common options that were shared by all the alignment methods; otherwise, we used the default values or slight modifications of the defaults to achieve the best performance for each method. For each aligner, we allowed up to three or four base pair differences within 100-bp BSRs. For GNUMAPbs and Novoalign, this resulted in cutoff values of -a 0.90 and -t 75, respectively. Considering that repeated genomic sequences often appear in mammalian genomes, we allowed up to 20 multiple mapping locations for each read. This approach to handling multiple alignments has only a small impact on the sensitivity of the more heuristic algorithms, but has moderate sensitivity gains in more traditional approaches like GNUMAP-bs. The Bismark software only reports up to two valid alignments. The parameter values that were used for each alignment method are given in Table 1.
We used two datasets to evaluate the performance of each of the alignment methods. The first was a simulated BSR dataset, which was carefully designed to mimic a typical human methylation experiment. The second was a real human BS-seq dataset, which was compared with an experimentally-derived human methylome profile. LAST -Q1 -j1 -d120 -n20 -f1 | last-map-probs.py -s150 -m0.95

BRAT-BW default settings for single-end reads
Parameters in parentheses were used for the human BSR dataset. Up to 3 to 4 sequencing errors or mutations were allowed in a 100-bp long BSR. Because the human genome contains many repeated genomic regions, 20 valid mapping locations were allowed for a given read when the alignment method supported this feature. http://www.biomedcentral.com/1471-2105/14/337

Simulated bisulfite sequencing experiment
For the simulation study, we randomly assigned 20% of the CGs in the whole human genome (NCBI build37/HG19) to represent unmethylated cytosines (Cs) by changing them to thymidines (Ts), thereby simulating complete BS conversion. For the remaining 80% of the CGs, 75% was randomly assigned to be fully methylated and, therefore the Cs remained as Cs. The remaining 5% was assigned to be methylated in proportions between 0.1 and 0.9. In this dataset, we assumed that all non-CG sites remained unmethylated, so these Cs were all changed to Ts for read generation. We used the dwgsim (http://sourceforge.net/ projects/dnaa/files/dwgsim) simulation tool with parameters -e 0.001-0.008 -1 100 -y 0.05 -r 0.002 -R 0.2 -C 10, to generate 100-bp BSRs with a 10? read depth across the genome. This simulation produced a BS-seq dataset that contained approximately 180 million (M) reads with a sequencing error rate that ranged from 0.001 to 0.008 and increased from the 5' to 3' ends, plus 5% randomly generated sequence, and a mutation rate of 0.002 (in which 0.2% were indels). We found that the mapping sensitivities of the probabilistic aligners GNUMAP-bs (97.0%), Novoalign (96.3%), and LAST (96.9%), were higher than the sensitivities of BSMAP (93.9%) and Bismark (93.2%), as shown in Table 2. The Bismark-Bowtie2 algorithm had the lowest false positive rate (0.1%), but at the cost of an approximately 7% decrease in its sensitivity. Overall, the Bismark-Bowtie2 algorithm had 1.2 M fewer false positives than GNUMAPbs, but approximately 10 M fewer BSRs were aligned correctly to the genome. The Bismark-Bowtie algorithm performed better than Bismark-Bowtie2 in that it aligned 5 M more reads with only a slight increase in the number of incorrectly aligned reads. The LAST algorithm aligned nearly all the BSRs (99.2%) to the genome, but although the proportion of correctly mapped reads was similar to that of GNUMAP-bs, the error rate was much higher. BRAT-BW was not included in this comparison because the application software removed the original sequence header/name which contained information pointing to the correct alignment location of the BSRs.
Differences in sensitivities of these methods were much more pronounced for BSRs that contained at least one sequencing error or genome variant, and the GNUMAPbs sensitivity was clearly better than the sensitivities of the other approaches (Table 2). We also evaluated how accurately each aligner predicted the CG methylation levels when the true methylation level ranged from 10% to 90%. All the alignment methods tested performed well with mean absolute errors less than 0.01 and with standard deviations less than 0.07.
When the computational performances of each of the methods were compared (Table 2), we observed that GNUMAP-bs required the most RAM (44.8 GB) and Bismark required the least (5.9 GB). Because some of the software applications presented here support Simulation study of 160 million (M) simulated BSRs generated from the human genome reference sequence. The GNUMAP-bs algorithm was the most sensitive aligner, especially for reads with ≥ 1 sequence variant (sequencing errors or mutations). The Bismark algorithm had the smallest error rate with 1.2 M fewer erroneously assigned reads than GNUMAP-bs, however GNUMAP-bs correctly aligned 6 to 10 M more reads. The BSMAP algorithm had the fastest total run time, however its sensitivity was less than the sensitivity of the GNUMAP-bs algorithm. LAST mapped nearly all reads with a sensitivity that was comparable to that of GNUMAP-bs, but the mapping error rate for LAST was much higher than it was for GNUMAP-bs. http://www.biomedcentral.com/1471-2105/14/337 computation on multiple threads and some do not, we presented two different measures of computational speed: 1) the total run time on a 16 CPU linux server and 2) the number of reads processed per CPU per second. GNUMAP-bs required approximately 40 hours of total compute time to process the 180 M BSRs, while BSMAP, the fastest algorithm, was nearly nine times faster than GNUMAP-bs. The LAST application does not support parallel computing, so LAST had the longest total alignment time. However, the LAST algorithm aligned the most reads per second per CPU (753), while Bismark-Bowtie2 aligned the least (26).

Human BS-seq dataset
We also evaluated the performances of the alignment methods using a BS-seq dataset generated from samples from a healthy human donor collected at the Andrology Laboratory at the University of Utah (Salt Lake City, UT, USA). The BS-seq data were generated by coupling BS conversion and the Illumina HiSeq2000 platform, which generated 101-bp BSRs for analysis. We aligned 283.6 M reads from three lanes of BS sequencing data containing 85 M to 100 M sequencing reads each, to the recent build of the human genome (NCBI build37/HG19). The human BSRs were processed for quality control as suggested previously [25]. Briefly, the quality control involved masking low quality bases or trimming consecutive lowest quality bases at the 3' ends of the reads. We used the same parameters on these data as were used for the simulation experiment (Table 1)  We compared these mapping results with the Sangerbased BS sequencing control data available from the Human Epigenome Project (HEP) [26] by selecting the data that were obtained using the same type of tissue as was used in our BS-seq dataset. The HEP data provides a natural gold standard for algorithmic evaluation. For this comparison, we focused on the 13,563 HEP CG sites on chromosome 22 (chr22), because these data provide the most comprehensive chromosomal CG coverage that is available in the HEP database. The data for the other chromosomes in the HEP showed similar profiles as chr22, but the coverage was much sparser.
Overall, the mapping results produced by the different alignment methods for the human BS-seq dataset were consistent with the results for the simulation dataset. Compared with GNUMAP-bs, LAST aligned almost twice as many reads to chr22 (3.02 M compare with 1.57 M); however, the CG read coverage for LAST (7.6 reads/CG) was only about 20% higher than for GNUMAP-bs (6.3 reads/CG) ( Table 3). While Novoalign aligned more reads to and covered more CG sites on chr22 than GNUMAP-bs, GNUMAP-bs had better CG read coverage (6.3 reads/CG) than Novoalign (6.1 reads/CG) and the other aligners, BSMAP (5.8 reads/CG), Bismark (5.8 reads/CG), Bismark-Bowtie2 (5.7 reads/CG), and BRAT-BW (5.0 reads/CG) ( Table 3).
To determine if increased coverage resulted in decreased alignment quality, we computed the correlation coefficient between the estimated HEP methylation levels and the methylation levels estimated from each alignment method. The correlation coefficients were 0.887 for GNUMAP-bs and ranged from 0.895 for Bismark-Bowtie2 to 0.882 for LAST (Table 3). We also computed a concordance statistic for each of the methods as defined previously [27]. Briefly, the concordance is the fraction of sites for which the methylation levels (aligner vs. HEP) differ by less than a predefined cutoff (we used 0.25). Based on this statistic, we found that the performances of most of the aligners were similar, with concordance values of 0.869 for GNUMAP-bs and ranged from 0.872 for Bismark-Bowtie2 0.872 and BRAT-BW to 0.858 for LAST (Table 3). These results showed that the increased coverage levels produced by LAST result in lower consistency with the HEP data in terms of both correlation and concordance. Pairwise comparisons of the concordance for the CGs covered/not covered between GNUMAP-bs and the other aligners are shown in Figure 1. For example, Figure 1 shows 577 CGs that were covered by GNUMAPbs but not covered by Bismark; the concordance of the GNUMAP-bs methylation estimates for these sites was 0.78. In contrast, only six CGs covered by Bismark were not covered by GNUMAP-bs, and the concordance for these sites was extremely low (0.33). These results clearly showed that the increased mapping coverage (for both numbers of CGs covered and reads per CG) obtained using GNUMAP-bs did not result in decreased alignment quality or methylation estimates.

Conclusions
BS sequencing presents difficult challenges to researchers attempting to process the sequencing reads from BSseq experiments. In this work, we present GNUMAPbs, a highly accurate and effective alignment algorithm http://www.biomedcentral.com/1471-2105/14/337 that is specifically designed to estimate DNA methylation levels with base-level resolution in BS-seq data. GNUMAP-bs uses a probabilistic approach to align BSRs to a reference genome. GNUMAP-bs was developed to achieve higher coverage and accuracy than other published BS-seq approaches by integrating base quality and alignment quality information in the mapping process. We have shown that the GNUMAP-bs probabilistic mapping approach results in an improved unbiased estimation of DNA methylation across the human genome. In simulated and real datasets, we showed that GNUMAPbs outperforms other BS-seq alignment methods when both coverage and consistency were balanced with Sanger based BS sequencing controls.
In addition, GNUMAP-bs provides many high-demand features needed for constructing a high quality methylome from BS-seq data. First, GNUMAP-bs incorporates quality sequencing data into a dynamic programming framework. This feature gives GNUMAP-bs the best balance between sensitivity and specificity of the tested BSR aligners, especially for reads that contain short polymorphisms. Second, GNUMAP-bs adaptively assigns an optimal mapping stringency based on an effective read length after the original read is trimmed. Third, GNUMAP-bs where the corresponding original bisulfite read (BSR) is known (e.g. CGGTGTC), a query k-mer (e.g. TGGTG) is searched against the BS-converted hash table. Locations with a k-mer match to both the genome and the read are aligned using the original read and original genome and the GNUMAP-bs probabilistic Needleman-Wunsch (N-W) algorithm. If the alignment score passes a quality threshold, the location is considered a match and recorded on the genome for future output. The posterior probability is computed for each mapped location using Equation (1). Finally, all the mapped BSRs are used to quantify the fraction of methylation at each CG location across the reference genome. http://www.biomedcentral.com/1471-2105/14/337 not only relies on the maximal score alignment but also probabilistically considers suboptimal alignments; that is, the alignment score is converted to a posterior probability and the probabilistic scores quantify the likelihood of the true source location for each read across the reference genome. As a result, for both the simulated and real datasets, we showed that GNUMAP-bs was more effective that the other methods in detecting read locations in the presence of sequencing errors. GNUMAP-bs also displayed the highest consistency with a known HEP methylation database. Because GNUMAP-bs supports Message Passing Interface (MPI) processing, the computational burden of the dynamic programming can be alleviated. Moreover, with computer resources becoming cheaper, computing clusters with large numbers of nodes and cores, and more computing clouds are becoming available. Therefore, memory and CPU running times are less of a bottleneck, which is especially useful for GNUMAPbs alignments. For this reason, accuracy should currently be a more important concern in BS-seq data analysis.

Methods
The GNUMAP-bs alignment algorithm is a modification of the GNUMAP algorithm, which consists of three main steps, all of which needed to be modified to align BSRs to a reference genome. A flow chart of the GNUMAPbs algorithm is displayed in Figure 2. The first step is the construction of a hash table using all genomic subsequences, where the k nucleotide (nt) long (k-mers) are the keys and the hash table values store the genomic locations of the k-mer. In addition, k-mers from the BSRs are incrementally referenced in the reads in the genomic hash table. In GNUMAP-bs, the genome and the reads are artificially 'BS-converted' by changing all Cs to Ts before the hashing step. This process ensures that all BSRs can be referenced into the hash table regardless of whether or not they contain methylated bases.
In the second step of the GNUMAP algorithm, all the BSRs are aligned to the genome at the hashed locations. The alignment is performed using a novel probabilistic alignment algorithm that uses base quality information. All the regions that meet a predefined alignment score threshold are retained for the final step. In GNUMAP-bs, the probabilistic alignment algorithm matches the unconverted (original) reads to the unconverted genome while allowing matches between Ts in the reads and Cs in the genome. Using dynamic programming in probability space for mapping BSRs to a genome has several benefits: 1) the Needleman-Wunsch algorithm is guaranteed to find the optimal alignment for a BSR; 2) by incorporating base qualities into the probabilistic algorithm, true DNA methylation can be more accurately identified in the alignment, especially in areas where the reads have low base quality; and 3) by making only a small change in the alignment scoring matrix, namely by removing the T (read) to C (genome) 'mismatch' penalty and scoring these alignments as a 'match' , the probabilistic Needleman-Wunsch algorithm can accurately account for the BS changes in BSRs without losing reads that are partially converted, or that have genome variations or sequencing errors.
During PCR amplification after BS conversion, the C to T conversions in the BSRs lead to G to A conversions on the PCR-synthesized strand as described previously [14]. These G to A changes provide evidence of nonmethylation on the original strand. The phenomenon also occurs for C to T changes on the PCR-synthesized strand. Therefore, the GNUMAP-bs algorithm maps each BSR twice, once looking for only C to T changes and then looking for only G to A changes. When a G to A change is observed, the read is considered to be a PCR-synthesized strand, and the non-methylation event is attributed to the original strand.
In the final step, the GNUMAP-bs algorithm quantifies the significance of an alignment (e.g. a posterior probability) by assigning a mapping score that is proportional to the relative quality of alignments and inversely proportional to the number of matching sites in the reference genome. Given a BSR, r, assume that there are J "best match" locations (potential mapping locations) for which the alignment score is q[d] for alignment to genomic position d, which is also aligned to the first base pair location of r, matches some user-defined threshold. Then, the posterior probability at each location P(r d ) can be computed as: It is worth noting that the calibrated posterior probability is a special case of the LAST [21] mapping probabilities equation where a scaling factor for a bit score in the simple E-value statistics is not explicitly computed.
This probability provides an intuitive value between 0 and 1 representing the methylation signal for each read at each CG site. Other BS mapping programs use 1 #mappings instead of Equation (1), which may not accurately capture the true posterior probability. The posterior probability indicates the relative significance of mapping a read to a particular position. The probabilistic Needleman-Wunsch computes a log-likelihood score; therefore, if a particular BSR has multiple possible mapping locations in the reference genome, the significance of each mapping decreases. After GNUMAP-bs computes P(r d ) for each read, the algorithm combines the probabilities into one methylation profile vector and infers a true methylation ratio at each CG site.
For example, consider a nucleotide, n i at genomic location i, and assume that n i is a C residue in a reference genome. Let C (m) be a methylated C residue and L(r) be the read length. For a set of aligned reads covering i, the methylation ratio c[i] can be computed as: Intuitively, the reads that cover each C location across the reference genome provide evidence of true map-ping so that a reliable methylation percentage can be obtained.

Software implementation and availability
The GNUMAP-bs algorithm is integrated into the is GNUMAP software suite, and is freely available for download at http://dna.cs.byu.edu/gnumap. The GNUMAP-bs pipeline can accommodate single-end or pair-end reads in either FASTA or FASTQ file formats. A reference genome file (or multiple reference files) in FASTA format is also required. However, to increase the efficiency of the workflow, users can opt to write the reference genome hash table to a file, which can be used in future runs. GNUMAP outputs read alignments in standard SAM file format. The GNUMAP-bs software suite contains an addition application function (sam2gmp) that summarizes a SAM file and writes it into a text file that contains one line for each cytosine in the genome. Each line in this file contains the chromosome number, the location of the cytosine on the chromosome, the number of reads, and the numbers of As, Cs, Gs, Ts, and Ns covering the location. In addition, the text file gives a likelihood ratio p-value that indicates whether there is a significant number of Cs at the location (i.e. methylation significance).
In addition to the adaptations for BSRs, the GNUMAPbs software also contains several modifications to reduce the computational time and memory needed for the alignments. For example, the initial read and genome conversion step reduces the genome alphabet to three bases, leading to an increase in the number of genome locations identified in the hash table. To reduce this effect, multiple seeds from each read are referenced by GNUMAP-bs, and the alignment is only conducted on the locations with the top two most k-mer hash references. Furthermore, although the original GNUMAP software supports multi-threaded computing within the same node using pthread, GNUMAP-bs is fully enabled to support MPI processing. This feature allows a large-scale alignment to be spread across multiple nodes in a cluster or supercomputing facility [28]. In this implementation, the genome is divided into equal parts across nodes, and then the same batches of reads are aligned by each node to their individual portion of the genome. Once the batch is completed, the nodes communicate via MPI to calculate the posterior probability scores. Because most of the CPU time is spent on the alignments, the communication overhead is relatively small, resulting in a highly efficient parallel algorithm. http://www.biomedcentral.com/1471-2105/14/337 and BRC generated the human BS-seq datasets. CH and WEJ wrote the paper. All authors read and approved the final manuscript.


extfieldsofstudy: Biology Medicine

s2fieldsofstudy: Biology Medicine

text:
Cold Atmospheric Plasma for Selectively Ablating Metastatic Breast Cancer Cells

Traditional breast cancer treatments such as surgery and radiotherapy contain many inherent limitations with regards to incomplete and nonselective tumor ablation. Cold atomospheric plasma (CAP) is an ionized gas where the ion temperature is close to room temperature. It contains electrons, charged particles, radicals, various excited molecules, UV photons and transient electric fields. These various compositional elements have the potential to either enhance and promote cellular activity, or disrupt and destroy them. In particular, based on this unique composition, CAP could offer a minimally-invasive surgical approach allowing for specific cancer cell or tumor tissue removal without influencing healthy cells. Thus, the objective of this research is to investigate a novel CAP-based therapy for selectively bone metastatic breast cancer treatment. For this purpose, human metastatic breast cancer (BrCa) cells and bone marrow derived human mesenchymal stem cells (MSCs) were separately treated with CAP, and behavioral changes were evaluated after 1, 3, and 5 days of culture. With different treatment times, different BrCa and MSC cell responses were observed. Our results showed that BrCa cells were more sensitive to these CAP treatments than MSCs under plasma dose conditions tested. It demonstrated that CAP can selectively ablate metastatic BrCa cells in vitro without damaging healthy MSCs at the metastatic bone site. In addition, our study showed that CAP treatment can significantly inhibit the migration and invasion of BrCa cells. The results suggest the great potential of CAP for breast cancer therapy.

Introduction
Breast cancer is the second leading cause of cancer deaths in women. It is estimated that 232,340 new cases of invasive breast cancer will be diagnosed in the United States and 39,620 women will die of the disease in 2013, according to the American Cancer Society. Breast cancer exhibits an affinity to metastasize to bone, resulting in debilitating skeletal complications associated with significant morbidity and poor prognosis. Roughly, 85% of individuals eventually develop bone metastases in advanced breast cancer [1]. The growth of disseminated tumor cells in the skeleton requires tumor cells to inhabit the bone marrow wherein metastatic breast cancer cells colonize the skeleton and interrupt normal bone remodeling processes. Current breast cancer treatment options such as surgery and radiotherapy contain severe limitations with regards to nonselective and incomplete tumor ablation. Thus, new treatments which can completely and selectively ablate solid tumors and transient breast cancer cells and tissues while keeping surrounding healthy cells and tissues intact is highly desirable.
Plasma is an ionized medium that contains numerous active components including electrons and ions, free radicals, reactive molecules, and photons [2][3][4]. It can be categorized as either thermal or non-thermal plasma [4]. Thermal plasma has been widely used to modify material surfaces, which is generally conducted in a vacuum [5,6]. Cold atmospheric plasma (CAP), a non-thermal plasma, shows high electron temperatures but very low gas temperatures due to a weak ionization rate [5]. Because of the varied mass of CAP particles and constituents, thermodynamic equilibrium of electron self-collision occurs much faster than equilibrium between electrons and larger particles, such as ions [6,7]. Thus, the overall plasma temperature is much lower than the electron temperature, which is close to room temperature. As an emerging technique for biomedical applications, CAP exposure has been shown to be highly effective in germicidal and sterilization, wound healing, blood coagulation, material surface modifications and crosslinking and treatment of various diseases, including cancer [8][9][10][11][12][13][14]. In particular, CAP can potentially offer a minimally-invasive surgical approach allowing for specific cancer cell or tumor tissue removal without influencing surrounding healthy cells and tissues, thus making it a promising technology for cancer therapy. Two possible underlying mechanisms for CAP's high selectivity towards cancer cells can be attributed to: (1) complex CAP composition and (2) the different nature of cancer cells and normal cells. Firstly, the variable cellular effects of CAP may be explained by the complex chemical composition of the CAP plume governed by the way the CAP is generated. Its presence can promote specific chemical reactions between charged particles and living cells triggering intracellular biochemical reactions that would elicit desired breast cancer therapeutic effects [15]. For example, reactive oxygen species (ROS, such as O, OH) are believed to be a major reason for cancer cell lysis [16] and bacterial inactivation [17], although the contribution of each type of the reactive species is still not well understood. In addition, the selective effects of CAP may be attributed to the significant difference of cancer cells and normal cells, thus rendering cancer cells more susceptible to CAP [18]. The objective of this research is to investigate a CAP-based therapy for the first time to selectively ablate metastatic breast cancer cell in metastatic bone site.
Traditionally, CAP has been developed to kill cells by utilizing deposition of different ROS active particles and excited radical, including UV photons [19]. Low dose of cold plasma treatment may have little to no effect on cancer cells, while high doses may ablate both cancer cells/tumor tissue and healthy cells/tissues. In this study, we will optimize CAP treatment parameters to selectively kill human metastatic breast cancer (BrCa) cells without minimally influencing healthy human bone marrow mesenchymal stem cells (MSCs).
For this purpose, three key CAP factors (i.e. treatment condition, exposure time, and treatment frequency) were tuned and evaluated. Firstly, we evaluated a series of CAP treatment conditions' (such as generator's output voltage and resistance) influence on MSC viability in vitro. After selecting the best treatment condition (i.e., minimum negative influence on MSCs), human metastatic MDA-MB-231 BrCa cells and MSCs with the same cell density were separately treated for 0, 30, 60, 90, and 120 s under this condition. Furthermore, BrCa cells and MSCs were cultured for 1, 3 and 5 days with daily CAP treatments (0, 30, 60 and 90 s). Cell growth and migration under designed CAP treatment were compared and evaluated.
Primary MSCs were derived from healthy consenting donor's iliac crest (Female; Age: 27). They were obtained at Tulane University under an IRB-approved protocol with informed consent, distributed from the Texas A&M Health Science Center, Institute for Regenerative Medicine and thoroughly characterized. We had the fully executed Material Transfer Agreement that was needed to obtain the cells. This provides verification that we in fact did not obtain the cells from the donors themselves. MSCs (passage #3-6) were cultured in a complete media comprised of Alpha Minimum Essential Medium (a-MEM, Gibco) supplemented with 16.5% fetal bovine serum, 1% (v/v) L-Glutamine (Invitrogen), and 1% penicillin streptomycin solution and cultured under standard cell culture conditions as described above.

CAP characterization
The CAP device can produce diverse chemically active species. In our study, the CAP device consists of four blocks ( Figure 1). Block 1 is a DC power supply 60 V/6A (SPS-606 Insteak). Block 2 is configured with a centrally powered electrode 1 mm in diameter and a ground outer electrode wrapped around a 4.5 mm diameter quartz tube which is the part of cold plasma production. Block 2 mainly controls the positional three degrees of freedom, allowing for rotational, vertical, and horizontal translation of the jet. Block 3 is made up of a capacitor, (VISHAY SPRAGUE-TVA1346-E3-100 mF, 100 V, axial) a transistor (STGW38IH130D, 33 A-1300 V-very fast IGBT), and a LM555 timer. LM555 timer contains two capacitors and two adjustable potentiometers R A and R B . Block 4 is the helium gas supply. The helium flow rate was set to approximately 4.6 L/min. Electrical measurements were conducted with a Tektronix TDS3014C Digital Phosphor  Oscilloscope. Optical measurements were made with a 100 mm probe 3 cm away from the front of the nozzle (Figure 2a). In order to measure the radius of the plasma jet, a holder with a horizontal trench (Figures 2a and 2b) was used to move the probe horizontally. Spectra of the jet were recorded at 0, 1, and 2 mm from the jet center. In addition, a radical experiment was conducted to measure the radius of the plasma jet.  MSC and BrCa cells were seeded at a density of 50,000 cells/ well on two 48 well-plates, respectively. 300 ml of media was added into each well and all of the cells were cultured for 24 hours before CAP treatment. The distance between the CAP nozzle and plate bottom was 3 cm during all treatment conditions. All treatments were conducted at about 25uC (room temperature). After 30 s CAP treatment, media were replaced and cells were cultured for additional 24 hours. Firstly, MSCs were treated with different CAP conditions, including output voltage and resistance value as listed in Table 1. These parameters can alter the CAP discharge power and affect cell behaviors. The attached cell numbers after 24 hours of culture were compared to untreated controls.
Two conditions (4 and 5 in Table 1) without negative inhibited influence to MSCs were prescreened, selected and applied to BrCa cells for 0, 30, 60, and 90 s. After the CAP treatments under the selected conditions, separately treated MSCs and BrCa cells were cultured for an additional 24 hours under standard cell culture conditions. After required culture time, the MTS assay reagent CellTiter96VR Aqueous One solution (Promega) was used to quantify attached cell number. In addition, a live/dead viability/ cytotoxicity kit for mammalian cells (Molecular Probes) was used to evaluate cell viability after CAP treatment (condition 5 for 0, 30, 60, 90 and 120 s). Briefly, after 24 hours, each well was rinsed twice with phosphate buffer saline (PBS). A working solution of 4 mM ethidium homodimer (EthD-1) and 2 mM calcein AM was prepared as described previously [20]. After rinsing, 80 mL of working dye solution was added to each well and incubated at 37uC 30 min before imaging each well with a fluorescent microscope (Olympus BX 60). A mosaic image was collected in Image Pro H to show and compare the distribution of live and dead BrCa cells and MSCs under 0, 30, 60, 90 and 120 s CAP treatments.

BrCa and MSC proliferation under CAP treatments
BrCa and MSC were pre-seeded at 50,000 cells/well with 300 mL of media per well in a 48 well-plate and cultured for 24 hours. Similar to the previous CAP treatment section, CAP was used to treat both of the BrCa cells and MSCs. The optimal conditions for selectively ablating BrCa cells with minimally influence of MSCs were applied daily for 0, 30, 60 and 90 s. The cells were cultured in respective media for 1, 3, and 5 days.
At predetermined time points, MTS assay was used to quantify cell number. Briefly, cells were lift with Typsin EDTA, and added to a new 96 well plate with 100 mL of suspension per well. 20 mL of the MTS reagent was added to each well and incubated for 1 hour. Absorbance was measured at 490 nm using a UV-vis spectrophotometer.

BrCa cell migration study after CAP treatment
In order to evaluate metastatic BrCa cell migration behavior after CAP treatment, the Matrigel Transwell Migration assay using a transwell chamber in conjunction with 8 mm membrane filter inserts was chosen. The pre-seeded BrCa cells were treated by CAP daily with different predetermined treatment times (0, 30, 60 and 90 s, condition 5). Then CAP treated cell were lifted and concentrated via centrifugation and further suspended in serumfree DMEM medium. 50,000 cells were added to the upper transwell chamber, whereas the lower chamber was filled with complete medium. After 24 hours of incubation in 5% CO 2 at 37uC, BrCa cells which had migrated through the membrane to the lower surface were fixed and stained with a Differential Quik Stain Kit. The numbers of cells on the lower surface of the membrane were counted by a light microscope (10x).

2D wound healing assay of BrCa cells after CAP treatment
The wound healing assay was used to study BrCa cell migration velocity and cell interactions after 0, 30, 60 and 90 s CAP treatment. Treated BrCa cells (500,000 per dish) were seeded in a petri dish and allowed to adhere for 24 hours. Confluent monolayer cells were straightly scratched by a pipette tip, to simulate a wound environment, and then washed three times with PBS to remove cell debris and suspended cell. Fresh medium was added, and the cells were allowed to close the ''wound'' for 24 hours. Images of cells migrating into the wound space were captured at 0 to 24 hours via a confocal microscope. The rate of wound closure, BrCa cell migration velocity and migration distance of the first line of bottom of scratch in first 9 hours were calculated.

Statistics
The cellular experiments were run in triplicate and repeated three times for each group. Data are presented as the mean value 6 standard error mean (SEM) and were analyzed with student's t-  test for pair-wise comparison. Statistical significance was considered at p,0.05. Figure 3 shows a typical spectrum of helium CAP jet interacting with the ambient air. The identification of emission lines and bands were mainly according to Volotskova et al [21] and Georgescu et al [22]. In the significant range of 250-300 nm there are very weak emission lines, which are interpreted as NO [23]. In addition, the inactivation effect on bacteria by ultraviolet radiation is mostly related to biological damage caused by UV wavelengths of 200-280 nm. Thus it can be concluded that UV photons are not a major CAP species with our setup. In the wavelength range of 300-550 nm, their magnitudes are at most a few thousandths of the highest peak, N 2 + (391 nm). The species lines between 300-500 nm are still undetermined [2,4]. Other apparent CAP species were interpreted as He lines between 550-750 nm. Atomic oxygen (O) (emission lines at 777 nm) is believed to have significant biological effects, and was present in our generated CAP. The CAP jet is a highly complex environment which combines the comprehensive effects of different ions and neutrals. From the emission spectra data shown here, O (777 nm), OH/O 3 (309 nm), N 2 + (391 nm), N 2 or NO lines (316, 337, 358, 427.5 nm) are the domain species of the spectra. Repeated measurements (not shown) were performed to confirm the initial results. Radical measurement of the CAP jet was conducted with the holder shown in Figure 2. Helium flow rate was kept at 4.6 L/min. Input and output voltage was maintained at 11.8 V, and 3.76 kV, respectively. With a spectra of different positions recorded, the emission drops of major CAP species were plotted in Figure 2 (c). We can see that the radius of the jet is roughly 2 mm (i.e., ,4 mm diameter).

CAP characterization
In addition, it was observed that the plasma between the two electrodes formed by optical emission is much stronger than the plasma plume. Therefore, the spectra measured in our experiment may not be the entirely emission of the plasma jet. Furthermore, it should be pointed out that the probe placed in front of the jet may cause the perturbation of the helium plasma jet.

Different BrCa cell and MSC responses under varied CAP conditions
Given its potential to interact with tissue or cells without a significant temperature increase, CAP holds great potential for minimally invasive cancer treatment. As mentioned previously, possible underlying mechanisms with regards to CAP's high selectivity towards cancer cells may be attributed to complex CAP compositions. In this study, we varied two CAP jet parameters (i.e., voltage and resistance) in order to achieve a desirable CAP composition and cellular responses. Table 1 shows that different CAP voltage and resistance conditions lead to either a positive or negative effect on healthy MSCs. The 'positive effect' means non-inhibited or increased cell growth after CAP treatment. In other words, MSC can keep and improve their proliferation after the treatment. 'Negative' means a decreased cell growth number and changed morphology after CAP treatment. Among all conditions, conditions 4 and 5 reveal desirable positive effects on MSCs. Moreover, Figure 4 shows the selective effect of CAP treatment on MSC and BrCa cells under the condition 5 (RA = 3000 V, RB = 3000 V, and 3.68 KV output voltage). Within 30 s treating time, the MSC growth was greatly improved while BrCa cell growth was inhibited. In contrast, Figure 5 showed the effects of condition 4 on both MSCs and BrCa cells. Although MSC growth was greatly enhanced initially, there is not significant selectivity shown here. After 60 s treatment, both MSC and BrCa cell densities were decreased. Thus, the condition of RA = 3000 V, RB = 3000 V and 3.68 KV output voltage is the selected optimal treatment condition for the later experiments. Figure 6 illustrates the viability of metastatic BrCa and MSCs under 0, 30, 60, 90 and 120 s of selected CAP treatment condition 5. The images show typical and representative areas of the treated cells. Our results show that at short durations of CAP treatment (30 s), MSCs grew well when compared to both controls and the same treated BrCa cells. After 60 s of CAP treatment, significant differences between BrCa and MSC responses were observed. Greater amounts of metastatic BrCa cells were detached and lysed by CAP while MSCs kept growing normally with fewer dead cells ( Figures 6C-E and 6H-J). Moreover, our results demonstrated that, for the first time, CAP treatment may selectively kill BrCa cells while leaving MSC growth unaffected. Although the underlying mechanism is not completely clear, we believe that a unique ROS (plasma killing species) and reactive nitrogen species (RNS, plasma healing species, such as NO, N 2 , N 2 + ) environment   created by CAP may contribute to the adverse and beneficial cell responses as shown in our study. Cell adhesion and viability data after different CAP treatments showed that the CAP device's voltage, resistance, and treatment duration played a significant role in modulating different cellular responses. Furthermore, in addition to single CAP treatments, MSCs and BrCa cells were periodically exposed to daily CAP treatment for up to 5 days. Both MSC and BrCa cell proliferation during the daily CAP treatments were evaluated and compared. The results of 1, 3, and 5 day MSC proliferation studies are shown in Figure 7. Under short (30 s and 60 s) daily CAP treatments, MSCs remained intact and proliferated well when compared to untreated controls (0s treatment). Under 90 s daily treatment, MSC cell density was lessened when compared to short term CAP treatment. In contrast, Figure 8 reveals a significantly different BrCa proliferation profile under the same daily CAP treatment regimen. Cancer cells were more sensitive to CAP when compared to untreated BrCa cells (0 s), all CAP-treated BrCa cells showed a significant inhibited proliferation rate after 3 and 5 days. Especially after 60 s and 90 s daily CAP treatment, cancer cells virtually stopped expanding. According to the results, selective ablation time for significantly lysing BrCa with the smallest possible effect on MSC is between 30-60 s.

Inhibited BrCa cell migration after CAP treatment
In addition, we investigated the CAP effects on the metastatic BrCa cell migration behavior via a transwell migration assay. Cell migration and invasion were significantly decreased in all three CAP treated BrCa cell groups (Figures 9 and 10). The number of migrated 90 s CAP treated cells decreased by nearly 3.5 folds when compared with control. The number of migrated 60 s CAP treated cells was reduced by 2 folds when compared to control.  These results exhibited that CAP play an important role in mediating migration of metastatic BrCa cells.
To further analyze cell migration and velocity after CAP treatment, a wound healing assay was performed to confirm the CAP inhibitory effects in BrCa cell migration. When started, the mean wound length was about 400 mm. After 20 hours incubation, the distance of untreated control was reduced to 0 ( Figure 11).
Consistent with the transwell migration results, CAP treated BrCa cells for 30, 60 and 90 s showed a significantly slow closure rate when compare to the control group. For 60s and 90s CAP treated BrCa cells, there was little to no wound closure after 20 h. 60 s CAP treated BrCa closed by 90.7%, whereas 90 s CAP treated BrCa only achieved 63.7% closure ( Figure 12). Consistent with cell closure percentage data, Figure 13 shows the single cell migration tracks of the first line of the scratch showed the greatly slow migration of 30 s, 60 s and 90 s CAP treated BrCa cells after 9 hours migration. In addition, CAP treated cells presented a decrease in the number of high-velocity movements ( Figure 14A), and exhibited a significant decrease in total migration distance ( Figure 14B).

Discussion
Tumor metastasis is the movement of tumor cells from a primary site to a secondary site (such as bone), which occurs in a series of steps including cancer cell migration, adhesion, proliferation, invasion and vessel formation, etc. All of those steps are regulated by complex mechanism including microenvironment interaction, gene expression and enzymes regulation. Our experiments demonstrated that through adjustment of resistance and output voltage of CAP jet, with the change of plasma emission intensity and power, a mild dose with suitable discharge power can be acquired to efficiently kill metastatic BrCa cells while keeping healthy MSC alive. Due to the observed significantly inhibited migration of BrCa cells under CAP treatment, one of the possible effects of moderate intensity and power plasma treatment may contribute to alter BrCa cell surface receptor functions [24] (such as integrin mediating cell adhesion and migration) [25]. Furthermore, the present research demonstrated that CAP treatment plays an important role in tumor therapy via immediate cell lyses [26]. Cell death in response to CAP treatment can be correlated to all aspects of the treatment condition, especially in the dose range. This is primarily a result of many likely factors including the presence of charged electrons, ROS, and RNS [27]. Furthermore,  plasma-induced DNA damage may likely initiate cellular change. Although the possible underlying mechanisms of cell death induced by CAP remain unclear, previous research has elucidated several potential mechanisms for diverse cell responses under CAP [20,28,29]. For example, ROS in CAP is reported to play a large role in cancer cell and bacterial ablation [30][31][32]. CAP induced DNA-damage is purported to be attributed to ROS-induced DNA lesions caused by point mutations [33]. Recent reports also indicated that lipid peroxidation and oxidative DNA damage in E. coli was caused by intracellular ROS formation of CAP [34]. In our study, we observed that high doses of CAP treatment can significantly decrease cancer cell viability, which are believed to be mainly caused by ROS. In addition to DNA lesions caused by ROS, DNA precursor oxidation also has the potential to influence cell apoptosis [35,36]. Also, it has been reported that CAP induced pH change of cell culture media may lead to cell apoptosis. However, our study showed that our medium pH change under the current CAP set-up was considered to be negligible. Furthermore, thermal effect on cells is virtually negligible in our study, which corresponds with our previous research [15]. Physically, CAP exposure can create temporal openings of the cell membrane, usually over a microsecond timescale, allowing for CAP specie transportation. This results in different cell responses [37]. In addition, low energetic UV photos may indirectly affect genes and cause cell death [37]. Moreover, since cancer cells proliferate much faster than normal cells, it has been reported that CAP can selectively target cancer cells by interfering with the mitotic cell cycle [18]. Thus, the observed selective effects of CAP in our study may be attributed to the significant difference in the distribution of BrCa cells and MSCs within the cell cycle, thus rendering cancer cells more sensitive to CAP treatment. Under fluorescence microscopy, specifically within the context of a live dead assay, we can clearly see that treated BrCa cells exhibited a significantly decreased viability with CAP treating.
Current results show that in vitro wound closure assay results were in line with the transmigration assay results. CAP treated BrCa cell wound areas close more slowly than untreated ones during a 20 h culture period. Interestingly, after 60 s and 90 s CAP treatment, when wounds reached close to full closure, BrCa cells in the first line of the bottom of the wound migrated in a horizontal direction instead of a vertical one. Consistent with the cell adhesion and proliferation studies, it is postulated that CAP exposure may also affect genes related to cell migration and metastasis, which play an important role in either horizontal (wound healing assay) or vertical (transmigration assay) directional migration.
Our current study focused on the initial interactions between CAP, MSCs and BrCa cells in order to optimize the parameters. The next step will be focused on devising an in vitro co-cultured model with a bone-tumor microenvironment and utilizing an in vivo model to further study CAP's effects. For future clinical applications, the CAP can be administered to any tumor locations in breast cancer patients via an endoscopic tube for endoscopic treatment. If the size of the solid tumor is large as is common in middle/late stage breast cancer patients, before conventional surgery, we expect that direct endoscopic exposure of CAP to the tumor can highly efficiently shrink the size of tumor; after surgery removal procedure, CAP can be further employed to kill existing/ remaining cancer cells and prevent further metastases.

Conclusions
In summary, this study investigated CAP generation conditions (voltage, resistances, treatment duration and frequency) for selectively killing BrCa cells. BrCa cells were found to be more sensitive to these CAP treatments than MSCs under plasma dose conditions tested. It has been demonstrated for the first time that CAP can selectively ablate metastatic BrCa cells in vitro without damaging healthy MSCs. The results suggest the great potential of CAP to be highly selective towards breast cancerous cells, resulting in tumor remediation with near complete ablation while maintaining healthy cells and tissue intact.