text:
Phylodynamic Analysis of Ebola Virus Disease Transmission in Sierra Leone

We generated genome sequences from 218 cases of Ebola virus disease (EVD) in Sierra Leone (SLE) during 2014–2015 to complement available datasets, particularly by including cases from a period of low sequence coverage during peak transmission of Ebola virus (EBOV) in the highly-affected Western Area division of SLE. The combined dataset was utilized to produce phylogenetic and phylodynamic inferences, to study sink–source dynamics and virus dispersal from highly-populated transmission hotspots. We identified four districts in SLE where EBOV was introduced and transmission occurred without onward exportation to other districts. We also identified six districts that substantially contributed to the dispersal of the virus and prolonged the EVD outbreak: five of these served as major hubs, with lots of movement in and out, and one acted primarily as a source, exporting the virus to other areas of the country. Positive correlations between case numbers, inter-district transition events, and district population sizes reaffirm that population size was a driver of EBOV transmission dynamics in SLE. The data presented here confirm the role of urban hubs in virus dispersal and of a delayed laboratory response in the expansion and perpetuation of the EVD outbreak in SLE.

Introduction
Between 2013 and 2016, the West African countries of Guinea, Liberia, Sierra Leone, Nigeria, and Mali reported the largest outbreak of Ebola virus disease (EVD) in recorded history, with 28,646 documented cases, of which 11,323 were fatal [1,2]. Sierra Leone (SLE) was one of the worst affected, with 14,124 cases and 3956 fatalities [1]. Poor case management and ineffective outbreak responses in a country previously devastated by civil war, with an inadequate healthcare system, under-resourced health facilities, unsafe burial practices, and poor infection control likely fueled the epidemic.
Molecular studies revealed that the SLE portion of the outbreak, which occurred largely in isolation with concurrent outbreaks in Liberia and Guinea, originated from the introduction of two distinct lineages (SL1 and SL2) of the Makona variant of Ebola virus (EBOV) from Guinea into eastern SLE [3][4][5]. The SL3 lineage responsible for the majority of cases in SLE, evolved later in the background of SL2, and was initially characterized by a single mutation at genome position 10,218 (G>A in the untranslated region of VP24). By July 2014, the virus had reached the capital city Freetown, although sequence data is only available from this region starting in August 2014 [6]. The Western Area (WA) division of SLE is composed of two districts: the capital Freetown (henceforth referred to as WA Urban District-WAU) and WA Rural District (WAR) ( Figure S1). By November 2014, the virus had diversified considerably and was co-circulating across multiple districts of SLE, with Waterloo (a transport hub located in the district of WAR) harboring the highest phylogenetic diversity [5,7].
A phylodynamic study of the entire West African outbreak, based on 1610 sequences, found that EBOV tended to disperse between areas in close geographical proximity, rather than spreading over large distances, and international transmission did not occur frequently into SLE [6]. Apart from the initial introductions, late back-and-forth movements between Kambia (SLE) and Forécariah (Guinea) helped prolong the outbreak in these areas, but these transmission chains did not spread further within SLE [6]. Repeated introductions into different geographic locations and widespread transmission in densely populated urban environments were the major contributors to the scale of the epidemic, with onward transmission to other administrative regions playing a role in the longevity of the outbreak [6]. This study; however, focused primarily on broad, regional patterns, and a comprehensive analysis of country-level dynamics has not yet been published for SLE. Furthermore, there are some important temporal and spatial gaps in the previous dataset analyzed for SLE, which we aimed to address in this study.
By August 2014, the WA had reported 84 probable or confirmed cases [1], at which time laboratory support was focused in the eastern part of the country, where the outbreak began, resulting in delayed reporting of results and possible underreporting of cases in western SLE. The WA eventually became the worst affected administrative district in SLE. The National Institute for Communicable Diseases Field Ebola Diagnostic Laboratory (NICD-FEDL), based in WAU, provided EVD diagnostic support to the severely stricken western part of SLE from August 2014 to June 2016 [8]. The NICD-FEDL operational period encompassed December 2014, a period for which there are only two published sequences (0.17% out of 1203 cases) from the WA of SLE, despite a peak in transmission during this time [1,6].
Here, we utilize additional genomes generated from samples collected at NICD-FEDL, along with all publicly available sequences from SLE, to conduct a complementary analysis of the dynamics of EBOV transmission between districts within SLE, and to test the hypothesis that administrative regions served either as sinks, sources, or hubs for virus dispersal relative to their population size and density.

Samples
Clinical samples in this study were derived from 7267 diagnostic submissions, consisting of buccal swabs and blood/serum, to the NICD-FEDL between August 2014 and March 2015, operating from Lakka Hospital, Freetown [8]. Sample information and metadata are summarized in Table S1.

Sample Processing, Sequencing, and Genome Assembly
Samples were stored at −20 • C or −70 • C shortly after diagnostic testing. Samples were subsequently transferred on dry ice to the NICD Biosafety Level 4 facility in Johannesburg, South Africa and stored at −70 • C.
A total of 532 samples were re-tested in South Africa by an in-house qRT-PCR, targeting the polymerase L-gene, to confirm suitability for sequencing (suitability criteria established as Ct value < 25). Viral nucleic acid was extracted from serum, or swabs immersed in Eagle's Minimum Essential Media (EMEM) culture medium (Lonza, Basel, Switzerland), using the QIAmp viral RNA kit (Qiagen, Hilden, Germany) in accordance with the manufacturer's instructions. Host ribosomal RNA was depleted as described before [9]. Double-stranded cDNA was prepared using unbiased amplification, and sequenced on the Illumina MiSeq platform [10,11]. A subset of samples with low viral loads were enriched for EBOV nucleic acids using the TruSeq RNA Access kit (Illumina) modified with probes specific for EBOV, prior to sequencing [12].
For 33 samples with insufficient volume for direct sequencing, virus isolates (one passage in VeroE6 cells) obtained in an earlier study were used for sequencing [16]. No evidence of cell culture-induced sequence changes, compared to clinical samples, were found in a short pilot experiment (results not shown). RNA was extracted from culture supernatant and processed as described above. Only assembled sequences with 95% or more coverage of the complete EBOV genome were included for further analyses in this study (n = 218).
BEAST v1.8.4 [18] was used to generate a time-structured Bayesian phylogeny. This analysis included only 1062 of the EBOV genomes from SLE; 168 genomes were excluded because they lacked district-level metadata and 16 sequences were excluded because they are known to have resulted from a re-introduction from Guinea. We used an uncorrelated relaxed molecular clock with lognormally distributed rate categories [19], along with the nonparametric Bayesian SkyGrid tree prior with 50 parameters [20]. Sites were partitioned into non-coding intergenic regions and codon positions 1, 2, and 3. The evolution of all four partitions was modelled by independent HKY substitution models with Γ4-distributed rate heterogeneity. To infer viral migration between districts (within SLE), an asymmetric continuous time Markov chain (CTMC) approach was chosen [21]. Bayesian stochastic search variable selection (BSSVS) was employed to identify strongly supported (Bayes factor > 3) migrations. An uninformative reference prior was used on the migration rate. The analysis was run five independent times, each with 100 million Markov chain Monte Carlo (MCMC) steps, sampling parameters, and trees every 10,000 generations. Between 20 and 60 million MCMC steps were discarded from each run as "burn-in", and the sampled trees were combined using Logcombiner v1.8.4.
The BEAST phylogenetic tree was visualized using Figtree v1.4.3 (http://tree.bio.ed.ac.uk/ software/figtree/). Baltic (https://github.com/evogytis/baltic) was used for parsing and visualizing maps and results from BEAST trees (code available at https://github.com/jtladner/Manuscripts/ tree/master/vanVuren_SLE-EBOV). For sub-lineage specific analyses, data was parsed from the trees generated from the BEAST analysis with the full data set. We defined a set of required members for each sub-lineage (Table S3), but we allowed for up to two additional, unspecified sub-lineage members per sampled tree in order to integrate over uncertainty in sub-lineage membership.

Statistical Analyses
Statistical analyses were performed in R and Microsoft Excel. We fit a linear model by robust regression in R v.3.4.3 using the rlm function contained within the MASS package v.7.3-47. The predict function was then used to construct 95% prediction intervals. Correlation factors (R 2 ) and p-values (one-tailed distribution t-test with two-sample unequal variance) were determined in Excel.

Ethical Statement
Permission for this study was provided by the Government of Sierra Leone, Office of the Sierra Leone Ethics and Scientific Review Committee, Directorate of Training and Research, Ministry of Health and Sanitation (version 12/10/2015, approved 30/10/2015). Clearance for the export of samples from Sierra Leone to South Africa was granted under export permit numbers PBSL/061/02/2015 and PBSL/063/02/2015 by the Pharmacy Board of Sierra Leone.

Correlation Between Cases and Population Size
We observed a strong positive correlation between population size and the number of laboratory-confirmed EVD cases per district in SLE (Table S4; Figure S2) [1,22]. Based on the best-fit slope, there was an increase of one case per 333 people. We observed one outlier district, WAR, where the actual number of confirmed cases (n = 1381) exceeded the expected number of cases within a 95% prediction interval. This finding supports previous observations that viral transmission within the transport hub and capital of WAR, Waterloo, played an important role in the amplification of the outbreak in SLE [5].

Sequencing of EBOV Genomes from SLE
A total of 218 full or near full-length EBOV genomes from individual patients from the western half of SLE were generated. The sequences were primarily derived from cases reported from WAU (n = 109), WAR (n = 34), Port Loko (n = 33), and Bombali (n = 29), but also included cases from Kambia (n = 2), Moyamba (n = 1), and Tonkolili (n = 10). Previously available sequence data from the densely populated WAU covered only 6.03% of laboratory confirmed cases, with our additional data increasing this coverage to 10.36% (Table S4).
Sequence data generated in this study included samples collected during August 2014 to March 2015. The number of sequences generated in this study from WA per month of the outbreak is summarized in Table S5. We provide additional sequence data from WAR and WAU from August   [1,6], particularly in WAU ( Figure S3). However, previously available sequence data for December 2014 from this part of SLE amounted to a single sequence from WAU, and one specified only as being from WA. Our data have increased the available number of sequences for this period 23-fold (Table S5).

Districts as Sources, Sinks, or Hubs
Using BEAST, we inferred the most likely district within SLE for each ancestral node in the phylogeny, and these reconstructions were used to infer the number of import and export events for each of the 14 districts (Table 1 and Figure 1). Overall, we detected widespread movement of EBOV within SLE, with an average of 225 movement events per sampled tree [95% highest posterior density interval (HPD95): 212-238]. This equates to approximately one inter-district movement event for every 4.7 sampled genomes. We also observed strong positive correlations between the number of imports and exports with the total number of laboratory-confirmed cases per district ( Figure S4a,b). Based on the slopes from regression analyses, there was, on average, at least one importation event per 43.5 confirmed cases, and at least one exportation event per 54.3 confirmed cases. Bombali, Kono, Port Loko, Bo, and Kenema districts all exhibited higher than average number of cases per import event, consistent with the maintenance of longer than average transmission chains. Whereas, Kambia, Moyamba, and Pujehun exhibited lower than average number of cases per import event. Kailahun stood out as a statistically significant outlier with 658 EVD cases, but zero import events from within Sierra Leone. As expected, based on the positive correlation between district population size and number of EVD cases, there was also a positive correlation between the number of import events and district population size or density (R 2 = 0.53 and 0.64 respectively) and a weak correlation between number of export events and district population size or density (R 2 = 0.18 and 0.13 respectively).
Source-sink dynamics have been used before to describe maintenance of human influenza A virus between tropic and temperate regions [23]. The source-sink model infers that there is continuous gene flow between two populations or locations, where export/import ratios can be used to characterize these as sources or sinks, or neutral where that ratio is close to one. In the context of our analysis, we defined sources as those districts from which only export events occurred; sinks as those districts where mostly import events occurred (export to import ratio of < 0.5); and hubs as those districts where significant import and export events occurred. Four SLE districts acted only as sinks, with no exportation events (Bo, Bonthe, Koinadugu, Pujehun); four districts acted mainly as sinks but some exportation events were detected (Kambia, Kono, Moyamba, Tonkolili); five districts acted as "hubs", with a substantial number of both imports and exports (Bombali, Kenema, Port Loko, WAR, WAU; and one district, Kailahun, acted only as a source with no importation events (from within Sierra Leone) detected. Geospatially, the districts in the southeast and northwest were the main sources or hubs for spread ( Figure 1D), early and late in the outbreak, respectively ( Figure 1E). EBOV was moved throughout SLE, but the high-density districts in the northwest of the country were the primary destinations for introductions. (Figure 1B, Figure S1). The southeastern districts of Kailahun and Kenema acted as the major sources of dispersal during the early stages of the outbreak ( Figure 1E). In fact, more than 75% of the inter-district movement events prior to September 2014 can be traced back to one of these two districts. Kailahun is where EBOV first entered from Guinea [3] and, based on regression analysis, this district stood out as a low import outlier relative to the number of EVD cases ( Figure S4a). This suggests that the relatively large number of cases in Kailahun was driven primarily by initial importations from outside of SLE, followed by extensive transmission within the district, as opposed to within-country movement from other districts. As expected, based on previous studies [3], Kailahun's exportations to other SLE districts occurred the earliest, starting in June, with the last exportation occurring in October 2014 ( Figure 1E). From Kailahun, the virus dispersed mostly to neighboring Kenema, but also to more distant districts, such as Bo and Kambia, to a lesser extent ( Figure 2).
Kenema also served as an important source for viral spread early in the outbreak, resulting in introductions of the virus into at least eight SLE districts (Kono, Pujehun, Bo, Bombali, Port Loko, Tonkolili, Moyamba, and WAR) between July and October 2014 (Figures 1 and 2). In fact, our regression analysis suggests that viral transmission within Kenema resulted in a significantly larger number of export events than expected based on the number of reported EVD cases ( Figure S4b). Therefore, Kenema can be regarded as an important early hub for EBOV movement in SLE, facilitating transmission between Kailahun and the rest of SLE early in the outbreak. Although most importations into Kenema occurred early in the outbreak from Kailahun, there were two late re-introduction events from WAR (October 2014) and Kono (January 2015) as the virus moved back east ( Figure 2).   Later in the outbreak, as the incidence rate in the eastern portion of the country decreased ( Figure  S3), viral transmission within four high density northwestern districts, Bombali, Port Loko, WAR, and WAU, was primarily responsible for amplification and spread of the outbreak within SLE ( Figure  1). In fact, from September 2014 to March 2015, more than 90% of all inter-district movement events originated in one of these four districts. Approximately 70% of all inter-district movement events involved the Western Area Division (WAR + WAU); however, 29% of all movement events occurred between WAU and WAR. Discounting the within-WA moves, ~59% of all other inter-district movement events involved WA as either the source or the sink.
During the peak period of EBOV transmission in SLE, WAR acted as the main hub for the spread of the virus throughout the country (Figure 1d,e). In fact, regression analysis revealed that significantly more export events occurred from WAR than expected based on the number of Later in the outbreak, as the incidence rate in the eastern portion of the country decreased ( Figure S3), viral transmission within four high density northwestern districts, Bombali, Port Loko, WAR, and WAU, was primarily responsible for amplification and spread of the outbreak within SLE (Figure 1). In fact, from September 2014 to March 2015, more than 90% of all inter-district movement events originated in one of these four districts. Approximately 70% of all inter-district movement events involved the Western Area Division (WAR + WAU); however, 29% of all movement events occurred between WAU and WAR. Discounting the within-WA moves,~59% of all other inter-district movement events involved WA as either the source or the sink.
During the peak period of EBOV transmission in SLE, WAR acted as the main hub for the spread of the virus throughout the country ( Figure 1D,E). In fact, regression analysis revealed that significantly more export events occurred from WAR than expected based on the number of confirmed EVD cases ( Figure S4b). Notably; however, nearly 70% of the exportations from WAR were into WAU (Figure 2). Waterloo city in WAR is the country's major transport hub, with major highways converging from other districts of SLE, thus serving as a gateway to and from the capital [5].
The role of WAU as an important hub for virus dispersal is not surprising, considering the city's capital status, high population density, and better access to healthcare facilities. Importations into WAU occurred as early as July 2014, with a peak in October, and these importations continued at least until March/April 2015 ( Figure 1C). The timing of the earliest peak of introductions into WAU (October 2014) suggests that high case numbers during the peak period of transmission (November/December 2014) was the outcome of extensive intra-district virus transmission following earlier introductions, rather than a direct effect of repeated introductions at that time. Late in the outbreak, in early 2015, the WAU was the main source for virus dispersal in SLE despite a decreased number of introductions; therefore, this region appears to have played an important role in prolonging the SLE portion of the outbreak ( Figure 1E). Excluding within-WA movements, the WAU served as a source for cases in four districts: Kambia, Port Loko, Bombali, and Bonthe.
Bombali and Port Loko both served as hubs for virus dispersal (Figure 1). Both experienced importation events over extended periods during the outbreak (Bombali: July 2014 to February 2015; Port Loko: July 2014 to June 2015). Interestingly, very little transmission seemingly occurred between these two neighboring districts (Figure 2). Although the first laboratory confirmed case of EVD in Bombali was detected in July 2014, a dedicated laboratory and Ebola Treatment Unit was only established in the district in December 2014 [24]. Prior to this, patient samples were being tested by laboratories in other districts of SLE. Time from sample collection to laboratory result decreased from 2.5 to 1.2 days before and after establishment of the ETU and laboratory, leading to a decrease in time from symptom onset to laboratory result from 7.3 to 5.2 days, and a subsequent sharp decrease in EVD incidence from December onwards [24]. Analysis of a database of samples submitted to NICD FEDL for EVD testing (including but not limited to the samples for which sequences were generated in this study), for which the relevant metadata was available, revealed that the average time from patient sample collection in Bombali to laboratory testing at the NICD laboratory in Freetown was 3.0 days (n = 288 submissions), compared to the overall average for all districts of 1.44 days (n = 2984 submissions). This difference was statistically significant (p < 0.001). The observations by Gleason et al. [24] and our data suggest that a delayed public health response in Bombali might have contributed to within-district transmission and, hence, its establishment as a hub for virus dispersal.
Although Kambia primarily served as a sink for EBOV transmission chains, this district played an important role during the tail end of the outbreak in SLE. More than 50% of the inter-district movement events after March 2015 were inferred to have originated in Kambia District. This included multiple re-introductions to the neighboring Port Loko District (Figure 2).

Temporal and Geographical Distribution of Lineages in SLE
The SLE EVD outbreak was dominated by seven primary sub-lineages (i.e., distinct transmission chains), all of which originated relativey early in the outbreak (between mid-July and early August 2014), and all of which circulated for at least four months (Figure 3; Figures S5-S8; Table S6). Each of these sub-lineages are defined by one or more basal substitution events (Table S7), though in rare cases sub-lineage members were found to lack these basal substitutions, either as the result of homoplasy or genome assembly errors. Lineage numbers were assigned according to the previously identified lineage nomenclature: SL1, SL2, SL3.1 (SL3. 1.1 and 3.1.2), and SL3.2 (SL3.2.1 to 3.2.5) [5]. Although there is, as of yet, no proof of phenotypic effects of these nucleotide substitutions, grouping in clades or lineages, which indicates genetic relatedness, is useful for inferring transmission chains.
The most recent common ancestor (TMRCA) for the most abundant SLE sub-lineage (SL3.2.4; TMRCA date median 16 July 2014), which accounted for 373 (29.94%) of the total number of characterized virus genomes, most likely occurred in Kenema district, and transmission was sustained over a long period (median length 383.27 days) (Table S6) responsible for back-and-forth transmission between Kambia (SLE) and Forécariah (Guinea) [6]. Only SL3.2.5 was sustained over a longer period than SL3.1.1, but this was due to detection of a single SL3.2.5 genome in September 2015, despite the last genome from this lineage, prior to that, being detected in March of the same year. This outlier genome from September 2015 (KU296454-Bombali) was closest related to a genome last detected in November 2014 (KR65322-Bombali), and might represent a case of sexual transmission. This statement warrants not only a deeper investigation into the evolution rate of this particular genome compared to expected rate, but also into the accuracy of the metadata available for this patient. The most recent common ancestor for SL3.1.2 (TMRCA date median 3 August 2014; median length 246.34 days); however, likely occurred in the Western Area, while SL3.2.1 (TMRCA date median 30 July 2014; median length 258.21 days) likely originated in Bombali (Table S6). All seven sub-lineages were detected in WAU, but were not equally abundant (Figure 4). The highest diversity of sub-lineages in the capital occurred between October 2014 and January 2015, which corresponds to the most intense transmission period based on the number of laboratory-confirmed cases [1].
2014; median length 132.47 days), and SL3.1.1 (TMRCA date median 6 August 2014; median length 396.26 days) likely occurred in Kenema. A low probability (20.7%) was found for SL3.1.1 to have originated in the WAR. Although SL3.1.1 accounted for only 8.11% of sequenced genomes, it was sustained for the second longest period. Transmission was sustained until September 2015, with all lineage SL3.1.1 sequences in August and September 2015 detected in Kambia. Viruses from this lineage were likely responsible for back-and-forth transmission between Kambia (SLE) and Forécariah (Guinea) [6]. Only SL3.2.5 was sustained over a longer period than SL3.1.1, but this was due to detection of a single SL3.2.5 genome in September 2015, despite the last genome from this lineage, prior to that, being detected in March of the same year. This outlier genome from September 2015 (KU296454-Bombali) was closest related to a genome last detected in November 2014 (KR65322-Bombali), and might represent a case of sexual transmission. This statement warrants not only a deeper investigation into the evolution rate of this particular genome compared to expected rate, but also into the accuracy of the metadata available for this patient. The most recent common ancestor for SL3.1.2 (TMRCA date median 3 August 2014; median length 246.34 days); however, likely occurred in the Western Area, while SL3.2.1 (TMRCA date median 30 July 2014; median length 258.21 days) likely originated in Bombali (Table S6). All seven sub-lineages were detected in WAU, but were not equally abundant (Figure 4). The highest diversity of sub-lineages in the capital occurred between October 2014 and January 2015, which corresponds to the most intense transmission period based on the number of laboratory-confirmed cases [1].

Discussion
EVD outbreaks before the 2013-2016 West Africa outbreak were limited to rural areas, and spread beyond affected villages was rarely reported. One of the unique features of the West Africa outbreak was the spread of EVD to densely populated urban centers, resulting in a major challenge to public health [25]. This outbreak was also characterized by extended transmission over a large geographical area, including all 14 administrative districts within SLE. We show a strong positive correlation between population size per district and number of laboratory-confirmed cases, consistent with high levels of mixing within SLE, with local population size being a primary driver of case numbers. This observation is to be expected for a well-distributed outbreak [6].
The wealth of publicly available viral genome sequences and associated metadata, along with the new sequences we generated, offered us the opportunity to perform a comprehensive, descriptive analysis of the SLE portion of the outbreak, which is complimentary to previous analyses that have either focused only on portions of SLE or on West Africa as a whole. In particular, we provide additional sequences from a peak transmission period in the Western Area, for which very little sequence data were available previously. We evaluated the role of different administrative districts as sources, sinks, or hubs for virus dispersal.
Six districts were major sources or hubs for virus dispersal within SLE. Dispersal from districts in the southeast of the country occurred early in the outbreak, after which the hubs in the northwest were responsible for further maintenance and perpetuation of the outbreak. Understanding the reasons why these particular districts acted as hubs or sources of dispersal could be valuable for the management of future outbreaks. Kailahun clearly served as a source because this is the district where the virus was initially introduced from Guinea, and subsequently transmitted mostly to the neighboring Kenema. In turn, Kenema served as an intermediate hub between Kailahun and the rest of SLE, likely due to its geographical location but also the presence of the Kenema Government Hospital. This hospital is well known for its work on another hemorrhagic fever, Lassa, and this status likely contributed to an influx of patients from Kailahun when the outbreak started.
The major role of the Western Area division, including WAU and WAR, is highlighted by the extensive within-WA transmission, and subsequent role of WAU in the dispersal of the virus late in

Discussion
EVD outbreaks before the 2013-2016 West Africa outbreak were limited to rural areas, and spread beyond affected villages was rarely reported. One of the unique features of the West Africa outbreak was the spread of EVD to densely populated urban centers, resulting in a major challenge to public health [25]. This outbreak was also characterized by extended transmission over a large geographical area, including all 14 administrative districts within SLE. We show a strong positive correlation between population size per district and number of laboratory-confirmed cases, consistent with high levels of mixing within SLE, with local population size being a primary driver of case numbers. This observation is to be expected for a well-distributed outbreak [6].
The wealth of publicly available viral genome sequences and associated metadata, along with the new sequences we generated, offered us the opportunity to perform a comprehensive, descriptive analysis of the SLE portion of the outbreak, which is complimentary to previous analyses that have either focused only on portions of SLE or on West Africa as a whole. In particular, we provide additional sequences from a peak transmission period in the Western Area, for which very little sequence data were available previously. We evaluated the role of different administrative districts as sources, sinks, or hubs for virus dispersal.
Six districts were major sources or hubs for virus dispersal within SLE. Dispersal from districts in the southeast of the country occurred early in the outbreak, after which the hubs in the northwest were responsible for further maintenance and perpetuation of the outbreak. Understanding the reasons why these particular districts acted as hubs or sources of dispersal could be valuable for the management of future outbreaks. Kailahun clearly served as a source because this is the district where the virus was initially introduced from Guinea, and subsequently transmitted mostly to the neighboring Kenema. In turn, Kenema served as an intermediate hub between Kailahun and the rest of SLE, likely due to its geographical location but also the presence of the Kenema Government Hospital. This hospital is well known for its work on another hemorrhagic fever, Lassa, and this status likely contributed to an influx of patients from Kailahun when the outbreak started.
The major role of the Western Area division, including WAU and WAR, is highlighted by the extensive within-WA transmission, and subsequent role of WAU in the dispersal of the virus late in the outbreak. Desire for better healthcare, attendance of burials/ceremonies, visiting family, and fleeing from stigmatization, combined with high population density may have contributed to the intense and protracted transmission of the virus into and within WAU, and its role as a hub. An influx of patients to the capital for better healthcare likely overwhelmed the health system and may have been one of the drivers of its role as a transmission hub.
A poor and delayed public health response likely contributed to Bombali becoming a major hub for virus transmission and dispersal. During a peak of transmission of four months between July and November 2014 in Bombali, there was no local laboratory or Ebola Treatment Unit capacity [24], necessitating the testing of samples from suspected EVD cases in Freetown laboratories, and resulting in delays, on average, just below four days. Although drivers of explosive EVD transmission are certainly complex, a poor laboratory response seems to have been a contributing factor in the timely implementation of outbreak control measures.
Interestingly, the less densely populated districts located in the northeast, interior, and south-west of SLE were almost exclusively sinks. This suggests that importations, although resulting in local transmission within the particular districts, did not result in extensive onward transmission to other districts. The sink dynamics of these rural interior districts highlights the complexity and importance of understanding human behavior in EVD dispersal. As opposed to desire for better healthcare in urban centers, this observation might be suggestive of a proportion of the rural population actively avoiding urban hospitals or control efforts.
Two recent outbreaks in Democratic Republic of Congo (DRC), in Equateur and North Kivu provinces, in 2018, raised fears of a repeat of the West Africa crisis. The most obvious similarity between these outbreaks was the relatively early introduction of cases into dense populated urban centers with ports and transport links to other major cities and neighboring countries [26]. However, the difference in the public health response was glaring, in terms of speed and approach, mostly due to hard lessons learned in West Africa, and eventually the Equateur Province outbreak was rapidly contained [27][28][29]. In addition to traditional public health measures during EVD outbreaks, the 2018 outbreaks were characterized by the early introduction of the experimental rVSV-ZEBOV vaccine for ring vaccination [26,30,31].
At the time of writing, the North Kivu EVD outbreak was still ongoing. On the 4th of September 2018, less than a month after the declaration of the outbreak, the first EVD death was recorded in the DRC city of Butembo [32]. Butembo, with an estimated population of more than one million, is an important commercial center, with an airport and large markets. Rapid mobilization of traditional public health measures, combined with new technologies, such as point-of-care testing and experimental vaccines, in affected areas that could potentially act as a source for virus dispersal, would likely represent an effective way of preventing rural-urban-rural transmission of EVD. However, even though better control measures are now available, their effective implementation is jeopardized by political unrest in the affected areas of the DRC [33].
The ultimate goal in widespread EVD outbreaks would be to skew the sink-source-hub dynamics of urban centers towards a sink role, by breaking transmission chains through traditional and next generation intervention efforts, and thereby limiting onward transmission. Our data suggest that, among other intervention efforts, reducing the time from identification of suspected cases or potential contacts in transmission chains to diagnosis could contribute to preventing highly populated urban centers from becoming virus dispersal hubs. Our data confirms that highly populated urban centers would naturally serve as sources or hubs for dispersal, but this should be minimized at all costs through prevention and control strategies. Further investigation into the factors that resulted in sink-source-hub dynamics of SLE districts could aid future intervention efforts to curb extensive geographical spread and prolonged maintenance of EVD outbreaks.

Supplementary Materials:
The following are available online at http://www.mdpi.com/1999-4915/11/1/71/s1; Table S1: Metadata of genomes generated in this study. Table S2: Metadata for complete dataset, Table S3: Required members for each sub-lineage, Table S4: Geographical distribution of EVD cases and sequences used in this study. Table S5: Temporal distribution of cases and sequences from the Western Area for the period August 2014 to March 2015. Table S6: The most recent common ancestors (TMRCA) for the seven main sub-lineages, and information on duration of circulation. Table S7: Mutations defining the seven sub-lineages of SL3. Figure S1. Map of SLE showing total case counts (laboratory confirmed as per WHO situation report) per administrative district. The districts are colored according to population density, with darker colors representing more densely populated districts. (https://www.citypopulation.de/SierraLeone-Cities.html) (http://apps.who.int/ebola/current-situation/ebolasituation-report-30-march-2016). Figure S2. Correlation plots of number of laboratory-confirmed EVD cases and district population size. The solid line represents the robust regression best-fit line and the dotted lines indicate the 95% prediction intervals. Figure S3. Laboratory-confirmed cases per district over time. Figure S4. Correlation plots of imports per cases (a) and exports per cases (b). The solid line represents the robust regression best-fit line and the dotted lines indicate the 95% prediction intervals. Figure S5. Phylogeographic reconstruction of the movement of SL3.2.4 in SLE over time. Districts are each assigned a unique color and the color of the arch is according to the source district. The circle indicates the probable location of the common recent ancestor. Figure S6. Phylogeographic reconstruction of the movement of SL3.1.2 in SLE over time. Districts are each assigned a unique color and the color of the arch is according to the source district. The circle indicates the probable location of the common recent ancestor. Figure S7. Phylogeographic reconstruction of the movement of SL3.2.5 in SLE over time. Districts are each assigned a unique color and the color of the arch is according to the source district. The circle indicates the probable location of the common recent ancestor. Figure


extfieldsofstudy: Medicine Biology

s2fieldsofstudy: Biology Medicine

text:
Helicobacter pylori upregulates TRPC6 via Wnt/β-catenin signaling to promote gastric cancer migration and invasion

Background Helicobacter pylori infection is recognized as a major risk factor for gastric cancer (GC) progression; however, the underlying molecular mechanisms have remained to be fully elucidated. Methods qPCR and Western blot were used to detect mRNA level and relative protein expression. Wound healing assay and transwell were used to determine migration and invasion of cells. Calcium imaging was used to determine calcium signaling in cells. Luciferase reporter assay and immunohistochemistry were performed. Results In the present study, it was demonstrated that H. pylori infection in GC is closely associated with the depth of tumor invasion, lymph node metastasis, tumor-nodes-metastasis stage, and distant metastasis. Migration and invasion assays indicated that H. pylori infection enhanced the migration and invasion of GC cells in a Ca2+-dependent manner. Calcium imaging was applied to detect intracellular Ca2+ and revealed that H. pylori induced an increase of intracellular Ca2+ in GC cells through release from Ca2+ stores and extracellular Ca2+ influx. Further study indicated that H. pylori infection led to an upregulation of the expression of transient receptor potential cation channel subfamily C member 6 (TRPC6) and induced an increase of Ca2+ through the TRPC6 channel. Furthermore, H. pylori increased TRPC6 transcription through the Wnt/β-catenin pathway, and Wnt/β-catenin/TRPC6 signaling was identified to be at least in part responsible for H. pylori-induced GC migration and invasion. Finally, it was observed that TRPC6 expression was significantly associated with the H. pylori infection status in GC tissues, and H. pylori infection was associated with metastasis and poor prognosis for GC patients. Conclusion The present results indicate that H. pylori causes an upregulation of TRPC6 expression through the Wnt/β-catenin pathway to promote GC progression, and this interaction may serve as a promising target for GC therapy.

Introduction
Helicobacter pylori persistently colonizes in the stomach, and has been classified as a type-I carcinogen that induces the development of gastric cancer (GC) by the World Health Organization and the International Agency for Research of Cancer. 1,2 The ability of H. pylori to colonize in the human stomach and to promote GC development may be attributed to the production of multiple virulence factors, 3,4 which modulate the pathogen-host interaction and determine the degree of bacterial pathogenicity. H. pylori exerts multiple effects on cell proliferation, apoptotic pathways, cytokine secretion, and DNA damage. [5][6][7][8] Numerous studies have indicated that H. pylori infection drives GC progression, [9][10][11] while the associated molecular mechanisms implicated in GC migration and invasion remain to be fully elucidated.
For the progression of cancer, it is required that tumor cells acquire migratory and invasive abilities, and Ca 2+ is a critical regulator of cell migration and invasion. 12 Calcium is responsible for the transmission of signals from the microenvironment into cells, and calcium signaling is an important step in cancer cell migration and invasion. 13,14 Numerous studies have indicated that Ca 2+ channels are essential for the progression of various tumor types. [15][16][17][18] Among these, transient receptor potential (TRP) channels are non-selective Ca 2+ entry channels. Several families of TRP channels exist in humans, among which TRP cation channel subfamily C member 6 (TRPC6) is greatly associated with the progression of GC. TRPC6 is significantly upregulated in human GC, and TRPC6 inhibition was reported to lead to growth suppression of GC cells. 19 However, to date, the functions of TRPC6 in H. pylori infection in GC and the underlying and mechanisms have remained largely elusive.
In the present study, H. pylori was identified to upregulate TRPC6 expression through the Wnt/β-catenin pathway, and TRPC6 then mediates the extracellular Ca 2+ influx in GC cells to finally promote GC cell migration and invasion in a Ca 2+ -dependent manner. The present study indicates that the Wnt/β-catenin/TRPC6 signaling cascade is responsible for the mediation of H. pylori-induced human GC progression, providing novel insight that may be utilized to develop strategies for preventing GC development.

Cell culture
The human GC cell line MKN45 was obtained from the Chinese Academy of Sciences, and the AGS cell line was purchased from the American Type Culture Collection (ATCC). MKN45 and AGS cells were cultured in RPMI 1640 medium (Hyclone, Waltham, MA, USA) supplemented with 10% fetal bovine serum in a humidified incubator containing 5% CO 2 at 37°C.

Analysis of TRPC6 expression in GC specimens
Tissue samples were collected from GC patients at the Department of General Surgery, Chinese PLA General Hospital. The written informed consent was obtained from all patients. This study was conducted in accordance with the Declaration of Helsinki. The use of the clinical specimens was approved by the Clinical Research Ethics Committee of Chinese PLA General Hospital. Tissues were obtained from each patient and blocked with 2.5% hydrogen peroxide in methanol. All the GC patients enrolled in our study were tested for H. pylori infection using a rapid urease test before the resection. Patients were regarded as positive for H. pylori infection if they were positive by the rapid urease test. Then, TRPC6 expression in gastric tissues was detected using the immunohistochemical staining. A polyclonal rabbit anti-TRPC6 antibody (Santa Cruz Biotechnology, CA, USA) was used for the staining. Then, the semi-quantitative analysis was performed. Briefly, images of all tissue samples were acquired with a constant set of imaging parameters on the microscope and imaging software. The images were then subjected to optical density (OD) analysis by Image-Pro Plus software. Adjustments to the background and color intensity range were performed. The intensity range selection was based on histograms, with the intensity (I) and saturation (S) set at maximum hue (H) was set at a range in which most of the brown DAB color was selected. These settings were saved and subsequently applied to all images analyzed. After defining the area of interesting, the mean OD of the selected area [integrated optical density/unit area] was determined by the software and represents the immunoreactivity of the candidate protein within tumor tissue. The acquired score of the OD was standard normalized and subjected to ImagePro Plus (Media Cybernetics, MD, USA) for further analysis. Then, patients were evaluated in terms of sex, age, tumor differentiation, depth of invasion, lymph node metastasis, TNM stage, distant metastasis, and survival.

H. pylori culture and infection in GC cells
The CagA-and Vac A-positive H. pylori strain NCTC11637 was obtained from ATCC, which was grown on Columbia agar (Oxoid, Basingstoke Hampshire, UK) plates containing 5% sheep blood and incubated at 37°C under microaerophilic conditions for 48-72 hrs. H. pylori was collected from the plates with PBS and then resuspended in RPMI 1640 medium without antibiotics. The densities of bacteria were measured by the OD at 660 nm [1 OD660=1×10 8 colony-forming units (CFU)/mL]. Cultured AGS and MKN45 cells were seeded on plates and grown to 80% confluency. Then, the cells were infected with the bacteria at a bacteria-to-cell ratio of 100:1 in the culture media.

Real-time quantitative PCR
Total RNA was extracted by TRIZOL Reagent Kit (Invitrogen, Carlsbad, CA, USA). RNA was reverse transcribed to cDNA using PrimeScript RT-polymerase (Takara, Japan). Quantitative real-time PCR was performed using SYBR Green Mix (Takara, Japan) and primers specific for TRPC6 and GAPDH. GAPDH was used as an internal control. All the experiments were performed in triplicate.

Western blotting
Cells were lysed with a lysis buffer, and protein concentration was measured using a BCA Protein Assay Kit (Beyotime, Beijing, China). Protein expression levels were detected using the Western blotting analysis. Primary antibodies: anti-TRPC6, anti-GAPDH were purchased from Cell Signal Technology (CST, USA). Then, secondary antibodies conjugated with horseradish peroxidase were used. The relative levels were quantified and normalized against GAPDH in the same sample.

Luciferase reporter vectors and Luciferase assays
TRPC6 promoter was amplified from human cDNA, and cloned into the KpnI and BglII sites of a PGL3-basic vector (Promega, USA). Similar plasmids containing mutations in the putative TCF/LEF-binding sites, named Mut, were constructed. AGS and MKN45 cells were co-transfected with the constructed luciferase vectors. Then, H. pylori was added to the cells. After 24 hrs, the cells were harvested in 1×Reporter Lysis Buffer (Catalog no.: E1960, Promega), and luciferase activities were measured using the Dual-Luciferase Reporter Assay System Kit (Promega, WI, USA). The luciferase activity for each lysate was normalized to the renilla luciferase activity.

Scratch wounding assay
GC cells were seeded into 24-well plates as a monolayer at 24 hrs. Cells were scraped with a 10 μL pipette tip. After scratching, the cells were washed to remove the detached cells. H. pylori was added to the cells. The extent of cell migration was observed, and images were captured at 0 and 24 hrs under a microscope. In each group, at least three parallel wells were utilized for testing.

Cell invasion assay
The invasive potential of GC cells treated by H. pylori was assessed using a transwell invasion assay. Matrigel Basement Membrane Matrix (BD Biosciences, San Diego, CA, USA) was added to each well and incubated at 37°C for 4 hrs to assess tumor cell invasion. Cells (approximately 1×10 4 ), in 200 μL of serum-free RPMI-1640 medium, were added into the upper compartment of the chamber. A total of 500 μL of RPMI-1640 medium with 10% fetal bovine serum were placed in the bottom compartment of the chamber. Then, H. pylori and the inhibitors were added to the upper chamber. After 24 hrs of incubation at 37°C with 5% CO 2 , the medium was removed from the upper chamber, and invading cells were fixed and stained with crystal violet. Then, the number of invasive cells were analyzed using a microscope (Olympus, Japan). Triplicate experiments were performed.

Measurement of [Ca 2+ ] cyt by calcium imaging
GC cells were cultured on coverslips and loaded with Fura-2AM (Invitrogen, NY, USA) at 37°C for 45 mins and then washed with PBS for 30 mins. After loading, cells were transferred to a standard perfusion chamber on the stage of a TMD inverted microscope (Nikon, Japan). [Ca 2+ ] i measurements were made at 340/380 nm excitation and 510 nm emission wavelengths from an SLM-Aminco spectrophotometer (Rochester, NY). F340/380 ratios were used to represent [Ca 2+ ] cyt . For the Ca 2+ -free solution, the CaCl 2 was replaced with NaCl. H. pylori and the correspond inhibitors or treatments were applied as described.

TOP/FOP flash reporter assay
A TOP/FOP flash assay was used for Wnt signaling pathway analysis. Briefly, the TCF/LEF-responsive luciferase construct was made under the control of minimal TK promoter and tandem repeats of the TCF/LEF transcriptional response element. GC cells were seeded in a 24-well plate and co-transfected with the pRL-TK plasmid (10 ng/ well) and either TOP flash plasmid or FOP flash plasmid (200 ng/well), and then treated by H. pylori and inhibitors as indicated. At 48 hrs later, the cells were analyzed using the Dual-Luciferase Reporter Assay System (Promega, E1910). The ratio of TOP Flash activity and FOP Flash activity represented the results of the activity.

Statistical analysis
Data are expressed as means ± standard error (SE). A twotailed Student's t-test or Mann-Whitney test was used to test the difference. The Kaplan-Meier curves and log-rank test were applied to assess survival. The associations between the protein expression and clinicopathological parameters were analyzed with Pearson's χ 2 test. SPSS Statistics software (version 19.0, SPSS Inc., Chicago, IL, USA) was used for statistical analysis, and the graphs were generated using GraphPad Prism 6.0 (Graphpad Software Inc, CA, USA).

H. pylori infection is correlated with the progression of GC.
To investigate the correlation between H. pylori infection and the clinicopathological features of GC patients, H. pylori colonization was examined in 360 GC tissues. As presented in Table 1, H. pylori colonization was present in 44.2% (159/ 360) of GC tissues. Analysis of the potential correlation between H. pylori infection and the clinicopathological features indicated that H. pylori infection in GC was not associated with sex, age, or degree of differentiation, while it was associated to the depth of invasion, TNM stage and lymph node or distant metastasis (Table 1). These results suggest that H. pylori infection has an important role in GC metastasis and is correlated with an advanced stage of GC.
H. pylori enhances GC cell migration and invasion in a Ca 2+ -dependent manner. To assess the effect of H. pylori infection on the migratory and invasive capacities of GC cells, AGS and MKN45 cells were infected with H. pylori NCTC 11637. The results indicated that H. pylori infection increased the migration and invasion of GC cells ( Figure 1A and B). Since Ca 2+ signaling is important in tumor cell migration, it was further explored whether Ca 2+ was involved in H. pyloriinduced GC cell migration and invasion. The results suggested that pre-treatment of GC cells with the Ca 2+ chelator BAPTA-AM significantly attenuated H. pylori-induced GC cell migration and invasion ( Figure 1A and B). Taken together, it was demonstrated that H. pylori infection promotes the migration and invasion of GC cells in a Ca 2+ -dependent manner.
H. pylori induces intracellular Ca 2+ signaling through intracellular Ca 2+ release and Ca 2+ influx. The present study further assessed whether H. pylori infection increases the intracellular Ca 2+ concentration of GC cells. Various amounts of H. pylori were added to the cell culture and the intracellular Ca 2+ concentration, [Ca 2+ ] i , was measured. As presented in Figure 2A, H. pylori induced a dose-dependent increase in [Ca 2+ ] i in GC cells, comprising a biphasic [Ca 2+ ] i increase ("peak" phase), followed by a descent to a [Ca 2+ ] i baseline ("plateau" phase).
The Ca 2+ source of the H. pylori-evoked increase in [Ca 2+ ] i in GC cells was then determined. One of the major intracellular Ca 2+ sources is the sarcoendoplasmic reticular Ca 2+ store (SERCA). GC cells were pretreated with thapsigargin, a Ca 2+ -ATPase inhibitor in the SERCA, followed by addition of H. pylori and measurement of [Ca 2+ ] i . The result indicated that thapsigargin significantly reduced the [Ca 2+ ] i peak levels, while it had no effect on the [Ca 2+ ] i plateau change ( Figure 2B).
The role of calcium channels in the cell membrane in the changes in [Ca 2+ ] i was then investigated. When H. pylori was first applied in Ca 2+ -free solutions and with pre-treatment using thapsigargin, a marked increase in [Ca 2+ ] i was observed in GC cells, and re-addition of extracellular Ca 2+ further induced an increase in [Ca 2+ ] i ( Figure 2C). Furthermore, the H. pylori-induced rise in [Ca 2+ ] i was inhibited by 2-aminoethoxydiphenyl borate (2-APB, 100 µM; Figure 2D), a non-specific inhibitor of TRP channels. Taken (Figure 3A and B). However, SKF96365 (10 µM), a blocker for TRPC channels, attenuated the H. pylori-induced rise in [Ca 2+ ] i in GC cells ( Figure 3C). Overall, these results suggest that TRPC channels may be the target calcium channels mediating H. pylori-induced Ca 2+ entry in GC cells. Since several TRPC channel subtypes are present in the cell membrane, it was attempted to identify the subtype of TRPC channel involved in H. pylori-induced Ca 2+ influx in GC cells. GC cells were infected with the H. pylori strain, and the expression of different subtypes of TRPC channel in response to H. pylori treatment was examined. TRPC6 expression was significantly increased in H. pyloriinfected GC cells ( Figure 3D and E)    known that H. pylori infection activates a series of signaling pathways, among which the Wnt/β-catenin pathway is involved in H. pylori-mediated tumor progression. [20][21][22] Therefore, the present study aimed to investigate whether H. pylori promotes TRPC6 expression through the Wnt/β-catenin pathway. The results of the present study also suggested that CagA positive H. pylori enhanced the TOP/FOP activity in GC cells ( Figure 4A and B), while CagA negative H. pylori could not enhance TOP/FOP activity ( Figure S1A). GC cells were then pretreated with XAV, an inhibitor of Wnt/β-   catenin signaling. The results indicated that XAV significantly attenuated CagA positive H. pylori-induced expression of TRPC6 in GC cells ( Figure 4C and D).
Next, it was examined whether Wnt/β-catenin controls the expression of TRPC6 transcription. A luciferase reporter assay revealed that CagA positive H. pylori enhanced the promoter activity of TRPC6, whereas XAV attenuated the H. pylori-induced promoter activity ( Figure 4E). Nevertheless, CagA negative H. pylori could not enhance TRPC6 promoter activity and expression ( Figure S1B and C). Furthermore, the gene encoding TRPC6 was analyzed using the University of California Santa Cruz (UCSC) genome browser (https://genome.ucsc.edu). One putative transcription factor (TCF) binding site was identified within the TRPC6 promoter ( Figure 4F). The chromatin immunoprecipitation assay showed than TCF could directly bind to the promoter of TRPC6 ( Figure 4G). Luciferase reporter vectors including the wild type or mutant binding site were then constructed and it was demonstrated that after infection with H. pylori, the luciferase activity of the vector with the mutant binding site was decreased compared with that of the vector with the wild-type promoter ( Figure 4H). These results clearly suggest that H. pylori increased the transcription of TRPC6 in GC cells through the Wnt/β-catenin pathway.
Wnt/β-catenin/TRPC6 signaling is responsible for H. pylori-induced GC progression. The role of TRPC6 in H. pylori-induced GC migration and invasion was further examined by performing TRPC6 knockdown, leading to a significant reduction in H. pylori-induced migration and invasion ( Figure 5A and B). Taken together, the present results indicate that Wnt/β-catenin/TRPC6 signaling is involved in H. pylori-induced GC progression. H. pylori infection increases TRPC6 expression in GC tissues and is associated with poor prognosis of GC patients. Since H. pylori infection caused an increase in TRPC6 expression in GC cells in vitro, TRPC6 expression was further examined in H. pylori-infected GC tissues. qPCR analysis of freshly resected GC tissues indicated that TRPC6 expression in H. pylori-positive GC was markedly higher than that in H. pylori-negative GC ( Figure 6A). The expression of TRPC6 was then assessed in 360 GC tissues by using immunohistochemistry. The results suggested an obvious upregulation of TRPC6 expression in H. pylori-positive GC tissues compared with that in H. pylori-negative tissues ( Figure 6B, Table S1). Furthermore, survival analysis demonstrated that H. pylori infection was highly associated with poor survival ( Figure 6C). Taken together, the present results indicated that H. pylori infection was associated with increased TRPC6 expression and poor prognosis of GC patients.

Discussion
H. pylori infection-associated chronic inflammation is an important factor for the development of GC. 23,24 Several previous studies have also suggested that H. pylori infection is closely associated with GC progression. [23][24][25][26][27] The present study indicated that H. pylori infection is associated with invasion and lymph node metastasis of GC in vivo, and it was demonstrated that H. pylori promoted GC cell migration and invasion in vitro. However, the molecular mechanisms of how H. pylori causes GC progression remain largely elusive. The present study revealed a novel mechanism, namely that H. pylori activates Wnt/β-catenin to transcriptionally upregulate TRPC6 expression and eventually induce Ca 2+ -dependent GC cell migration and invasion.
Calcium signaling is an important regulator for cell migration and invasion. 13 In the present study, calcium signaling was demonstrated to be involved in H. pyloriinduced GC cell migration and invasion. A previous study revealed that H. pylori promoted a dose-dependent [Ca 2+ ] i increase in human gastric mucous epithelial cells. 28 However, the mechanisms and exact dynamics of how H. pylori infection induces calcium mobilization in GC cells remain elusive. The present study indicated that H. pylori induced a biphasic [Ca 2+ ] i increase (transient "peak") and a [Ca 2+ ] i baseline value ("plateau" phase). The SERCA Ca 2+ store is an important intracellular Ca 2+ store. 29   There are several major classes of plasma membrane Ca 2+ -permeable channels, which mediate Ca 2+ influx in response to various activating stimuli. 8 In the present study, a panel of membrane Ca 2+ -permeable channels was screened, which have been previously reported in GC cells, and it was revealed that only TRPC channel blocker SKF96365 attenuated the H. pylori-induced increase in [Ca 2+ ] i in GC cells. Since the family of TRPC channels has seven members, 31,32 the expression of each subtype in response to H. pylori infection was further evaluated by qPCR. The results suggested that H. pylori infection induced TRPC6 expression and that H. pylori may induce the increase in [Ca 2+ ] i by activating the expression of TRPC6. The present study then attempted to investigate the possible mechanisms of how H. pylori induces TRPC6 expression.
Previous studies have indicated that H. pylori activates various intracellular pathways. [33][34][35][36] The Wnt/β-catenin pathway is critical to tumorigenesis and cancer development, 37,38 and is also responsible for H. pylori-induced gastric carcinogenesis. 11,39,40 In the present study, it was demonstrated that H. pylori activated Wnt/β-catenin signaling and induced the transcription of TRPC6 in GC cells. A search with the UCSC genome browser indicated that the TRPC6 promoter had a TCF binding site, and mutation of this binding site significantly attenuated H. pylori-induced TRPC6 promoter activity. Furthermore, inhibition of Wnt/β-catenin or TRPC6 markedly inhibited H. pylori-induced migration and invasion of GC cells. In addition, the clinical analysis revealed that TRPC6 levels in GC samples from H. pyloripositive patients were higher compared with those from H. pylori-negative patients.
In conclusion, the present study indicated that H. pylori promotes GC migration and invasion through the Wnt/βcatenin/TRPC6/Ca 2+ signaling pathway ( Figure 6D). It provides a novel molecular target for developing effective strategies to prevent GC progression.


extfieldsofstudy: Chemistry Medicine

text:
A VDAC1-Derived N-Terminal Peptide Inhibits Mutant SOD1-VDAC1 Interactions and Toxicity in the SOD1 Model of ALS

Mutations in superoxide dismutase (SOD1) are the second most common cause of familial amyotrophic lateral sclerosis (ALS), a fatal neurodegenerative disease caused by the death of motor neurons in the brain and spinal cord. SOD1 neurotoxicity has been attributed to aberrant accumulation of misfolded SOD1, which in its soluble form binds to intracellular organelles, such as mitochondria and ER, disrupting their functions. Here, we demonstrate that mutant SOD1 binds specifically to the N-terminal domain of the voltage-dependent anion channel (VDAC1), an outer mitochondrial membrane protein controlling cell energy, metabolic and survival pathways. Mutant SOD1G93A and SOD1G85R, but not wild type SOD1, directly interact with VDAC1 and reduce its channel conductance. No such interaction with N-terminal-truncated VDAC1 occurs. Moreover, a VDAC1-derived N-terminal peptide inhibited mutant SOD1-induced toxicity. Incubation of motor neuron-like NSC-34 cells expressing mutant SOD1 or mouse embryonic stem cell-derived motor neurons with different VDAC1 N-terminal peptides resulted in enhanced cell survival. Taken together, our results establish a direct link between mutant SOD1 toxicity and the VDAC1 N-terminal domain and suggest that VDAC1 N-terminal peptides targeting mutant SOD1 provide potential new therapeutic strategies for ALS.

INTRODUCTION
Amyotrophic lateral sclerosis (ALS) is a progressive and fatal neurodegenerative disease caused by the death of upper and lower motor neurons in the brain and spinal cord (Cleveland and Rothstein, 2001). The age of onset is typically between 50 and 60 years, followed by progressive paralysis and death 2-5 years after diagnosis (Mulder et al., 1986;Dorst et al., 2019). Most cases of ALS are sporadic and lack any apparent genetic linkage, although in 10% of cases, the disease is inherited in a dominant manner. About a fifth of these familial cases have been attributed to mutations in the gene encoding cytoplasmic Cu/Zn superoxide dismutase (SOD1) (Rosen et al., 1993).
To date, more than 180 different human SOD1 mutations have been identified throughout the length of the SOD1 protein that are directly linked to familial ALS (fALS), 1 including active dismutase mutants, such as SOD1 G93A and SOD1 G37R , and inactive dismutase mutants in which the mutation affects the metal-binding region, such as SOD1 G85R and SOD1 H46R (Abu-Hamad et al., 2017). The latter group of mutants are more unstable than are the former. Moreover, wild type human SOD1 (SOD1 WT ) can become misfolded and toxic, thus sharing an aberrant conformation with SOD1 mutants, when oxidatively modified (Tiwari et al., 2009;Bosco et al., 2010;Guareschi et al., 2012;Lim and Song, 2016;Medinas et al., 2018;Xu et al., 2018). The exact mechanism which drives motor neuron degeneration and disease progression remains unknown, although multiple hypotheses have been proposed to explain mutant SOD1-dependent toxic effects (Ilieva et al., 2009). These include ER stress, oxidative stress, glutamate-mediated excitotoxicity and mutant SOD1 misfolding and aggregationinduced pathology. Indeed, aggregation of misfolded SOD1 proteins is a common pathological observation among subjects with different SOD1 mutations and is, therefore, believed to be central to ALS pathogenesis (Bruijn et al., 1998;Wang et al., 2005;Prudencio et al., 2009;Abu-Hamad et al., 2017;Shvil et al., 2018).
Although predominantly a cytosolic protein, SOD1 is also found localized in other cellular compartments, including mitochondria. Both in mouse and rat models of ALS and post-mortem tissue samples from ALS patients, mutant SOD1 was found in fractions enriched for mitochondria derived only from affected but not unaffected tissues (Mattiazzi et al., 2002;Liu et al., 2004;Vijayvergiya et al., 2005;Bergemalm et al., 2006;Deng et al., 2006;Vande Velde et al., 2008). Moreover, a clear temporal correlation between disease progression and mitochondrial association was shown for different SOD1 mutants in rodent models (Liu et al., 2004). In addition, we have recently reported a clear inverse correlation between mutant SOD1 mitochondrial association in motor neuron-like NSC-34 cells and disease duration in patients carrying mutations in SOD1 (Abu-Hamad et al., 2017). 1 http://alsod.iop.kcl.ac.uk/als/ Highly purified floated mitochondria coupled with protease accessibility has demonstrated deposition of mutant SOD1 on the cytoplasmic-facing surface of spinal cord mitochondria (Liu et al., 2004;Vande Velde et al., 2008). Sensitivity to proteolysis and immunoprecipitation with specific antibodies for misfolded SOD1 further demonstrated that misfolded species of SOD1 are associated with the outer mitochondrial membrane of the spinal cord (Vande Velde et al., 2008). In addition, mutant SOD1 was proposed to interact with other components of the outer mitochondrial membrane, including Bcl-2 (Pedrini et al., 2010) and the protein import machinery (Li et al., 2010), thus affecting the corresponding functions. Importantly, this was seen only for spinal cord mitochondria but not for mitochondria isolated from unaffected tissues Li et al., 2010). More specifically, direct binding of misfolded SOD1 to the voltage-dependent anion channel-1 (VDAC1) was previously shown, causing reduction of VDAC1 conductance and channel instability, leading to inhibition of VDAC1 transport of adenine nucleotides across the outer mitochondrial membrane Magri et al., 2016). VDAC1, also known as the mitochondrial porin, is located at the outer mitochondrial membrane, where it assumes a crucial position controlling the metabolic cross-talk between the mitochondria and the rest of the cell, thus regulating the metabolic and energetic functions of mitochondria. VDAC1 is also a central player in mitochondria-mediated apoptosis and has been implicated in apoptotic-related functions, given its role as the target for pro-and anti-apoptotic Bcl2family of proteins (Shimizu et al., 1999;Arbel and Shoshan-Barmatz, 2010) and due to its function in the release of apoptotic proteins from the mitochondrial inter membrane space (Tajeddine et al., 2008;Abu-Hamad et al., 2009). VDAC1, the main VDAC isoform, is composed of 19 transmembrane β-strands forming a membrane-embedded β-barrel and a flexible amphipathic 26-residue-long N-terminal domain lying inside the pore but able to translocate from within the pore to the channel surface (Geula et al., 2012). This mobility is important for controlling channel gating but also for interactions with pro-and anti-apoptotic proteins (Abu-Hamad et al., 2009;Arbel and Shoshan-Barmatz, 2010;Shoshan-Barmatz et al., 2010Arbel et al., 2012;Geula et al., 2012). Importantly, cells expressing an N-terminally truncated form of VDAC1 are resistant to apoptosis (Abu-Hamad et al., 2009). These findings suggest that the VDAC1 N-terminal domain is required for interaction with VDAC1-associated proteins and apoptosis.
Here, we demonstrate the direct interaction of VDAC1 with mutant SOD1 and show that this interaction requires the VDAC1 N-terminal domain. Moreover, SOD1-mediated toxicity was prevented by synthetic VDAC1-N-terminal peptides. Finally, we show that a VDAC1 N-terminal peptide enhanced the survival of mutant SOD1 G93A motor neuron-like NSC-34 cells and mutant SOD1 G93A mouse embryonic stem cellderived motor neurons. These findings point to VDAC1 N-terminal peptides as offering possible novel therapeutic strategies for ALS.

Peptides
The peptides used in this study (listed in Table 1) were synthesized by GL Biochem (Shanghai, China). The peptides were first dissolved in DMSO as a 40 mM solution and then diluted 20-fold in the appropriate buffer. Peptide concentrations were determined as described previously (Shteinfer-Kuzmine et al., 2018). The final concentration of DMSO in control and peptidecontaining samples was ≤0.5%.

Protein Purification
Recombinant hSOD1 wt , hSOD1 G93A , and hSOD1 G85R were expressed in sf-9 cells and purified by hydrophobic interaction chromatography using phenyl-Sepharose 6 Fast Flow high sub (Amersham Biosciences), followed by ion exchange chromatography using a HiTrap Q-Sepharose anion exchange column (Amersham Biosciences), as described previously (Hayward et al., 2002).

VDAC1 Purification
VDAC1 was purified from rat liver mitochondria using celite:hydroxyapatite and CMC chromatography, as previously described (Gincel et al., 2001). DNA sequences encoding full-length murine VDAC1 and N-terminally truncated VDAC1 ( N-VDAC1) lacking residues 1-26 were cloned into the pET21a vector (Novagen) using the NheI/XhoI sites. Escherichia coli BL21(DE3) cells were transformed with plasmid pET21a harboring the VDAC1 or N-VDAC1 genes. Protein expression was induced for 3 h using 1 mM isopropyl-β-Dthiogalactopyranoside (IPTG; Sigma). Proteins were purified on agarose-packed nickel-nitrilotriacetic acid resin (Ni-NTA; Qiagen) in the presence of 8 M urea. Refolding of the eluted protein was performed essentially as described previously (Hiller et al., 2008). The refolded protein was further purified as above for mitochondrial VDAC1.

VDAC1 Channel Reconstitution, Recording and Analysis
The reconstitution of recombinant WT or N-VDAC1 into a planar lipid bilayer (PLB) prepared from soybean asolectin, and subsequent single channel current recordings and data analysis were carried out as described previously (Gincel et al., 2001). Currents were recorded under voltage-clamp conditions before and 5 min after the addition of 40 µg of recombinant hSOD1 wt , hSOD1 G93A , or hSOD1 G85R to the cis compartment The VDAC1-based peptides used in this study with the amino acid sequence number, molecular mass, calculated molar extinction coefficient and purity are presented. The cell-penetrating sequence is underlined. The molar extinction coefficient was calculated based on amino acid composition using the following link: http://www.biomol. net/en/tools/proteinextinction.htm.
using a Bilayer Clamp BC-535B amplifier (Warner Instrument, Hamden, CT, United States). Current amplitude histograms were prepared using AxoGraph X software. Relative conductance was determined as the average steady-state conductance at a given voltage normalized to the conductance at 10 mV, the maximal conductance. Relative conductance-voltage plots were prepared using Microsoft Excel software.
Cell Treatment With VDAC1-Based Peptides, Cell Death and XTT Analyses Apoptotic cell death was also analyzed using Acridine Orange (AcOr)/ethidium bromide (EtBr) staining (McGahon et al., 1995). Cells in 24-well plates were washed with 200 µl PBS and 10 µl of a solution containing 100 mg/ml AcOr and 100 mg/ml EtBr in PBS was added. The cells were then visualized by fluorescence microscopy (ZOE fluorescence cell imager, Bio-Rad), images were recorded and cells at early and late apoptotic stages were counted.

Immunostaining
For immunostaining, SH-SY5Y cells (4.5 × 10 4 ) were grown on sterilized coverslips in 24-well plates. 30 h post-transfection, cells were fixed using 4% paraformaldehyde (PFA; diluted in PBS) for 15 min, and then washed 3 times with PBS (5 min each wash). Cells were then permeabilized with 0.3% Triton X-100 in PBS for 5 min followed by washing with PBS. Cells were then blocked for 1 h with blocking buffer (1% BSA free fatty acids diluted in PBS). Anti-VDAC1 polyclonal antibody (ab15895) and mouse anti misfolded SOD1 (B8H10, Medimabs) were incubated at room temperature for 1-2 h in a buffer of 1% BSA free fatty acids and 0.3% Triton-X100 in PBS. Following incubation with primary antibodies, cells were washed with PBS and incubated with fluorescent conjugated secondary Alexa Flour 488 anti-rabbit and Alexa flour 647 anti-mouse antibodies. The coverslips were carefully dried and mounted on slides using Immumount (Immumount TM , Thermo). After overnight drying, images were acquired on an Olympus IX81 confocal microscope.

Mouse Embryonic Stem Cell (mESC) Cultures
Mouse embryonic stem cell lines harboring human mutant SOD1 (SOD1 G93A ) were a kind gift from Dr. Kevin Eggan (Harvard Stem Cell Institute). This cell line carries green fluorescent protein (GFP) under the control of the promoter for the motor neuron (MN)-specific transcription factor HB9 (HB9:GFP cells) (Di Giorgio et al., 2007). We used the SOD1 G93A mESC line to derive GFP + MNs.
For an in vitro differentiation assay, EBs were enzymatically dissociated after 7 days in culture with TrypLE Express (Gibco) and seeded on 0.01% poly-l-ornithine (Sigma) precoated coverslips followed by laminin (10 µg/ml, Sigma). Cells (5 × 10 4 ) were seeded on coverslip in 24-well plates with ADFNB cell medium supplemented with of Ciliary neurotrophic factor (10 ng/ml, CNTF) and Glial cell-derived neurotrophic factor (GDNF; Miltenyi Biotec). During plating, 1-20N-Ter-Antp VDAC1-peptide was added to the cultures. Half of the medium was replaced with fresh medium every second day and at the indicated time points. The cultures were fixed in 4% paraformaldehyde in PBS.

Determination of Neurite Outgrowth, Intersections Between Neurites and Survival
Stereological estimation of neurite lengths to evaluate neurite outgrowth in cultured cells was carried out as described previously (Ronn et al., 2000). The total neuritic length per cell was estimated by counting the number of GFP + soma and neurite intersections with test lines of an unbiased counting frame superimposed on images of cell cultures obtained using a 20× objective (NA 0.75) of a Nikon Eclipse E800 epifluorescence microscope equipped with a Nikon DXM1200F CCD camera. The absolute length, L, of neurites per cell was subsequently estimated from the number of neurite intersections, I, per cell by means of the equation L = (πd/2)I describing the relationship between the number of neurite intersections and the vertical distance, d, between the test lines used.

Statistical Analyses
Statistical comparisons between groups were performed by a two-tailed unpaired Student's t-test. The mean ± SEM of results obtained from at least three independent experiments are presented. The significance of differences was calculated by a two-tailed Student's t-test and is reported as * * p < 0.01. Statistical comparisons between conditions in the mESC MN assays was performed by one-way ANOVA followed by Dunnett's Multiple Comparison test against the control condition ( Figure 6B) or Tukey's Multiple Comparison test (Figures 6C,E). The confidence interval was stated at the 95% confidence level, placing statistical significance at p < 0.05. GraphPad Prism 6 was used for plotting data and statistical analysis.

VDAC1 N-Terminal Domain and a VDAC1-Derived Peptide Specifically Interact With Mutant but Not Wild Type SOD1
The interaction of purified WT and the SOD1 mutants SOD1 G93A and SOD G85R with purified VDAC1 (Figure 1A) was assayed by MST (Figures 1B-E). MST measures any variation in the thermal movement of a fluorescently labeled binding partner. The subsequent fluorescence depletion in a heated spot of the protein solution is measured as a function of increasing interacting partner concentration, with dissociation constants (K D ) values being derived from the depletion curves (Wienken et al., 2010; Figure 1D). Fluorescently labeled VDAC1 incubated with increasing concentrations of WT or mutant SOD1 (0-100 µM) showed that mutant SOD1 G93A and SOD1 G85R but not SOD1 WT bound to VDAC1 ( Figure 1B).
Next, to identify the binding site for mutant SOD1 in VDAC1, we took advantage of different VDAC1-based peptides which we have developed and previously tested (Arzoine et al., 2008;Arbel and Shoshan-Barmatz, 2010). As the N-terminal domain of VDAC1 has been the shown to interact with several proteins, such as hexokinase, Bcl-2 and Bcl-xL (Arzoine et al., 2008;Arbel and Shoshan-Barmatz, 2010;Arbel et al., 2012), we tested whether a synthetic N-terminal peptide interacts with SOD1. Accordingly, fluorescently labeled mutant SOD1 G93A or SOD1 G85R protein was incubated with increasing concentrations of the synthetic VDAC1 N-terminal peptide and changes in fluorescence were monitored ( Figure 1C). By plotting the percentage change of normalized fluorescence ( F Norm %) as a function of peptide concentration, a fitted curve yielded dissociation constants (K D ) for the three versions of SOD1 ( Figure 1D). The results showed that the VDAC1 N-terminal peptide bound both mutant SOD1 G93A and SOD1 G85R , but not to SOD1 WT (Figures 1C,D), indicating that these mutants interact specifically with the N-terminal region.
The specificity of the N-terminal peptide to mutant SOD1 was demonstrated by testing the binding of another VDAC1derived peptide, LP3. This peptide, representing the sequence of a VDAC1 loop facing the cytosol (Arzoine et al., 2008), showed significantly lower binding to mutant SOD1 than did the (1-26)-N-terminal peptide ( Figure 1E). These results show that VDAC1 and the N-terminal VDAC1 peptide specifically interact with mutant SOD1 G93A and SOD1 G85R .

Mutant but Not Wild Type SOD1 Interacts With VDAC1 and Reduces Its Channel Activity
To test whether mutant hSOD1 G93A and hSOD1 G85R binding to VDAC1 affects VDAC1 function, as well as the requirement of the N-terminal domain for such binding, full length and N-terminally truncated VDAC1 were expressed in E. coli, purified ( Figure 2M) and reconstituted into a PLB as described previously (Gincel et al., 2001).
Single-channel conductance under voltage-clamp conditions was measured as a function of time, reflecting ions passing through the channel in response to an applied voltage gradient. Current-time traces recorded at −10 or 10 mV from purified recombinant full length VDAC1 showed a stable full open state that was maintained for extended periods. Addition of hSOD1 WT , even at the highest concentration (60 µg/ml), had no effect on the current (Figures 2A,B). VDAC1 showed a bellshaped relative conductance curve as function of the voltage, with hSOD1 WT having no effect on VDAC1 channel conductance at all tested voltages, i.e., −60 to +60 mV ( Figure 2C). This was also revealed in the current amplitude histograms (Figure 2D), which showed a single channel conductance of 32 pA, at 10 mV.
In contrast to hSOD1 WT , both mutant hSOD1 G93A and hSOD1 G85R reduced the channel conductance of bilayerreconstituted VDAC1 at all tested voltages and decreased the current amplitude histograms (Figures 2E-L). hSOD1 G93A was found to be more effective in reducing VDAC1 conductance than was hSOD1 G85R when added at the same concentration, resulting in 57 and 40% inhibition of channel conductance, respectively ( Table 2).
As reported previously (Prezma et al., 2013;Shteinfer-Kuzmine et al., 2018), VDAC1 lacking the N-terminal 26 residues ( N-VDAC1) showed no voltage-dependent gating, with hSOD1 WT having no effect on channel conductance at all voltages tested or on the channel current amplitude histograms (Figures 3A-D). Moreover, and in contrast to what was observed upon SOD1 mutants interaction with VDAC1 (see Figures 2E-L), hSOD1 G93A or hSOD1 G85R had no effect on N-VDAC1 conductance at all tested voltages (Figures 3E-L). These results thus suggest that VDAC1 N-terminal is important for binding of hSOD1 G93A and hSOD1 G85R to VDAC1.
The summary of the effects of SOD1 WT , hSOD1 G93A , and hSOD1 G85R on the channel conductance of VDAC1 and N-VDAC1 is presented in Table 2.

Cell-Penetrating VDAC1 N-Terminal Peptides Inhibit Cell Death of NSC-34 Cells as Induced by Mutant SOD1 G93A
In our previous study (Shteinfer-Kuzmine et al., 2018), several novel VDAC1 N-terminal-derived peptides were designed and tested for their ability to induce cell death in cancer cells. Here, we sought to determine whether any of these peptides does not induce cell death yet can interact with misfolded mutated SOD1 and protect against SOD1 G93A -mediated cell death in neuronal cultures. Accordingly, we assessed the effects of synthetic cell-penetrating VDAC1-N-terminal-derived peptides in inducing cell death in NSC-34 cells, a mouse motor neuronlike hybrid cell line, and in the U-87MG and A549 cancer cell lines (Figures 4A-C). We considered five cell-penetrating peptides derived from the VDAC1-N-terminal domain in these studies: (1) (Table 1). FIGURE 2 | Mutant SOD1 but not wild type SOD1 interacts with VDAC1 and inhibits channel conductance. (A,B) Recombinant full length VDAC1 purified from E. coli was reconstituted into a PLB and channel currents through VDAC1, in response to a voltage step from 0 to 10 mV (A) or to -10 mV (B), before and 15-20 min after addition of 40 µg/ml (final concentration) of SOD1 were recorded. (C) VDAC1 relative conductance as a function of voltage in a 60 to -60 mV step before (•) and after addition of SOD1 WT (•). Relative conductance (conductance/maximal conductance) was determined as the average steady-state conductance at a given voltage normalized to the conductance at 10 mV, considered the maximal conductance.  3.65 ± 0.14 3.92 ± 0.08 Results are the conductance (1 M NaCl, pH 7.4, at 10 mV), calculated from experiments as in Figures 2, 3 presented as mean ± SE (n = 3).
peptides, yielded peptides that could not induce cell death (Figures 4A-C).
To test in vitro whether the different shortened versions of the VDAC1 N-terminal peptide could prevent mutant SOD1 toxicity in neurons, NSC-34 cells were transfected to express SOD1 WT or mutant SOD1 G93A protein and then treated with or without the above indicated VDAC1 N-terminal peptides. Cell viability was quantified using the XTT assay and apoptosis using acridine orange and ethidium bromide staining. Whereas SOD1 WT had not effect on cell viability, expressing SOD1 G93A reduced cell viability by 25-30%. The presence of VDAC1 N-terminal peptides reduced the toxic effect of SOD1 G93A in a concentration-dependent manner (Figures 4D-G). Importantly, this effect was observed using the three non-cell death-inducing VDAC1-derived N-terminal peptides, i.e., the (1-20)-N-Ter-Antp, (5-20)-N-Ter-Antp and (10-20)-N-Ter-Antp peptides.
In order to show that the rescue effect of the N-terminal peptide is not specific for mutant SOD1 G93A , but a general effect, we have expressed two different dismutase active mutants SOD1 G93A or SOD1 G37R to induce cell death, increasing it from 17% in the control plasmid-transfected cells to 67 and 70%, in cells expressing mutant SOD1 G93A or SOD1 G37R , respectively. Incubation with the (10-20)-N-Ter-Antp peptide, decreased SOD1 G93A and SOD1 G37R -mediated cell death by 68 and 81%, respectively ( Figure 4H).
In addition, to test whether the N-terminal VDAC1-based peptide could promote misfolded SOD1 detachment from the mitochondria and more specifically from VDAC1, SH-SY5Y neuronal cells were transfected to express two different mutated SOD1 G93A or SOD1 G37R proteins and then treated with or without the VDAC1 (10-20)-N-Ter-Antp peptide. Colocalization of misfolded SOD1 with VDAC1 was determined by immunostaining analysis using an anti-VDAC1 antibody that does not target the N-terminus region of the protein, and the B8H10 antibody specifically recognizing misfolded SOD1. In cells transfected with SOD1 G93A or SOD1 G37R , VDAC1 staining is punctuated as expected for mitochondrial localization, while misfolded SOD1 shows both punctuated staining but mostly diffused, indicating that part of the misfolded SOD1 protein is mitochondria bound (Figures 5A,C). Yet, a clear co-localization of misfolded SOD1 with VDAC1 can be observed (Figures 5A,C). However, this co-localization was greatly eliminated in cells subjected to treatment with (10-20)-N-Ter-Antp peptide (Figures 5B,D).

A VDAC1 N-Terminal Peptide Enhances the Survival of Mouse Embryonic Stem Cell-Derived Motor Neurons Expressing SOD1 G93A
The effects of increasing concentrations of the VDAC1-based (1-20)-N-Ter-Antp peptide on three characteristics of SOD1 G93Aexpressing ESC-derived motor neurons were analyzed in vitro. Specifically, we assessed neurite outgrowth, MN survival shortly after induction of final differentiation and survival of maturing SOD1 G93A -expressing MNs. To follow neuronal development and remodeling of neuronal extensions, we determined neurite outgrowth in cell cultures by tracing neurites and their branches using the stereological procedure, as described previously (Ronn et al., 2000). SOD1 G93A -expressing MNs grown in relatively pure cultures die rapidly over the 1st week (Aggarwal et al., 2017). Therefore, the first two assays were performed over the first 2 days of culture. This period is characterized by small to minimal loss of MNs, with the exception of those cells that fail to attach or which die immediately. The third assay was performed after the onset of cell death, when around 50% of cells had died.
Images showing the effects of the N-terminal peptide on neurite outgrowth of mESC-derived SOD1 G93A -expressing MNs (Figures 6A,D), as well as quantification of neurite outgrowth (Figure 6B), are presented. Cells were incubated with the indicated concentrations of the peptide and analyzed 24 h later. The results showed that adding the peptide at a 10 µM concentration significantly (p < 0.05) improved neurite outgrowth of SOD1 G93A -expressing MNs (Figures 6A,B). The effect of the peptide on the survival of mESC-derived SOD1 G93A -expressing MNs was analyzed 24 h after plating. The results showed that when added at concentrations of 5 or 10 µM, the peptide significantly (p < 0.05) improved MN density (Figures 6A,C). The reduced effect of the peptide at higher concentrations (25 µM) may result from other non-specific interactions. Moreover, the VDAC1 N-terminal peptide extended survival of the mESC-derived SOD1 G93Aexpressing mature MNs (Figures 6D,E). When five cultures were incubated with the indicated peptide concentration for 96 h, it was seen that the peptide significantly (p < 0.05) improved MN survival when administered at a 10 µM concentration (Figures 6D,E). Thus, the peptide improved survival of SOD1 G93A -expressing MNs by twofold.

DISCUSSION
The results of the present study contribute to the better understanding of mutant SOD1-mediated mitochondrial dysfunction and cellular toxicity with relevance to ALS pathogenesis. Our results demonstrate that SOD1 mutants bind FIGURE 3 | The N-terminal domain of VDAC1 is required for mutant SOD1 interaction with VDAC1 and inhibition of channel conductance. (A,B) Currents passing through bilayer-reconstituted recombinant N-terminally truncated VDAC1 ( N-VDAC1) were recorded in response to voltage step from 0 to 10 mV (A) or -10 mV (B) before and 15-20 min after the addition of 40 µg/ml (final concentration) of SOD1 WT . (C) Relative conductance of N-VDAC1 as a function of voltage in a step from 60 to -60 mV before ( ) and after addition of SOD1 WT ( ). Relative conductance (conductance/maximal conductance) was determined as the average steady-state conductance at a given voltage normalized to the conductance at 10 mV, taken as the maximal conductance. directly and selectively to the N-terminal domain of VDAC1. Both dismutase-active SOD1 G93A and dismutase-inactive SOD G85R mutants, but not wild type SOD1, bind to the VDAC1 N-terminal region. Moreover, we demonstrated that different versions of an N-terminal peptide suppressed mutant SOD1 toxicity in motor neuron-like NSC-34 cells expressing mutant The rescuing effects of the VDAC1 N-terminal-derived peptides are shown as a percentage of cell death. The significance of quantitative analysis of triplicates of different biological repeats (n = 3) was performed by Student's t-test; * * P < 0.01. (H) SH-SY5Y cells (4.5 × 10 4 cells/well in 24-well plates) were transfected with an empty plasmid or a plasmid encoding for mutant SOD1 G93A or SOD1 G37R . Twenty-four hours post-transfection, the cells were incubated for 5 h with (10-20)N-Ter-Antp peptide (20 µM) and then analyzed for apoptosis using acridine orange and ethidium bromide staining, as described previously (McGahon et al., 1995). Fluorescence microscopy images were analyzed and about 100 to 300 cells were counted for each treatment in representative microscopic fields. The significance of quantitative analysis of triplicates of different biological repeats (n = 3) was performed by one-way Anova; * * P < 0.01.

Mutant SOD1 Interacts With VDAC1 to Mediate Mitochondrial Dysfunction
The direct interaction of mutant SOD1 with VDAC1 was previously determined by immunoprecipitation using anti-SOD1, anti-VDAC1 and anti-misfolded SOD1 antibodies together with mutant SOD1 from rat spinal cord tissues or using purified proteins in a lipid bilayer system Magri et al., 2016). Now, we extended these findings and showed using both MST and VDAC1 channel conductance that the dismutase-active SOD1 G93A mutant or the dismutase-inactive SOD1 G85R mutant but not SOD1 wt specifically interact with VDAC1.
Binding of misfolded SOD1 species to mitochondrial membranes was shown to disrupt transport of metabolites required for oxidative phosphorylation, reduce membrane potential, and the activity of electron transport chain complexes (Liu et al., 2004), and to the generation of reactive oxygen (ROS) or nitrogen species with damaging effects on respiratory chain complexes (Martin et al., 2007). All these effects can be produced by interactions of mutant SOD1 with VDAC1, which mediates the transport of metabolites, ions (including Ca 2+ ) and ROS .

VDAC1 Interacts With Mutant SOD1 Through the VDAC1 N-Terminal Domain
We identified the binding sites mediating VDAC1-mutant SOD1 interactions by showing that SOD1 G93A and SOD1 G85R do not bind the VDAC1 N-terminally truncated protein, but interact with the VDAC1-derived N-terminal peptide, suggesting that mutant SOD1 interacts with the VDAC1 N-terminal domain to mediate its cell toxic effects. This finding is not surprising as the N-terminal domain of VDAC1 was shown to be the interaction site for many proteins (Arzoine et al., 2008;Abu-Hamad et al., 2009;Arbel and Shoshan-Barmatz, 2010;Arbel et al., 2012;Geula et al., 2012) and to possess an ATP-binding site (Yehezkel et al., 2007). The first eight amino acids of the VDAC1-N-terminal domain are hydrophobic in nature, thus providing a natural possible site of contact with misfolded SOD1. Wild type recombinant SOD1 remains soluble, whereas mutations in the SOD1 protein leads to exposure of certain hydrophobic residues normally buried with the protein core. These structural changes lead to misfolding and aggregation of mutant SOD1 proteins via exposure of hydrophobic residues, such that they interact with intracellular membranes, such as mitochondria, ER and others (Israelson et al., 2015).
The interaction of misfolded SOD1 with VDAC1 is further suggested by the co-localization of two different SOD1 mutants with VDAC1 at the mitochondria ( Figure 5). Furthermore, the N-terminal-derived peptide decreasing the extent of this co-localization, points to VDAC1-N-terminus as the misfolded SOD1 interaction site. Surprisingly, we have noticed a tendency for a nuclear localization of VDAC1 in cells accumulating misfolded SOD1. This phenomenon is not clear and should be further investigated.
The VDAC1 N-terminal region is proposed to move within the channel pore (Hiller and Wagner, 2009) and to translocate from the internal pore to the channel surface (Geula et al., 2012), allowing it to interact with cytosolic proteins. The multiple glycine residues ( 21 GlyTyrGlyPheGly 25 ) following this domain represent a GXXXG motif that connects the N-terminal domain to β-strand 1, and confer the flexibility required for N-terminal domain translocation out of the channel pore (Geula et al., 2012). The GXXXG motif has been shown to be involved in dimerization in proteins such as glycophorin A (Gerber and Shai, 2001), human carbonic anhydrase (Whittington et al., 2001), yeast ATP synthase (Saddar and Stuart, 2005), carnitine palmitoyltransferase (Jenei et al., 2009), and others. In VDAC1, this motif is not required for VDAC1 dimerization but it might be involved in interaction with VDAC1-associated proteins (Geula et al., 2012). Interestingly, SOD1 contains three GXXXG (residues 11-16, 32-36, and 36-41) and two GXXXXG motifs (residues 51-56 and 56-61), but their importance for the interaction of VDAC1 with mutant SOD1 proteins is unknown.
In a PLB, reconstituted VDAC1 but not the N-terminally truncated protein bound a misfolded SOD1 mutant but not wild type SOD1, leading to reduced VDAC1 conductance. It is very likely that the exposure of residues in mutant SOD1 proteins that are normally hidden leads to increased interactions between SOD1 mutants and the VDAC1 N-terminal domain. Such association would reduce the overall transit of ions through the VDAC1 pore, thereby leading to the observed reduction in VDAC1 conductance. This interaction is expected to also reduce ATP, ADP, metabolite and ROS transport, in turn leading to an inhibition of cell growth and induction of mitochondrial dysfunction and cell toxicity (Figure 7). FIGURE 7 | Mutant SOD1 binds to VDAC1 and inhibits VDAC1 activities, with addition of a VDAC1 N-terminal peptide preventing such inhibition. Schematic model showing mutant SOD1 binding to VDAC1, with the N-terminal peptide serving as a decoy. (A) Mutant SOD1 is proposed to bind to the VDAC1 N-terminal domain and inhibit VDAC1 conductance, thereby suppressing both influx and efflux of different mitochondrial metabolites and ions, including Ca 2+ , and ROS. This reduction in metabolite flux results in reduced energy production and increased oxidative stress, leading to mitochondrial dysfunction and cell death. (B) VDAC1 N-terminal-derived peptides bind mutant SOD1 and prevent its association with VDAC1, thereby preventing mitochondria dysfunction. The N-terminal peptide thus provides a new therapeutic approach for inhibiting mutant SOD1 toxicity in ALS. OMM and IMM indicate outer and inner mitochondrial membrane, respectively while IMS indicates, the intermembrane space. We have shown here that different versions of the N-terminal peptide are able to suppress the toxicity of motor neuron-like NSC-34 cells expressing mutant SOD1 G93A in a concentrationdependent manner. For these experiments, we used short versions of the VDAC1 N-terminal domain peptide which lack the GXXXG motif, as this domain is required for N-terminal peptide-induced cell death (Figures 4A-C) but not for the interaction with mutant SOD1 (Figures 4D-F).
Finally, experiments with mESC-derived SOD1 G93Aexpressing MNs showed that the (1-20)-N-Ter-Antp VDAC1-based peptide significantly improved neurite outgrowth and the survival of SOD1 G93A -expressing MNs in a concentration-dependent manner (with an optimal peptide concentration of 10 µM). This not only corroborates the insight gained from studies on NSC-34 cells but also demonstrates that the VDAC1-based peptide acts as a neuroprotectant against misfolded SOD1-mediated MN toxicity. Further, these results underline the importance of a direct misfolded SOD1-VDAC1 interaction for the appearance of misfolded SOD1 toxicity. It remains to be shown whether the VDAC1-based (1-20)-N-Ter-Antp peptide also has the same effect on fully mature human motor neurons. Indeed, the results justify more thorough in vitro and in vivo analysis of the therapeutic potential of the VDAC1-based (1-20)-N-Ter-Antp peptide in treating ALS. These include improving peptide stability, such as by introducing amino acids in the D-conformation and carrying out toxicological studies, as performed with another VDAC1based peptide specifically inducing cell death of cancer cells (Shteinfer-Kuzmine et al., 2018).
In summary, we have shown here that SOD1 mutants interact with VDAC1 through the VDAC1 N-terminal domain to exert their inhibitory effect on VDAC1 channel conductance. Moreover, we have shown that VDAC1 N-terminal-derived peptides specifically bind mutant SOD1 and inhibit mutant SOD1-induced toxicity in motor neuron-like NSC-34 cells expressing mutant SOD1 or in mouse embryonic stem cellderived motor neurons. Thus, we suggest that this VDAC1-based peptide represents a new strategy for interfering with mutant SOD1-mediated cell toxicity.