s2fieldsofstudy: Biology Psychology

text:
Epigenetic impacts of stress priming of the neuroinflammatory response to sarin surrogate in mice: a model of Gulf War illness

Background Gulf War illness (GWI) is an archetypal, medically unexplained, chronic condition characterised by persistent sickness behaviour and neuroimmune and neuroinflammatory components. An estimated 25–32% of the over 900,000 veterans of the 1991 Gulf War fulfil the requirements of a GWI diagnosis. It has been hypothesised that the high physical and psychological stress of combat may have increased vulnerability to irreversible acetylcholinesterase (AChE) inhibitors leading to a priming of the neuroimmune system. A number of studies have linked high levels of psychophysiological stress and toxicant exposures to epigenetic modifications that regulate gene expression. Recent research in a mouse model of GWI has shown that pre-exposure with the stress hormone corticosterone (CORT) causes an increase in expression of specific chemokines and cytokines in response to diisopropyl fluorophosphate (DFP), a sarin surrogate and irreversible AChE inhibitor. Methods C57BL/6J mice were exposed to CORT for 4 days, and exposed to DFP on day 5, before sacrifice 6 h later. The transcriptome was examined using RNA-seq, and the epigenome was examined using reduced representation bisulfite sequencing and H3K27ac ChIP-seq. Results We show transcriptional, histone modification (H3K27ac) and DNA methylation changes in genes related to the immune and neuronal system, potentially relevant to neuroinflammatory and cognitive symptoms of GWI. Further evidence suggests altered proportions of myelinating oligodendrocytes in the frontal cortex, perhaps connected to white matter deficits seen in GWI sufferers. Conclusions Our findings may reflect the early changes which occurred in GWI veterans, and we observe alterations in several pathways altered in GWI sufferers. These close links to changes seen in veterans with GWI indicates that this model reflects the environmental exposures related to GWI and may provide a model for biomarker development and testing future treatments. Electronic supplementary material The online version of this article (10.1186/s12974-018-1113-9) contains supplementary material, which is available to authorized users.

Background
A coalition of 34 countries deployed approximately 956,600 troops during the 1990-1991 Gulf War [1], with the majority,~700,000, from the USA. An estimated 25-32% of these veterans fulfil the requirements of a Gulf War illness (GWI) diagnosis [2]. GWI is an archetypal, medically unexplained, chronic condition characterised by persistent sickness behaviour, with neuroimmune and neuroinflammatory components.
Symptoms of GWI include fatigue, musculoskeletal pain, cognitive dysfunction, chemical sensitivities, loss of memory and sleep disruption, which can be characterised as 'sickness behaviour' [3,4]. 'Sickness behaviour' is normally a result of inflammatory response to illness or injury, which usually resolves itself over time after the initial insult is removed. Symptoms were reported within 6 months of the conflict [5][6][7].
Although the exact cause of GWI is still unknown, there is a consensus that exposure to environmental toxins is the most likely cause [8]. GWI symptoms are highly heterogeneous, and specific symptoms may be related to specific experiences: for example, Gulf War veterans exposed to nerve agents or oil well fires are at increased risk of brain cancer compared to other Gulf War veterans [9].
A leading hypothesis for the cause of GWI is that the high physical and psychological stress of combat interacted with exposure to acetylcholinesterase (AChE) inhibitors [1,4,[10][11][12][13][14][15][16][17][18][19]. Military personnel were exposed to a number of AChE inhibitors [1,20], including pyridostigmine bromide (PB), a reversible AChE inhibitor, as a prophylactic against nerve agents; sarin, soman, and related nerve agents, irreversible AChE inhibitors, which combatants were inadvertently exposed to after demolition of Iraqi supply depots, such as at Khamisiyah; organophosphate pesticides, irreversible AChE inhibitors, which were widely used to prevent pest-borne diseases and irritation [20]; permethrin, an insecticide which may inhibit AChE [10,21]; and DEET, an insect repellent and a weak AChE inhibitor which may enhance the activity of other AChE inhibitors [22]. For example, an estimated 95,000 deployed personnel were exposed to the plume from the Khamisiyah demolition, and approximately 250,000 may have been exposed to low levels of nerve agents during aerial bombardments earlier in the conflict [23]. Further, the number of nerve agent alarms heard is correlated with risk for GWI [24]. Accumulating research has indicated that deleterious health effects of exposures to psychophysiological stress [25,26] and environmental toxicants [27,28] involve epigenetic modifications that affect transcriptional regulation.
The overall objective of this study was to examine genome-wide epigenetic transcriptional modifications in the brain using an established mouse model of GWI [4,12,15,16]. Our previous research demonstrates that effects on neuroinflammatory pathways occur shortly after initial exposures. For example, pre-exposure with the stress hormone corticosterone (CORT) causes an increase in expression of specific chemokines and cytokines in response to diisopropyl fluorophosphate (DFP) [4], an irreversible acetylcholinesterase inhibitor [29] used here as a sarin surrogate. This corresponds well with the work in GWI study participants [30][31][32][33][34][35][36][37], which have shown immunological abnormalities, and a recent paper [38], which has shown specific immune-related biomarkers for GWI veterans. We hypothesized that epigenetic and transcriptomic changes upon initial exposures would identify gene pathways linked to poor health outcomes in GWI.

Animals
Adult male C57Bl/6J mice were purchased from Jackson Laboratory (Bar Harbor, ME, USA). A total of 79 animals were used for the analyses presented here. All procedures were performed under protocols approved by the Institutional Animal Care and Use Committee of the Centers for Disease Control and Prevention, National Institute for Occupational Safety and Health and the US Army Medical Research and Materiel Command Animal Care and Use Review Office. The animal facility was certified by AAALAC International. Upon receipt, the mice were housed individually in a temperature (21 ± 1°C) and humidity-controlled (50 ± 10%) colony room maintained under filtered positive-pressure ventilation on a 12-h light/ 12-h dark cycle beginning at 06:00 h. The plastic tub cages were 46 × 25 × 15 cm; cage bedding consisted of heattreated pine shavings spread at a depth of approximately 4 cm. Teklad 7913 irradiated NIH-31 modified 6% rodent chow, and water were available ad libitum.

Dosing
The dosing paradigm is presented in Fig. 1. Mice were given CORT in the drinking water (200 mg/L in 0.6% EtOH) for 4 days. This CORT regimen is known to be immunosuppressive as evidenced by decreased thymus weight [39]; thymus and spleen weights were confirmed to be decreased (> 20%) in similarly exposed animals [4,15]. On day 5, mice were given a single intraperitoneal injection of either DFP (4 mg/kg) or saline (0.9%).
Thus, there were four exposure groups: (1) saline: vehicle for 4 days, then saline injection on day 5; (2) CORT: CORT for 4 days with a saline injection on day 5;( 3) DFP: vehicle for 4 days with DFP injection on day 5; and (4) CORT + DFP: CORT for 4 days with a DFP injection on day 5.

Brain dissection and tissue preparation
Mice were killed by decapitation and the brains rapidly removed. The frontal cortex, consisting of the anterior portion of the cortex [4], and total hippocampus were dissected free-hand on a thermoelectric cold platform (Model TCP-2; Aldrich Chemical Co., Milwaukee, WI, USA) and immediately frozen at − 80°C.

Differential gene expression
Frontal cortex mRNA-seq data was generated on the Illumina HiSeq 2000 by Q Squared Solutions Expression Analysis LLC (Morrisville, NC, USA), paired-end, with a read length of 100 bp (four samples of saline, CORT and DFP, five samples of CORT + DFP). Hippocampus mRNAseq data was generated by Sickkids (Toronto, Ontario). Sequencing was carried out on the HiSeq 2500, paired-end, with a read length of 125 bp (n = 4 for all groups).
Trimmed files were aligned with the STAR aligner [43] (version 2.5.2a) in a two-pass mode. The GENCODE GRCm38.p4 assembly (mm10) and annotations were obtained from the GENCODE website [44,45] and used throughout. For the frontal cortex, there was an average of 33,547,161 reads per sample and 98.2% average mapped reads. For the hippocampus, there was an average of 35,429,647 reads per sample and 98.2% average mapped reads.
The resultant Bam files were analysed using the Genomi-cAlignments [46] and DEseq2 [47] (version 1.10.1) R packages, using the DEseq2 standard pipeline. A recent comparison study has identified this as an appropriate tool to use with replicates and when relatively large biological effects are expected [48]. Briefly, DESeq2 first fits a generalized linear model for each gene, with read counts modelled as a negative binomial distribution. An empirical Bayes approach is used for shrinkage of dispersion estimation, and the Wald test is used for significance testing, which is then adjusted for multiple corrections using the Benjamini and Hochberg method [49].
Samples were checked for similarity using Poisson dissimilarity matrix [50] with the R package PoiClaClu 1.0.2 and visualised with pheatmap 1.0.8.

Cell proportions
The R Bioconductor package DeconRNASeq [51] was used to estimate the proportion of different cell types within the sample from the RNA-seq data. Data enriched for specific CNS cell types were downloaded from the Gene Expression Omnibus (GEO) [52,53], Series GSE52564, which contains data from the Mus musculus cerebral cortex [54], to use as a reference of cell-typespecific gene expression. This RNA-seq data was trimmed and aligned and gene expression quantified as above. An expression signature for each of six cell types (astrocytes, neurons, oligodendrocyte precursor cells (OPC), myelinating oligodendrocytes (MO), microglia and endothelial cells) was obtained by finding those genes with a five-fold difference in expression in one cell type, compared to each of the others.
Using these expression signatures, the proportion of astrocytes, neurons, oligodendrocyte precursor cells (OPC), myelinating oligodendrocytes (MO), microglia and endothelial cells were estimated for each of the 17 cortex samples and 16 hippocampus samples using DeconRNA-Seq [51]. DeconRNASeq is based on a linear model of a sum of pure tissue or cell-type-specific reads of all cell types, weighted by the respective cell-type proportions. To estimate the proportions of known tissue types in a sample, DeconRNASeq solves a non-negative least-squares constraint problem with quadratic programming to obtain the globally optimal solution for estimated fractions. It is accurate down to cell types making up only~2% of the total cell populations [51].

DNA methylation modifications
DNA was extracted using the E.Z.N.A Tissue DNA Kit (VWR-Omega Bio-Tek). Bisulfite conversion was carried out using the Qiagen Epitect Fast Bisulfite Conversion Kit, and library preparation was performed using the Ovation NuGen RRBS Kit. Reduced representation bisulfite sequencing (RRBS) was carried out by the Princess Margaret Genomics Centre, part of the University Health Network, Toronto, on a NextSeq 500, using a single end, 70 base read length and multiplexed at 9-10 samples per flowcell. Samples for the cortex were 6 saline, 6 CORT, 6 DFP and 8 CORT + DFP, and for the hippocampus, 4 saline, 5 CORT, 3 DFP and 8 CORT + DFP.
RRBS fastq files were trimmed to remove adaptors and low-quality reads (q < 30) using TrimGalore version 0.4.1 [40] around Cutadapt [41] (version 1.9.1). Trimmed files were then aligned to the GENCODE GRCm38.p4 (mm10) assembly, using Bismark (v0.16.0) [55] wrapped around Fig. 1 Overview of exposure timeline. CORT + DFP exposed animals were given CORT in the drinking water for 4 days and injected with DFP on the 5th day, before being culled 6 h later bowtie2 (version 2.2.6) [56]. For the frontal cortex, the average reads per sample was 34,962,881 with 57% mapping efficiency, and for the hippocampus, the average reads per sample was 37,512,807, 63.6% mapping efficiency.
The resultant bam files were analysed with MethPipe (3.4.2) [57,58], using the suggested methods. The bisulfite conversion rate was > 98.9 for all samples. The methylation level for every cytosine site at single-base resolution is estimated as a probability based on the ratio of methylated to total reads mapped to that loci. Differential methylation was calculated by beta-binomial regression, with all batches and exposures included as factors, and CORT + DFP exposure set as the test factor. We examined both differentially methylated cytosines (DMCs) and differentially methylated regions (DMRs).

Chromatin accessibility
H3K27ac is a mark of active enhancers, strongly suggesting that genes with differential enrichment of H3K27ac will be differentially expressed [59,60]. Native chromatin immunoprecipitation sequencing (ChIP-seq) using an MNase digestion and alignment to the mm10 genome using Burrows-Wheeler Alignment was carried out by the Genome Sciences Centre, BC Cancer Agency. Samples were sequenced single end, 75 base read length on a Hiseq 2500 platform. The average read per sample was 131,727,189 with 98.9% average mapped reads. There were four samples per group, with each sample having immunoprecipitated and input DNA sequenced.
PePr (Python) and diffReps (Perl) packages were used for ChIP-seq analysis, as a recent paper by Steinhauser et al. [61] suggested that both are good tools when biological replicates are available. The results from each were compared and analysed to provide a conservative list of sites showing differential enrichment.
PePr (1.1.14) [62] was run according to the authors' suggested pipeline. PePr used a sliding window approach, shifting all reads toward their 3′ direction by half of the empirically estimated DNA fragment length and estimating window width based on the average width of the top pre-candidate peaks. The genome was then divided into consecutive widths that overlap by 50%, and the number of reads within each window was counted. This read count was then normalised based on total read count among ChIP and control sample and the relative average peak heights among ChIP samples. Read counts were modeled across replicates and between groups with a local negative binomial model. Genomic regions with less variable read counts across replicates were ranked more favourably than regions with greater variability, thus prioritizing consistently enriched regions [62]. Narrow peaks were assumed, in line with previous literature on H3K27ac (e.g. [63]). diffReps (1.55.6) [64] was run according to the authors' suggested pipeline. Bam files first had to be converted to bed files, using Bedtools (v2.26.0) [65]. Unlike PePr, diffReps used a set window size of 1000 bp for narrow peaks with a step size of 100 bp. The genome was pre-screened to remove regions with low read count to improve power and decrease computational time. Normalisation was carried out using the read count for a particular window over read count across all samples. An exact negative binomial test was used for differential analysis, which used biological replicates. p values were adjusted by the Benjamini-Hochberg method [49]. Peaks were annotated to genes using region_analysis [64].

Identifying genes unique to the CORT + DFP exposure
For several computational tools, e.g. DESeq2, PePr and dif-fReps, only direct comparisons between any two exposure groups (1v1) could be made, rather than the multifactorial comparisons available for differential DNA methylation modifications with MethPipe (RRBS). Therefore, in these cases, a series of 1v1 comparisons were made to conservatively estimate which genes were differentially expressed (RNA-seq) or enriched for H3K27ac (ChIP-seq). To use the RNA-seq data as an example, the 1v1 comparisons were carried out as:

Genes differentially expressed between CORT and
CORT + DFP and then include only the subset of genes which were not differentially expressed between saline and DFP.
(a) Changes between CORT and CORT + DFP may be due to DFP or the combination of CORT + DFP; removing those differentially expressed between saline and DFP removes those which are due to DFP alone 2. Genes differentially expressed between DFP and CORT + DFP and then include only the subset of genes which were not differentially expressed between saline and CORT.
(a) Changes between DFP and CORT + DFP may be due to CORT or the combination of CORT + DFP; removing those differentially expressed between saline and CORT removes those which are due to CORT alone 3. Intersect of the genes which appear in both list 1 and list 2 (a) Both should be genes differentially expressed due to the combined CORT + DFP exposure 4. Those genes which appear in list 3 and are differentially expressed between saline and CORT + DFP (a) Final check, as they should be different in the CORT + DFP exposure vs saline This provides a list of changes unique to the combined CORT + DFP exposure, which are not seen in either CORT or DFP exposure alone. Note that this is conservative, as genes with low but significant differential expression in either CORT or DFP alone, but a larger change in the combined CORT + DFP exposure, will be lost.

Enrichment analysis
For the sets of significant genes identified by RNA-seq, RRBS and ChIP-seq, as well as those genes that were significant in two or more of these analyses, gene set enrichment analyses were carried out. The R package clusterProfiler 3.4.4 [66] was used for Gene Ontology Biological Process (GO BP) [67,68] and KEGG pathway [69,70] enrichment analysis, with p and q value cutoffs of ≤ 0.05. Reactome pathway analyses were carried out using the ReactomePA 1.20.2 R package [71], with a p value cutoff of ≤ 0.05. All three packages reference the latest versions of their respective databases. These significantly enriched annotations were then visualised with the 'enrich-Map' function of DOSE 3.2.0 R package [72], with parameters altered to aid legibility with different numbers of enriched annotations.

Overlap between gene sets
The overlap between gene sets was visualised using the UpSetR R package [73]. This provides similar information to a Venn diagram, but in a way which makes proportions clear. Significance of overlap was determined by permutation analysis: random gene sets of the same size as our observed gene sets were taken from the same annotations, and the overlap between the two random gene sets was recorded. This was repeated 1,000,000 times, and the number of occurrences of an overlap equal to or larger than our observed overlap was divided by 1,000,000 to give an empirical p value.

Differential gene expression
Samples were clustered using a Poisson dissimilarity matrix to determine if samples from the same exposure group showed similar expression profiles. As can be seen in Additional file 1: Figure S1 and Additional file 2: Figure S2, the samples largely clustered by exposure group. The only sample that appeared to be out of place was a CORT + DFP sample in the frontal cortex that appeared intermediate between CORT and DFP alone.
In the frontal cortex, the RNA-seq analysis identified 206 GENCODE genes (204 with unique entrez IDs) that were uniquely differentially expressed in the CORT + DFP exposure group compared to all other groups (Additional file 3: Table S1). Enrichment analysis showed 12 enriched KEGG pathways (Additional file 4: Figure S3; Additional file 5: Table S2) and 24 enriched GO BP annotations ( Fig. 2; Additional file 6: Table S3). These annotations formed several broad groups related to immune response, including chemokine production, oxidative stress and steroid biosynthesis.
In the hippocampus, 667 GENCODE genes (637 with unique entrez IDs) were uniquely differentially expressed in the CORT + DFP exposure group (Additional file 7: Table  S4) compared to all other groups. Enrichment analysis showed 19 enriched KEGG pathways (Additional file 8: Figure S4, Additional file 9: Table S5) and 294 enriched GO BP annotations (Fig. 3, Additional file 10: Table S6). Similar to the frontal cortex, these annotations were grouped into several clusters (Fig. 3), including immune-related annotations (e.g. I-kappaB and NF-kappaB signalling), annotations related to nervous system differentiation, and development.
The two analyses revealed two overlapping KEGG annotations: cytokine-cytokine receptor interaction and rheumatoid arthritis and nine GO BP annotations (Additional file 11: Table S7). There are 32 genes differentially expressed under CORT + DFP priming and exposure found in both the cortex and hippocampus RNA-seq data (Additional file 12: Figure S5; Additional file 13: Table S8). Enrichment analysis of these 32 genes showed 31 enriched KEGG pathways, including annotations such as rheumatoid arthritis and cytokine-cytokine receptor interaction (Additional file 14: Table S9) and 333 GO BP annotations, including positive regulation of steroid biosynthetic process, positive regulation of chemokine production and regulation of I-kappaB kinase/NF-kap-paB signaling (Additional file 15: Table S10).

DNA methylation modifications
We next examined DNA methylation in the frontal cortex and hippocampus using RRBS to identify DNA methylation modifications associated with the exposures. The frontal cortex RRBS data showed 297 differentially methylated cytosines corresponding to 60 differentially methylated regions. Once these regions were annotated to genes (53 entrez IDs; 60 GENCODE; Additional file 16: Table S11), there was no significant enrichment for any KEGG or GO BP annotations. The hippocampus RRBS data showed 926 differentially methylated cytosines corresponding to 192 differentially methylated regions and annotated to 98 GEN-CODE genes (95 unique entrez IDs; Additional file 17: Table S12). Enrichment analysis was carried out for KEGG pathways or GO BP annotations, showing three significant GO BP enrichments: norepinephrine metabolic process (n = 3, p = 0.048, q = 0.045), cilium morphogenesis (n = 7, p = .048, q = 0.045) and cilium organization (n = 7, p = 0.049, q = 0.046). It is interesting to note in relation to the acetylcholinesterase action of DFP that a CpG site within the acetylcholinesterase gene (Ache) was significantly differentially methylated in the hippocampus (chr5: 137291317; adjp = 0.0296).
Given our hypothesis that DNA methylation modifications contribute to long-term changes in gene expression as a function of GWI exposures, this apparent lack of large, coordinated changes in DNA methylation was unexpected but could have at least two explanations. First, DNA methylation is thought to be relatively stable, and therefore, there may not have been an opportunity for substantial methylation changes to have occurred only 6 h after DFP exposure. A second possibility is that any changes were confounded by the number of different cell types within the brain. Consequently, methylation changes from any single cell type, especially cells making a small proportion of the tissue, may be lost in the 'noise'. To investigate this second possibility, we used RNA-seq data to estimate the proportions of cells in our two tissues.

Cell proportions
The estimated average proportion of each of our five cell types of interest for each exposure group is shown in Fig. 4.
In the rat cortex, neurons make up~40% of cells [74] and 44% of whole mouse brain [75]. In the human and mouse cortex, microglia make up~5% of cells [76,77]. These reports are in line with our estimates of~40-50% of cells being neurons and~4-6% of cells being microglia. This, therefore, may indicate that enriching for specific cell types, such as microglia, may enhance our ability to detect cell-type-specific methylation modifications due to these exposures. For example, currently only~1:25 RRBS counts will come from microglia.
An interesting incidental finding in the cortex was that CORT exposure, with or without co-exposure with DFP, was associated with an increase in the proportion of neurons and a decrease in the proportion of myelinating oligodendrocytes (MOs) in the frontal cortex (Fig. 4b). As we would not expect neurogenesis to occur in the frontal cortex, this suggests that the increase in the proportion of neurons is driven by a decrease in the absolute number of myelinating oligodendrocytes. A reduced number of Fig. 2 Frontal cortex RNA-seq significantly enriched gene ontology biological process annotations. Gene ontology biological process annotations significantly enriched in genes which were differentially expressed in the frontal cortex of CORT + DFP exposed mice, with groups of similar annotations highlighted oligodendrocytes would be in line with previous work in rats where CORT was shown to reduce the proliferation of oligodendrocytes [78,79]. We emphasize that these are estimated cell proportions; however, the results indicate that stereology to confirm this will be important in future studies.

Chromatin accessibility
H3K27ac is a mark of regions of the genome that are being actively transcribed [59,60]. Our RNA-seq data showed enrichment for genes involved in histone modification, suggesting that changes in chromatin accessibility may play a role in the response to the exposures. H3K27ac ChIP-seq provides an additional layer of epigenetic regulation, which may respond more quickly than methylation. It also allows an indirect examination of current transcription in the largest cell population, as H3K27ac is found at actively transcribed regions. In ChIPseq (and RRBS), every locus gives a single signal: either it is enriched for H3K27ac or it is not (or is methylation or is not). However, in RNA-seq, every locus could produce none, one or hundreds of RNA molecules, meaning that a small cell population with large changes in gene expression could mask the signal from a large population with small changes in gene expression. Therefore, using ChIP-seq will allow indirect examination of potential gene expression changes in neurons.
PePr identified 3294 GENCODE genes (3023 entrez IDs; Additional file 18: Table S13) with differential enrichment of H3K27ac, whereas diffReps identified 1518 GENCODE genes (1465 entrez IDs; Additional file 19: Table S14). The overlap between these two analyses was 563 GENCODE genes (557 entrez IDs; Additional file 20: Table S15) which were used for further analysis. However, gene annotation enrichment for each of the two gene sets (PePr and diffReps) demonstrated a large overlap in enriched annotations, suggesting that they are both finding changes in similar pathways but that the individual pathway members they find are different (85% of diffRep and 69% of PePr KEGG pathways (74) are found in both; 68% of diffReps and 54% of PePr GO BP annotations (521) are found in both).
The enrichment data indicated a clear bias towards neuronal-linked annotations, including neuronal morphology and synapse-related annotations (Additional file 21: Figure S6 and Fig. 5; Additional file 22: Table S16 and Additional file 23: Table S17). Of particular interest are the observed enrichment of GO BP annotations 'cognition' and 'learning or memory' , both of which are observed to be disrupted in GWI study participants, 'Circadian entrainment' which may relate to observed sleep disruption, and 'response to steroid hormone' , which likely relates to CORT (a steroid hormone involved in the response to stress).
These findings demonstrate that there are potential changes in neuronal-related gene expression in the frontal cortex, as was also seen in the hippocampus RNA-seq, Fig. 3 Hippocampus RNA-seq significantly enriched gene ontology biological process annotations. Top 50 gene ontology biological process annotations significantly enriched in genes that were differentially expressed in the hippocampus of CORT + DFP exposed mice, with groups of similar annotations highlighted highlighted by the fact that 33 genes were found in both the ChIP-seq frontal cortex analysis and the hippocampus RNA-seq analysis.

Overlap between genes found in different analyses
As shown in Fig. 6, there is not a large overlap in genes found between any of our analyses. However, this disparity may be partly explained from the aforementioned difference between mRNA and DNA, whereby one locus can produce many mRNA molecules, but DNA either has a modification or does not. This is reflected by the fact that the largest percentage overlap is between those genes found with RRBS and ChIP-seq, as these are both examining DNA modifications: 12% of  genes found in frontal cortex RRBS, and 16% in hippocampus RRBS, are also found in the frontal cortex ChIPseq, whereas this is only 1% and 5% for frontal cortex and hippocampus RNA-seq respectively. Similarly, 15.5% of genes found in the frontal cortex RNA-seq are also found in the hippocampus RNA-seq.
Because of the very few annotations enriched in genes differentially methylated in either of our tissues, we see very little overlap between annotations in methylation and ChIP-seq. However, we do see annotations and pathways enriched in both our ChIP-seq and RNA-seq data (Additional file 24: Figure S7 and Additional file 25: Figure S8).

Discussion
To our knowledge, this is the first transcriptome-and epigenome-wide study to examine evidence for transcriptional, chromatin and DNA methylation modifications in a model of GWI. A previous study [12] examined 84 microRNAs (MiRNAs) and global but not gene-specific changes in DNA methylation and hydroxymethylation in rats subjected to restraint stress and a protocol of PB, DEET and permethrin to simulate troops' chemical exposures. This study found differential expression of two miRNAs and in increase in global methylation in the hippocampus. Although this was an important step forward, our study was able to examine the whole transcriptome, a histone modification associated with chromatin accessibility and DNA methylation at a genome-wide, single cytosine level.
Overall, our results represent several interesting findings. First, as expected, there was a large change in the expression of immune-related genes in both the frontal cortex and hippocampus, building upon previous findings in this model [4]. Second, many genes associated with synaptic function are changed in their activity, as shown by our frontal cortex H3K27ac ChIP-seq data and hippocampus RNA-seq. These changes in gene expression seem to be subtler than those found for immune-related genes (lower expression), but differential expression of genes related to synaptic function in mice is associated with impaired memory and cognition, consistent with impairments reported by GWI suffers. Interestingly, long-term potentiation-and depression-related genes are enriched in the ChIP-seq data. Finally, we see evidence of not just a change in gene activity but of a suggested change in cell proportions. This is in line with previous work with CORT [78,79].
It is possible that microglia are responsible for this large change in immune-related gene expression, but we acknowledge that other cells, such as astrocytes, also express cytokines and chemokines (e.g. Kim et al. [80]). However, this is usually at a much lower level than that in microglia: in genes significantly differentially expressed uniquely with CORT + DFP in both the frontal cortex and hippocampus the six genes with the largest fold change are expressed in microglia (Additional file 26: Table S18). Therefore, although other cell types may be contributing to cytokine and chemokine expression, it is highly unlikely that microglia are not the cell type driving this change. Similarly, we attribute many of the transmembrane transporter-related annotations we see in the frontal cortex ChIP-seq data and hippocampus RNA-seq data to neurons; however, many of these transporters are also expressed in glial cells [54]. For this reason, future work should be carried out to isolate or enrich specific cell populations, allowing these predictions to be tested.
Our findings of altered cholinergic neurotransmitter expression after DFP exposure are in line with previous studies of low dose sarin exposure in rats [81,82]. Specifically, changes in M1 and M3 acetylcholine receptor expression were seen in the frontal cortex in a dose-dependent manner when the animals were maintained under hyperthermic conditions (i.e. a stressor; Henderson et al. [81,82]). In relation to this, we see choline-related annotations at all levels we analysed: our H3K27ac ChIP-seq analyses show changes in cholinergic synapse-related genes and differential binding related to Chrm2, the gene coding for M2; the KEGG annotation, 'choline metabolism in cancer' , is enriched in our hippocampus RNA-seq genes and our H3K27ac ChIP-seq genes; finally, in our hippocampus RRBS, we see altered methylation of a CpG site within Ache, the gene coding for acetylcholinesterase.
These findings not only link to rat models but also to GWI subjects, where differences in prefrontal cortex working memory have been previously found and have been attributed to the cholinergic system [83]. This again connects with our findings from ChIP-seq where both 'learning and memory' and 'cholinergic synapse' annotations are significantly enriched. Therefore, pre-exposure with CORT may not only be potentiating the immune response to AchE inhibitors but also causing longer term changes to the cholinergic signalling system in line with previous animal models. Human studies have been more varied, with evidence both for [84,85] and against [86] long-term changes in the cholinergic system, perhaps related to the heterogeneity of symptoms or to differences in tissues being examined. Recent work by Locker et al. [15] has shown that increased neuroinflammatory response in this model is not directly related to AChE inhibition but instead may result from changes in other aspects of signaling (e.g. epigenetic alterations in gene expression).
In relation to the potential reduction in myelinating oligodendrocytes in the cortex, this may have an effect on some of the phenotypes seen in GWI: reduced oligodendrocytes have been linked to major depressive disorder (MDD), functional consequences in neurons and mood-related symptoms in rats [87]. As this change in cell proportion would affect myelinating cells, it could also contribute to the reported alterations in white matter in GWI veterans [20,[88][89][90].
Our RNA-seq data suggests both immune dysfunction and oxidative stress. Previous reports have shown immune dysfunction and oxidative stress in other disorders with similar phenotypes to GWI such as myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) [91]. However, whether the immune response is causing oxidative stress or oxidative stress is causing an immune response is less clear. This is directly relevant, as oxidative stress contributes to the toxicity of AchE inhibitors [29] and has been suggested as a cause of GWI [1,92].
One of the six genes (Additional file 26: Table S18) we find strongly differentially expressed in both the frontal cortex and hippocampus, Tlr2, has previously been linked to sickness behaviour when expressed in hypothalamic microglia [93] and was reported to show increased expression in a model of GWI [94]. Therefore, it is of great interest in terms of the sickness-like behaviour that defines GWI. TLR activation by immune insult appears to increase Tlr2 expression [95]. We can speculate that the change in Tlr2 expression we see is upstream of the change in Il1b and TNF, as specific activation of TLR2 increases their expression [93,96]. Further, TLR stimulation increases SLC15A3 expression in dendritic cells, which in turn regulates TNF and IL-1β expression [97]. At the very least, these six, consistently, strongly upregulated genes make up a robust core of immune-and microglia-related genes for further investigation (Additional file 26: Table S18).
Our results correspond very well with those of Broderick et al. [98] in blood samples from GWI study participants, who found changes to NF-KappaB-related genes (which are enriched in our hippocampus RNA-seq data and in the genes significant in both the cortex and hippocampus RNA-seq), and in pathways under the broad theme of neuronal development and migration, which we see in our hippocampus RNA-seq and our cortex ChIP-seq. They also saw ligand-receptor interactions supporting neurotransmission, which again we see strongly in our cortex ChIP-seq data. Further, our KEGG pathway analysis of the transcriptomics data highlighted rheumatoid arthritisrelated gene enrichment. Interestingly, a study in veterans with GWI also highlighted the possibility of medication used to treat rheumatoid arthritis also being used to treat GWI [99]. The same study identified the TNF-alpha pathway (which we see in our transcriptomics) and the estrogen pathway (which we saw in our ChIP-seq) as potential drug targets [99]. Therefore, data from our model corresponds to data identified from veterans with GWI and could be used as a model to test these compounds as potential therapeutics.
A recent paper described eight potential blood biomarkers for GWI, the strongest being a 9.27-fold increase in CaMKII protein in the blood of these veterans [100]. In our study, Camk2b was found to be differentially methylated in the hippocampus and have differential H3K27ac in the cortex.
Further linking our study with human data was the finding that plasma CRP levels are increased in the blood of GWI veterans [38]. CRP is often used as a biomarker of IL-6-mediated inflammation [38,101], and IL-6 related annotations were enriched in the genes differentially expressed in both the cortex and hippocampus (Additional file 15: Table S10). This annotation was due to four genes, Il1b, Tlr2, Il1a and Tnf, four of our six most strongly and consistently differentially expressed genes.
We see minimal overlap in genes with significant differential methylation between the hippocampus and frontal cortex (two genes, Bcar3 and Tmem242). This is likely due, at least partly, to the two factors outlined above: cellular heterozygosity and a short time point after treatment. Of the two consistent genes, Bcar3 is involved in cell proliferation, which when overexpressed in breast cancer conveys estrogen resistance, and Tmem242 is a transmembrane protein with very little current annotation. We cannot detect any coordinated methylation changes (denoted by enriched annotations in our significant genes) in the frontal cortex and very few in the hippocampus. As mentioned above, we detect differential methylation of one CpG site within Ache, the gene coding for acetylcholinesterase. A number of the genes found to be differentially methylated can be linked to GWI. Examples include Lims1, which is known to be regulated by TNF and is involved in cell growth and survival [102]; Sesn1, Aplnr, Pxn and Actn1, which have been linked to ME/CFS [103,104]; Col5a3, which was identified in a rat model of GWI in the hippocampus [12]; and Slc1a2, which has been linked to ALS, a disorder GW veterans are at increased risk of [105,106]. This may indicate that networks of genes linked to GWI are only just beginning to be methylated, or that different groups of genes are methylated in different cell types. However, until these coordinated networks are elucidated, investigation of individual genes could lead to spurious associations.

Conclusions
We are able to show a range of changes in the transcriptome of this well-established mouse model, many of which reflect gene expression seen in veterans with GWI. Further, we see alteration in H3K27ac, showing potential chromatin configuration changes, which could lead to epigenetic effects with long-lasting implications. We also find differences in DNA methylation, although these are less easily interpretable than the transcriptome and H3K27ac changes.
Additional research is needed to assess whether effects of these epigenetic and transcriptional modifications on long-term health outcomes are cumulative and/or are potentiated by later exposures (e.g. infection). Notably, however, our findings reveal gene pathways known to be involved in long-term adverse health effects in GWI veterans. These results suggest that epigenetic and transcriptional regulation during the initial exposure period likely contribute to pathological outcomes in GWI. It would be important to examine these modifications in peripheral tissues from GWI veterans to ascertain whether biomarkers could be developed to predict future health outcomes.

Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. The datasets generated and/or analysed during the current study are not publicly available as they are being used for follow-up studies. They will be made publicly available after future manuscripts are published.
Authors' contributions POM, JPO, DBM and GB contributed to the conception and design of this study. JPO, DBM, LTM, KAK and JVM were involved in treatment of the mice and collection of tissue. BH, WCV, GB and DGA were involved in analysis of sequencing data and statistical analysis. DGA and POM drafted the manuscript, with all other authors critically revising the manuscript for publication. All authors read and approved the final manuscript.
Ethics approval and consent to participate All procedures were performed under protocols approved by the Institutional Animal Care and Use Committee of the Centers for Disease Control and Prevention, National Institute for Occupational Safety and Health and the US Army Medical Research and Materiel Command Animal Care and Use Review Office. The animal facility was certified by AAALAC International.

Competing interests
The authors declare that they have no competing interests.


extfieldsofstudy:
Medicine
Computer Science
Biology

s2fieldsofstudy:
Biology
Computer Science

text:
EncoMPASS: an online database for analyzing structure and symmetry in membrane proteins

Abstract The EncoMPASS online database (http://encompass.ninds.nih.gov) collects, organizes, and presents information about membrane proteins of known structure, emphasizing their structural similarities as well as their quaternary and internal symmetries. Unlike, e.g. SCOP, the EncoMPASS database does not aim for a strict classification of membrane proteins, but instead is organized as a protein chain-centric network of sequence and structural homologues. The online server for the EncoMPASS database provides tools for comparing the structural features of its entries, making it a useful resource for homology modeling and active site identification studies. The database can also be used for inferring functionality, which for membrane proteins often involves symmetry-related mechanisms. To this end, the online database also provides a comprehensive description of both the quaternary and internal symmetries in known membrane protein structures, with a particular focus on their orientation relative to the membrane.

INTRODUCTION
Protein structure determination and prediction, active site detection, and protein sequence alignment techniques all exploit information about relationships between protein structures. Hence, over the last few decades online projects for collecting and organizing structural data have proliferated, fuelled by the quasi-exponential growth of experimentally determined protein structures. However, membrane proteins, which constitute 20-30% of any given genome and >50% of all FDA-approved drug targets, are noticeably underrepresented in such databases. For example, membrane proteins account for only ∼2% of the Protein Databank (PDB) (1), of which ∼800 structures represent unique proteins as of September 2018 (http://blanco.biomol.uci.edu/ mpstruc/). Membrane proteins also appear to pose a challenge to common strategies for relating protein structures such as those used by SCOP (2) and CATH (3), whose classifications for membrane proteins are inconsistent (4). This inconsistency may reflect the distinct features that membrane proteins exhibit. In particular, while the broad functional diversity of membrane proteins is reflected in the wide divergence of their amino acid sequences (5,6), the fold space available to these proteins is restricted by their anisotropic environment. Therefore, identifying differences and commonalities in their architectures might require the development of novel strategies. In addition, the membrane orientation is a unique feature that can potentially be leveraged to improve upon such methods.
Another striking feature of membrane protein architectures is that they are abundant in symmetries and pseudosymmetries, on both intramolecular and quaternary levels (7)(8)(9). These symmetries often reflect evolution and function and can be used to predict active sites, conformational changes and mechanisms (9)(10)(11)(12). Despite their importance, the few online resources that address structural symmetries provide limited information. For example, the SymD webserver (13) can analyze a protein structure file for symmetries, but the reported results do not describe the repeating elements and the method can detect only one symmetry per structure. In the PDB database, the only information about symmetry refers to that found between chains in a complex (14), and not internal symmetries. Importantly, no available resource relates the symmetry axis to the membrane orientation, which can provide an important clue to protein function.

D316 Nucleic Acids Research, 2019, Vol. 47, Database issue
To address these issues, we created the Encyclopedia of Membrane Proteins Analyzed by Structure and Symmetry (EncoMPASS). Expanding on the offline, manuallycurated membrane protein structure classification database HOMEP (15), the EncoMPASS database is automated and uses accurate structure and sequence alignments to relate membrane protein structures. Instead of building a hierarchical classification in the style of e.g. SCOP, we constructed EncoMPASS around the structural and sequence similarities of all protein chains with similar transmembrane topologies. This results in networks of sequence and structural homologues for each protein chain. Uniquely, EncoMPASS also provides detailed information about both quaternary and internal structural symmetries, including the type of symmetry, the symmetry axes and their orientation with respect to the membrane, as well as a sequence alignment of the symmetry-related residues. The database curates results of several symmetry detection approaches and allows the user to select the analysis that best fits their objectives.
On the online version of the EncoMPASS database, the user is guided through the data by graphical and interactive interfaces, such as: 3dmol representations (16) of each membrane protein complex, each chain, and all symmetryrelated structural repeats; intuitive, interactive graphs that elucidate sequence and structure relationships among the database entries; and search bars that allow for either a straightforward or a more tailored exploration of the database. All data used to generate the displayed content are also available for download.

DATABASE CONTENT
To maximize the accuracy of our structural analysis, the En-coMPASS database collects crystallographic structures of membrane proteins with resolution ≤3.5Å. EncoMPASS uses the manually-curated Orientations of Proteins in Membranes (OPM) database (17) as the primary source of structure coordinates that are reliably oriented relative to the predicted lipid bilayer. However, the information in the OPM coordinate files is occasionally incomplete or inconsistently formatted, or the biological assembly can be in disagreement with the one reported in the PDB. Indeed, predicting the correct biological assembly of membrane protein structures is a major challenge (bioRxiv: https://doi.org/10. 1101/391961), leading to discrepancies between the assignments in the PDB (1), PDBTM (18) and OPM databases. To maximize the accuracy of the predicted biological assembly we follow the strategy implemented in PDBTM (18). First, each structure is classified as either potentially problematic or plausibly correct. Then, for each member of the former set, we check whether a structure in the latter set can serve as a template for the biological assembly. If no template is found, we assume the biological assembly described in the PDBTM web server. Thus, the coordinate file published on EncoMPASS is a revised version of the coordinate file in the PDB database, with its biological assembly optimized where necessary, and with membrane boundaries estimated by the algorithm underlying the OPM database, PPM (19). As of September 2018, EncoMPASS contains 2344 coordinate files, corresponding to 67% of all membrane protein entries present in the PDB; the majority of the excluded entries do not meet the resolution threshold.
All structural and sequence similarity networks presented in EncoMPASS are based on a large set of pairwise structure alignments between individual chains of proteins, carried out using Fr-TM-Align (20). However, aligning two structures with very different topologies can force one of the two structures to fragment excessively and thereby produce a fit with a biologically meaningless alignment. Therefore, we only compare structures with similar numbers of transmembrane regions. The quality of the structural alignment is measured using the root mean squared deviation (RMSD) of the aligned C␣-atoms, as well as the Template-Modeling score, or TM-score (21,22). The TM-score is the metric we use to establish a structural relationship between two proteins. Both the TM-score and RMSD are global scores, but we also wish to identify the most structurally conserved regions of a membrane protein chain. To this end, we also report the per-residue C␣-C␣ distances between a given chain and all others aligned to it.
Protein sequence similarity provides complementary information to structural similarity and can highlight, for example, proteins undergoing conformational changes during their function. Thus, independent of the structural alignment, the sequences of all the same pairs of topologicallyrelated proteins are aligned with MUSCLE (23).
For each protein complex and chain we also report the results from a comprehensive symmetry analysis. First, En-coMPASS curates the results of two standard algorithms for symmetry detection, CE-Symm 2.2 (bioRxiv: https:// doi.org/10.1101/297960) and SymD 1.61 (24). In addition, the database includes results from a multi-step symmetry detection (MSSD) procedure (bioRxiv: https://doi.org/10. 1101/391961). The MSSD method processes structures using CE-Symm and SymD with a range of customized parameters and then, among other considerations, filters the results based on the position of the detected structural repeats with respect to the membrane. MSSD can report multiple symmetries within either a complex or a chain; the corresponding symmetry axes and transformations; the residue ranges and multiple alignment of the structural repeats; and the orientation of the repeats and the symmetry axes relative to the membrane. Finally, taking advantage of the networks of structural and sequence similarities described above, the database reports so-called inferred symmetries. These are obtained by assuming the presence of repeated or symmetric elements in a given chain based on their detection with the MSSD method in related protein chains. The results of this final analysis allow for informative comparisons of the degree of symmetry in different conformations of the same protein or between different proteins in the same protein family.
A detailed description of the procedures used to create the dataset is provided in (bioRxiv: https://doi.org/10.1101/ 391961).

DATABASE ACCESS
The EncoMPASS web server is hosted at http://encompass. ninds.nih.gov and is equipped with a a simple PDB-code search as well as a more versatile search tool that allows the user to explore the many features of the database. The complete list of search criteria organized by category is reported in Supplementary Table S1.
Information is provided for each PDB entry on two levels: entire protein complexes ( Figure 1) and individual transmembrane chains (Figure 2). The pages for each protein complex are divided into four sections, containing a general description of the protein, a summary of its transmembrane chains and their structure relationships, and descriptions of the symmetries obtained by the standard symmetry detection algorithms and by the MSSD procedure.
The first section of each whole structure page reports general information about the complex and provides links to the corresponding entry in the PDB and OPM databases, as well as a graphical interface that allows real-time visualization of the structure from different perspectives. This section also provides a button for downloading the coordinates of the protein. Note that the displayed structure is not identical to that in the PDB: for example, the biological assembly might differ, and the membrane boundaries are highlighted by two planes of pseudo-atoms.
The second section contains a table summarizing the characteristics of each of the transmembrane chains in the complex. This includes: the number of transmembrane regions, the number of sequence and structural neighbors, and the symmetry order of any internal symmetry detected by either the MSSD, CE-Symm or inferred symmetry methods in that order of preference. By clicking on the transmembrane chain identifier, the user is redirected to the relevant page for that chain, while the sequence-, structure-and all-neighbors fields are linked to tables describing the corresponding structures.
The third and fourth sections of each complex page are dedicated to symmetry analyses. The third section summarizes the results of the two standard symmetry detection algorithms, CE-Symm and SymD. CE-Symm can detect multiple symmetries for the same structure, so the user can scroll through the different symmetries. Each symmetry is represented using different colors for each repeat and a black line for the symmetry axis. By comparison, SymD can detect only one symmetry axis in a structure and its output does not indicate the boundaries of each repeat. Hence, in the 3dmol visualization of the SymD results, all residues that are related by the detected symmetry axis are colored blue. The raw data and a PyMOL visualization script can be downloaded directly from this section. It is important to note that the output of the two symmetry recognition programs does not distinguish between quaternary and internal symmetries and, therefore, the results for complexes might include information on both. On the other hand, the subsequent section (the fourth section) contains results from our MSSD procedure focused exclusively on quaternary symmetries. The layout of the MSSD results section mirrors that of the CE-Symm results, but the reported data also includes classification of the topology of the repeats with respect to the membrane (antiparallel or parallel), as well as the angle of each reported symmetry axis to the membrane normal. Finally, the user can download a PyMOL script that allows visualization of the superposition of all repeats.
The web pages for individual chains are extended versions of the whole structure pages, with a total of five sections (Figure 2). The second section on each chain page is notable, in that it includes visual representations of the analysis of sequence and structure relationships. Three graphs are presented. The first graph illustrates the structural similarity between the chain of interest and its close structural neighbors (TM-score ≥ 0.6), where the distance between any two points is proportional to their similarity (i.e., inversely proportional to the TM-score). Since large structural differences could be explained simply by a difference in size, each point on this plot is colored in shades from blue to red to show the (greater or fewer, respectively) number of transmembrane regions relative to that of the chain of interest. Hovering over the points with the cursor brings up their PDB identifiers, while clicking on a point navigates to the EncoMPASS page describing the corresponding protein chain. All the information provided in this first plot, however, is focused on structural similarity. The user may also be interested in the sequence similarity between protein chains that have related structures, or conversely, the structural variability among entries with similar sequences. Thus, in the second graph, the structural similarity according to the TM-score is plotted against the sequence identity from the MUSCLE pair-wise sequence alignment, for every compared chain. The contour lines on the graph illustrate the abundance of pairs of compared structures across the entire database. Finally, to illustrate how different regions of the protein chain compare, the third graph provides the structural similarity (i.e. distance between C␣ atoms) as a function of residue number for each structure alignment. All similarity measures and structure alignment outputs for the protein chain of interest are also downloadable.
Sections three and four of the individual chain pages follow the same layout as those on the complex page, but obviously only report symmetries within the given chain. In the fifth section, each chain page includes results from an additional procedure that relies on the structural and sequence neighbors to infer symmetry. That is, if the chain has a neighbor, for which the MSSD analysis detected a more extensive symmetric relationship, the symmetry of that neighbor is used as a template for mapping out the repeats in the current chain. The symmetry axis and structural similarity between these repeats are reported in a similar fashion as in the MSSD results section.
The EncoMPASS dataset is created and updated automatically using a set of python libraries available at https://www.github.com/EncoMPASS-code/EncoMPASS. All data entries and processed PDB files can be downloaded as a single file from http://encompass.ninds.nih. gov/downloads. The online database for EncomPASS is designed using the Oracle Relational Database System, and utilizes the Java 2 Enterprise Edition platform, the Spring Framework, and JavaScript libraries. Representative page for the whole structure of PDB entry 3M75 (a TehA homologue), including a brief summary of the characteristics of each of its transmembrane subunits (left) and analysis of the symmetries within the complex (right) using two standard symmetry detection algorithms as well as the multi-step symmetry detection (MSSD) method. Note that on these pages the MSSD section contains results only if quaternary symmetries (i.e., involving multiple chains) have been detected. According to all three methods (CE-Symm, SymD, and MSSD), the TehA homologue forms a three-fold circularly (C3) symmetric assembly. Note that in the molecular viewer in the SymD section, the repeats are not differentiated by color, since SymD does not report residue ranges for the symmetry-related repeats. According to the MSSD results, the axis of this C3 symmetry is perpendicular to the membrane plane, i.e., the protomers have parallel topologies.

CASE STUDY: INVESTIGATING SODIUM-DEPENDENT TRANSPORT
As an illustration of how EncoMPASS can be used to generate function-related hypotheses about membrane proteins, we consider the example of the sodium-coupled betaine symporter BetP. BetP uses the translocation of two sodium ions to energize the transport of each betaine molecule (25). The available structural data for BetP, however, only indicate the location of one of the two required sodium ions, at the so-called Na2 site, for example in PDB identifier 4AIN, chain B) (26). Unfortunately, no evidence has been found for the other sodium ion (Na1) in the six available structures of BetP. Interestingly, however, EncoMPASS reports the presence of two inverted-topology structural repeats in the BetP protomer ( Figure 3A, 4AIN B). Hence, examining the region of the protein that is structurally related to the Na2 binding site may reveal residues contributing to the Na1 site. This task is made straightforward by using the structural alignment of the repeats provided in the MSSD section of EncoMPASS ( Figure 3B). Mapping the related residues onto the structure narrows the list of candidate residues that could interact with the ion to those within one or two helical turns and on the same face of the helix, which could contribute either side chain or backbone groups ( Figure 3C). Residues S376, F380, T246, T250 and S253 satisfy these requirements, and hence, might contribute to the binding site. Khafizov et al. applied similar reasoning to arrive at this prediction, which they then tested using molecular dynamics simulations, as well as biochemical, biophysical and electrophysiological experiments, to conclude that F380, T246 and T250 are indeed important for sodium binding (12).
To understand further the principles of sodium-coupled substrate transport, we can use EncoMPASS to investigate the conservation of the BetP sodium binding sites in related structures. In the polar plot provided for 4AIN chain B (Figure 3d), two clusters are immediately visible near the origin of the graph, and these clusters correspond to structures of two betaine/carnitine/choline transporter (BCCT) family members, BetP and CaiT. On the edges of the plot, we identify more distantly-related structures such as those of LeuT, ApcT, AdiC and Mhp1. Then, using the table of 'All Neighbors' of 4AIN, chain B, we can scroll through a list of all related structures (Figure 3E), select a representative for each protein and extract the relevant structural alignments from the downloadable Superpositions file. Us-Nucleic Acids Research, 2019, Vol. 47, Database issue D319 Figure 2. Example EncoMPASS online database chain page. Representative page for chain A of PDB entry 3M75. These pages follow the same outline as for the complex in Figure 1, but contain analysis specific to each chain, including sequence and structural relationships, and any detected internal symmetries. Chain A of 3M75 has 39 sequence and structure neighbors. (Visualization of the polar or scatter plots, shown here as clickable icons, reveals that these structures are all closely related, both structurally and by sequence.) In terms of the internal symmetry, all methods agree that chain A contains a five-fold circular pseudo-symmetry (C5), although there are some differences in the details. As in the pages for complexes, the MSSD section reports the relationship of the symmetry axis to the membrane plane. The last section ('Symmetry Inferred from Neighbors'), only available for chains, reports results based on inference from putative structural homologues, in cases where MSSD detected a more comprehensive symmetry than in the chain of interest. By inferring the symmetry information from another TehA-homologue structure, 3M72 chain A, this approach detected a symmetry relationship for 3M75 chain A with greater coverage (86% of the chain, compared to 77% when using MSSD) and a higher similarity (TM-score 0.70 rather than 0.68).
ing these alignments ( Figure 3F), we can examine which proteins contain sodium-coordinating residues (e.g., LeuT and Mhp1), and which do not (e.g., CaiT), and correlate those structural features with their known coupling mechanisms (12).
These analyses of BetP provide a clear demonstration of the value of EncoMPASS for delving into the structurefunction relationships of membrane proteins, including by leveraging information relating to their symmetry and their structural neighbors.

CONCLUSION
The online EncoMPASS database curates a wide variety of membrane protein structural data, which we expect will appeal to diverse communities. For example, users interested in finding template structures for homology modeling should find particularly useful the visual analysis of structure similarity networks; the analysis of symmetric regions provides complementary information to infer functionality; the complete set of accurate pairwise structure alignments and their corresponding structure-based sequence alignments constitutes a robust benchmark for membrane pro- The results shown were additionally filtered by the term 'Chain A' and sorted by TM-score to reveal the least similar structures. (F) Structure-based sequence alignments of representative structural homologues relative to BetP allow interrogation of putative sodium binding sites or sodium-independent mechanisms. Alignments were extracted from pairwise structure alignments provided in EncoMPASS, and combined using pyali (https://github.com/christang/pyali). The positions of Na1 and Na2 site residues in BetP are indicated by arrows colored by repeat.


Nattokinase: A Promising Alternative in Prevention and Treatment of Cardiovascular Diseases

Cardiovascular disease (CVD) is the leading cause of death in the world and our approach to the control and management of CVD mortality is limited. Nattokinase (NK), the most active ingredient of natto, possesses a variety of favourable cardiovascular effects and the consumption of Natto has been linked to a reduction in CVD mortality. Recent research has demonstrated that NK has potent fibrinolytic activity, antihypertensive, anti-atherosclerotic, and lipid-lowering, antiplatelet, and neuroprotective effects. This review covers the major pharmacologic effects of NK with a focus on its clinical relevance to CVD. It outlines the advantages of NK and the outstanding issues pertaining to NK pharmacokinetics. Available evidence suggests that NK is a unique natural compound that possesses several key cardiovascular beneficial effects for patients with CVD and is therefore an ideal drug candidate for the prevention and treatment of CVD. Nattokinase is a promising alternative in the management of CVD.

Introduction
Cardiovascular diseases (CVDs) are the most prevalent cause of deaths worldwide. In 2015, the number of CVD-related deaths represented 31% of all deaths globally (www.who.int/cardiovas-cular_diseases). To date, there are limited approaches available for the control and/or management of CVD-related mortality. 1 Natto, a cheese-like food made of soybeans fermented with Bacillus subtilis, has been consumed as a traditional food in Asian countries for more than 2000 years. Natto consumption is believed to be a significant contributor to the longevity of the Japanese population. 2 Recent studies demonstrated that a high natto intake was associated with decreased risk of total CVD mortality and, in particular, a decreased risk of mortality from ischaemic heart diseases. 2 Before the 1980s, very little was known about the mechanism by which natto consumption led to overall cardiovascular health. In 1987, Sumi et al 3 discovered that natto contained a potent fibrinolytic enzyme called nattokinase (NK). Since then, a considerable amount of NK research has been performed on NK in Japan, Korea, China, and the United States, and these studies confirmed that NK, an alkaline protease of 275 amino acid residues with a molecular weight of approximately 28 kDa, 4 is the most active ingredient of natto and is responsible for many favourable effects on cardiovascular health. First, NK has potent fibrinolytic/ antithrombotic activity. [3][4][5][6] In addition, in both animal and human studies, NK also has an antihypertensive, 7,8 anti-atherosclerotic, 9,10 lipid-lowering, 9,11 antiplatelet/anticoagulant, 12 and neuroprotective actions. 13,14 All these pharmacologic actions of NK have relevance to the prevention and treatment of CVD. Indeed, NK supplementation has shown to enhance markers of fibrinolysis and anticoagulation and to decrease blood pressure (BP) and atherosclerosis in human subjects. 8,9,[15][16][17] The most unique feature of NK is that, as a single compound, it possesses multiple CVD preventative and alleviating pharmacologic effects (namely, antithrombotic, antihypertensive, anticoagulant, anti-atherosclerotic, and neuroprotective effects). There are no other drugs or drug candidates with multiple pharmacologic properties similar to NK. In addition, NK is a natural product that can be administered orally, has a proven safety profile, is economical to use, and provides distinct advantages over other pharmaceutical products. It therefore has the potential to be developed as a new-generation drug for the prevention, treatment, and long-term care of CVD. 7 The present review aims to concisely summarise the key pharmacologic effects and mechanisms of action of NK with a focus on their relevance to, and potential for, the prevention and treatment of CVD. The advantages of NK as an agent for the management of CVD will be outlined. Some remaining issues surrounding NK as a drug candidate for CVD will be critically reviewed and analysed.

Pharmacologic Actions of NK
The pharmacologic effects and mechanisms of action of NK will be reviewed and summarised under the following subheadings: fibrinolytic/antithrombotic effects, anti-atherosclerotic

Potent fibrinolytic/antithrombotic effects of NK
Although natto has been consumed in Asia for thousands of years, the fibrinolytic property of NK was only discovered in 1987. 3 Since this initial discovery, multiple laboratory and human studies have consistently reported the potent antithrombolytic action of NK. 18,19 A considerable amount of work has been performed to evaluate the thrombolytic effects of NK in vitro and in animal models. Using a rat model, Fujita et al 5 examined the effect of NK on chemically induced thrombi in the common carotid artery (CCA) and found NK to be 4 times more potent than plasmin in thrombus dissolution. At a concentration of 2836 FU, NK lysed 88% of thrombi within 6 hours, 20 and NK exhibited significant prophylactic antithrombotic effects in vivo. The efficacy of NK against thrombosis was confirmed in a carrageenan-induced tail thrombosis model. 21 The survival rate of mice with pulmonary thrombosis was increased by NK and the formation of thrombosis in mice was remarkably inhibited by NK, demonstrating significant antithrombotic effects. 22 In addition, in a model of rat experimental pulmonary thrombosis, oral administration of NK led to a decrease in thrombus count and plasma euglobulin lysis time (ELT), as well as an increase in tissue plasminogen activator (tPA), indicating that NK is capable of activating plasma fibrinolysis in vivo. 23 Omura et al 24 further found that a purified protein layer, NKCP, which mainly consisted of NK, had both fibrinolytic and antithrombotic effects, which they described as being similar to that of heparin. In rats, Natto treatment shortened ELT and significantly prolonged partial thromboplastin time compared with a nontreated rat group. 25 It is now known that NK not only degrades fibrin directly and effectively but also increases the release of tPA with a subsequent increase in the formation of plasmin. 5,17,[26][27][28] Plasminogen activator inhibitor 1 (PAI-1) is the primary inhibitor of tPA and regulates fibrinolytic activity in the fibrinolytic cascade. 29 In a study investigating the mechanism by which NK exerted its fibrinolytic effect, NK enhanced fibrinolysis through cleavage and inactivation of PAI-1. 4 In this study, NK was shown to cleave active recombinant prokaryotic PAI-1 into low-molecular-weight fragments as well as enhance tissue-type plasminogen activator-induced fibrin clot lysis. The enhanced fibrinolytic activity observed in the absence of PAI-1 appeared to be induced through direct fibrin dissolution by NK. 4 Nattokinase also enhanced the production of clot-dissolving agents such as urokinase through the conversion of prourokinase to urokinase. 5,19 Furthermore, NK was shown to be capable of blocking thromboxane formation resulting in an inhibition of platelet aggregation without producing the side effect of bleeding. 12 Thus, NK was found to be a potent antithrombotic agent, and, by reducing thrombus formation, was able to slow the progression of plaque formation and reverse evolving atherosclerotic lesions. 10 Data from human studies also strongly support NK or natto as a potent and promising fibrinolytic agent. In an early human trial, oral administration of NK was shown to produce a gradual enhancement of fibrinolytic activity in the plasma, as indicated by plasma ELT and the production of tPA. 6 Following administration of natto bacillus (100 mg/kg) to healthy adult volunteers, ELT was reduced and tPA activity was increased significantly (P < .05). 23 In an open-label, self-controlled clinical trial, Hsia et al 17 found that after 2 months of administration, fibrinogen, factor VII, and factor VIII levels were decreased significantly implying a promising cardiovascular benefit from NK administration. Even after a single dose of oral NK at 2000 FU, blood fibrin/fibrinogen degradation  product levels were significantly increased 4 hours after NK administration (P < .05), confirming elevation of thrombolysis and anticoagulation profiles. 16 Hsia's study supported NK as a useful fibrinolytic/anticoagulant agent to reduce the risk of thrombosis and CVD in humans.
Thrombosis is a common pathology underlying ischaemic heart disease, ischaemic stroke, and venous thromboembolism. 30 Thrombolytic/fibrinolytic therapy is used in myocardial infarction, cerebral infarction, and, on occasion, in massive pulmonary embolism. Fibrinolytic therapy can, in many circumstances, reduce mortality. As a result, the use of such treatments has increased rapidly in recent years. 31 Commonly used commercial thrombolytic agents, such as plasminogen activators, streptokinase, and anisoylated plasminogen streptokinase-activator complex, all come with serious drawbacks and adverse effects. These thrombolytic drugs are usually expensive, have a short half-life after intravenous administration, and lead to uncontrollable acceleration of fibrinolysis and haemorrhaging. To overcome these risks, a safer thrombolytic agent such as NK is needed for the treatment of thrombolytic disorders or CVD. 5

Anti-atherosclerotic and lipid-lowering effects of NK
The underlying pathological change shared by many CVDs is atherosclerosis, which is the primary cause of heart disease and stroke. 32 Therefore, drugs that have anti-atherosclerotic effects would have broad clinical relevance; however, these types of drugs are difficult to come by. Nattokinase would appear to be one of those drugs that have such promising anti-atherosclerotic and lipid-lowering effects.
A number of animal studies exploring NK as an anti-atherosclerotic agent have demonstrated that dietary natto extract supplementation suppresses intimal thickening in rats when compared with the control group. 10,33 The suppression of intimal thickening after vascular injury may be attributable to the enhanced thrombolytic activities of NK. 10,33 There was a striking difference in intimal thickening between the control group and natto fed group where the animals were fed the natto extract or NK for 3 weeks before the injury. Both NK and the natto extract suppressed intimal thickening in rats with endothelial injury. Although Chang et al 34 believed that the natto extract suppressed intimal thickening through a synergistic effect attributed to its antioxidant and antiapoptotic properties, other study demonstrated that NK prevented arteriosclerosis by its direct antioxidant effect leading to reduced lipid peroxidation and improved lipid metabolism (inhibition of low-density lipoprotein [LDL] oxidation). 35 When used in combination with red ginseng, NK was found to reduce aortic plaque area in hypercholesterol diet-fed rabbits. 36 In a recent clinical study conducted in our laboratory, 9 daily NK supplementation was an effective way to suppress the progression of atherosclerosis in patients with atherosclerotic plaques. Following NK treatment for 26 weeks, there was a significant reduction in CCA intima-media thickness (CCA-IMT) and carotid plaque size when compared with the baseline before treatment. The carotid plaque size and CCA-IMT reduced from 0.25 ± 0.12 cm 2 to 0.16 ± 0.10 cm 2 and from 1.13 ± 0.12 mm to 1.01 ± 0.11 mm, respectively. The reduction in the NK group was more significant (P < .01) than that in the group treated with simvastatin (daily dose of 20 mg). Our data suggested that NK was a better alternative to statins, a commonly used drug to reduce atherosclerosis, and furthermore, NK could be a viable alternative therapy for cardiovascular attack and stroke in patients. 9 The underlying mechanisms by which NK suppresses atherosclerosis are not known. Early studies indicated that NK enhanced thrombolytic activities. 10,33 Available data suggest that the anti-atherosclerotic effect of NK is due to the collective effects of the combination of antithrombotic, anticoagulant, antioxidant, and lipid-lowering effects of NK or natto extract containing NK. 9,[35][36][37] In addition to its anti-atherosclerotic effects, NK or natto extract also has a favourable effect on lipids. Using NK or natto extract containing NK, studies from various laboratories 25,[36][37][38][39][40][41] confirmed that NK has a hypolipidaemic effect and can significantly reduce the increased serum triglycerides, total cholesterol, and LDL cholesterol (LDL-C) levels in animal models. Our studies found that in patients with hyperlipidaemia, NK treatment (26 weeks at 6500 FU) reduced total cholesterol, LDL-C, and triglycerides. In addition, NK increased the level of high-density lipoprotein cholesterol (HDL-C). 9 However, in a small pilot study, Wu et al 42

Antihypertensive effects of NK
Natto has been observed to have a beneficial effect on CVDs by lowering BP. In an in vitro study, it was found that natto contained relatively strong inhibitors that suppress angiotensin-converting enzyme (ACE), a key enzyme responsible for the production of a hypertensive peptide hormone, angiotensin II, in the renin-angiotensin system. 43 In 2008, Kim et al 8 conducted the first randomised doubleblinded placebo-controlled trial to investigate the effects of NK supplementation on BP in prehypertension or stage 1 hypertension subjects. Oral administration of NK for 8 weeks resulted in a reduction of both systolic and diastolic BP (net changes were −5.55 and −2.84 mm Hg, respectively, P < .05). The results suggested that NK could play a role in the prevention and treatment of hypertension. 8 Very recently, a study by Jensen et al 15 suggested that NK consumption for 8 weeks was associated with beneficial changes to BP in patients having hypertension. This is consistent with several laboratory studies demonstrating effective reduction of BP in spontaneous hypertensive rats by NK administration. 7,[44][45][46][47] The mechanism by which NK decreased BP in hypertensive conditions is not clear. Jensen et al 15 found that the decline of BP in hypertensive patients in their study was independent of plasma renin activity and the role of ACE in the antihypertensive action of NK is somewhat controversial. Although studies in humans showed that ACE concentrations did not demonstrate a statistically significant difference in patients receiving NK treatment, 8 other studies using animal models suggested that the antihypertensive action of NK is related to the inhibition of angiotensin I-converting enzyme. 46 In a study using spontaneously hypertensive rats, NK and its fragments reduced hypertension by different mechanisms. 7 Nattokinase may decrease BP through cleavage of fibrinogen in plasma, whereas the fragments of NK potentially prevented the elevation of plasma angiotensin II level to suppress hypertension in rats. 7 Ibe et al 47 found that natto significantly decreased BP 4 hours after oral administration in spontaneously hypertensive rats and the ACE inhibition of natto was dose dependent. Natto extracts from B subtilis-fermented pigeon pea also significantly increased ACE inhibitor activity in hypertensive rats. 44 Given the significant side effects associated with the longterm use of antihypertensive drugs, 48 NK could be a promising alternative for the management of hypertension in patients with CVD. Thus, daily use of NK could be a successful strategy for the treatment of hypertension. 45

Antiplatelet/anticoagulant effects of NK
Low-dose aspirin (85-100 mg daily), as a potent anticoagulant agent, is widely used for the prevention of heart attacks, stroke, and atherothrombotic diseases. Aspirin exerts its antiplatelet action by inhibiting cyclooxygenase (COX) and subsequently reducing the synthesis of thrombogenic thromboxane A2 (TXA2) in platelets. However, the long-term use of aspirin comes with serious gastrointestinal (GI) side effects and bleeding. 49 In a study comparing the antiplatelet effects of NK and aspirin, NK was shown to display excellent antiplatelet aggregation and antithrombotic activities in vitro and in vivo, inhibiting thromboxane B2 formation from collagen-activated platelets. 12 Wang et al 50 found that NK decreased fibrinogen levels in a cerebral ischaemic model and concluded that this was mediated by a pathway similar to that of aspirin. Natto showed an excellent inhibitory effect on platelet aggregation induced by adenosine 5ʹ diphosphate and collagen. 25 In addition, NK was found to have positive in vitro haemorheological effects by decreasing red blood cell aggregation and low-shear viscosity. 51 All the above data suggest that NK could be a good candidate, without any obvious adverse effects, for the improvement of blood flow and possibly superior to aspirin.

Neuroprotective effects of NK
Nattokinase is capable of degrading amyloid fibrils at neutral pH and normal body temperature, suggesting a role in the treatment of amyloid-related diseases such as Alzheimer disease (AD). 52 The ability of NK to dissociate amyloid suggested that NK was a potential drug candidate for amyloid-related disorders and this was confirmed in a recent study involving both in vivo and in vitro models. 53 Oral administration of NK in the rat model of AD demonstrated a positive effect in modulation of specific factors in the AD pathway. 13 In a rat model of cognitive deficits of AD induced by intoxification of colchicine, nano-nutraceuticals containing NK were demonstrated to enhance the impaired learning and memory capability and to be effective inhibitors in the suppression of amyloid-β and BACE-1 activity, thus suggesting a neuroprotective efficacy of NK. 54 Ahmed et al 55 demonstrated that NK, at a dose of 360 FU/kg, significantly decreased cholinesterase activity, TGF-β, IL-6, and p53 levels accompanied by a significant increase in Bcl-2 levels as compared with an untreated AD control group. Their data suggested that the neuroprotective effect of NK was due to its proteolytic, anti-inflammatory, and antiapoptotic effects.
In one of the early clinical studies describing the role of NK in the prevention of stroke progression in patients with acute ischaemic stroke, the authors demonstrated a clear neuroprotective effect in patients. 56 The neuroprotective effects of NK were demonstrated in a photothrombotic stroke mouse model. Ahn et al 57 suggested that the neuroprotective effect in the ischaemic brain was induced through improved blood flow by inhibiting platelet aggregation and thrombus formation by NK, as reported in a conference abstract. Oral intake of NK was also shown to confer neuroprotective effects against focal cerebral ischaemia as evidenced by reduced infarct volume in gerbils through enhanced fibrinolytic activity. 50 Regarding the mechanism by which NK protects the brain in ischaemic stroke, a study by Ji et al 14 showed that the neuroprotective effect of NK was associated with its antiplatelet activity, antiapoptotic effect, its ability to relax vascular smooth muscle, and its protection of endothelial cells through increased fibrinolytic activity and facilitating spontaneous thrombolysis.
Cardiovascular health is closely linked to brain health. There is a strong causal association between CVD risk factors with the incidence of cognitive decline and AD. 58 Therefore, neuroprotective effects of drugs used in patients with CVD have great potential and would add benefits to patients' overall outcome. Despite more than 100 neuroprotective agents demonstrating neuroprotective effects on focal ischaemic stroke in recent preclinical studies, none has proven to be beneficial in clinical studies. 59 Hence, more effective drugs, with proven Chen et al 5 neuroprotective actions, for the management of cerebrovascular diseases, such as ischaemic heart disease and stroke, are urgently needed.

Clinical Human Studies of NK
Although NK has recently gained popularity as a candidate drug for CVD, clinical investigations of NK in humans are relatively limited. The first NK clinical study conducted in 1990 aimed to assess its fibrinolytic activity in healthy subjects after oral administration. 6 A further Japanese study examined the efficacy of NK in the prevention of stroke progression in patients with acute ischaemic stroke and showed clear beneficial effects in these patients after oral administration. 56 As shown in Table   1, there are no more than 10 published clinical studies covering the use of NK or the natto extract. Other related NK clinical studies include investigations of antihypertensive effects in patients in North America 15 and Asia. 8 Further studies also include anti-atherosclerotic and lipid-lowering effects, 9,42 fibrinolytic activity and effects on coagulation, 6,16 and pharmacokinetics 60 and toxicology 61 in humans. The low number of NK clinical studies is probably related to the fact that NK is not registered as a drug, but, as a nutritional supplement, and, to date, clinical evidence has not been adequately conclusive.
It should be pointed out that clinical.gov has, to date, published around 10 clinical trials of NK or natto with the most recent one, entitled "Nattokinase Atherothrombotic Prevention Abbreviations: BP, blood pressure; CCA-IMT, common carotid artery; Cvd, cardiovascular disease; HdL-C, high-density lipoprotein cholesterol; LdL-C, low-density lipoprotein cholesterol; NK, nattokinase.
Study" (NCT02080520) being currently conducted by a team of researchers at the University of Southern California (USC; details can be found in https://clinicaltrials.gov). Among these studies, the one from USC (NCT02080520) is a considerably large-sized phase 2 clinical trial and we are awaiting the outcomes of this trial. More human studies are required to ascertain the efficacy of NK, the indications for its use, and possible drug interactions when used concurrently with other drugs used to treat patients with CVD.

Advantages of NK
As stated above, the advantages of NK include a proven safety profile with a long history of human consumption, the convenience of oral administration (many antithrombotic drugs are injectable), ease of mass production, and it possesses multiple key favourable cardiovascular effects.
The fact that multiple clinical trials on NK have been conducted, or approved to be conducted, provides a good indication of the safety of administration in human subjects. Recent toxicology studies (both in vivo and in vitro) have provided strong corroboration of the safety of oral consumption of NK and these studies concluded that, to-date, there is no toxicologic concern for NK human consumption. 61,62 The safety of NK has been comprehensively evaluated in several Good Laboratory Practice (GLP)-compliant studies in rodents and in human volunteers. 61 Nattokinase was found to be nonmutagenic and nonclastogenic in vitro, and no adverse effects were observed in 28-day and 90-day subchronic toxicity studies conducted in rats at doses up to 1000 mg/kg/day, which is 100 times higher than the usual dose used in humans. In human volunteers, no adverse effects were observed following 4 weeks of NK consumption at a daily dose of 10 mg/kg for 28 days. 61 No sign of toxicity was observed in the Ames test, the bone marrow cell micronucleus tests, or the mouse sperm abnormality tests following NK treatment. In an acute toxicology study conducted by Fu et al, 62 no sign of toxicity was observed at daily doses of NK of up to 4 g/kg for a period of 30 days. All these data collectively confirm NK to be a safe product with very little or no toxicologic concerns.
In Japan, China, Korea, and Western countries including Australia, United States, Canada, and Europe, NK products are widely used in the form of tablets or capsules as nutritional supplements for health promotion including blood thinning, blood clot prevention, and improved circulation. The cost of NK supplements is relatively inexpensive, with the cost of daily doses of 200 to 400 mg being approximately US $0.3 to US $0.8 depending on the source of the product. With the advancement of NK production and purification technologies, it is expected that the cost of NK will continue to be low, making it an economically viable drug. 28 One of the great advantages of NK is that, as an orally active compound, it possesses many key favourable cardiovascular effects for patients with CVD. Nattokinase has the potential to streamline drug use and to improve patient compliance in CVD prevention and treatment.

Remaining Issues
One key question that has not been adequately addressed is the mechanism by which NK is absorbed into the bloodstream after oral administration. There is no current convincing data available to demonstrate the bioavailability and metabolism of NK administered as an oral dose. The only related study is one showing that, in rats, intraduodenal administration of NK, at a dose of 80 mg/kg leads to degradation of fibrinogen in plasma. 63 It should be noted that the NK used was not subjected to digestive gastric fluids and the dose was high in pharmacologic terms. Thus, the relevance of the data is questionable. Recent pilot studies on the pharmacokinetics of NK have not provided an answer to the question of how NK is absorbed into the body. 16,60 In fact, the pharmacokinetic data are rather inconsistent and mismatched with pharmacodynamic activities. Ero et al 60 claimed that NK was measured directly in their study with a peak concentration at 13.3 hours after oral administration, whereas the thrombolytic activity of NK, measured in another study, was shown to peak 2 to 4 hours after oral dosing. 16 These data imply that measurements of NK in the published studies may not be those of the intact NK molecule. Although it has been assumed that NK is stable in the GI tract 19,28 , it is important to address the key question of the mechanism of NK absorption into the body, with more convincing direct evidence. In fact, several groups have pointed out that the molecular size of NK is considered to be generally too large for oral absorption through the GI tract. 19 It is recognised that NK may be susceptible to chemical oxidation and subsequent inactivation, or denaturation, in the GI tract. 64 Thus, further studies are required to fully understand the pharmacokinetics of NK.
Potential drug interactions and contraindications to the use of NK in humans are not currently known. There is one report of a patient concurrently using aspirin and NK (400 mg daily), experiencing an acute cerebellar haemorrhage. 65 In addition, multiple microbleeds were demonstrated on brain magnetic resonance images suggesting that NK may increase risk of intracerebral haemorrhage in patients who have bleedingprone cerebral microangiopathy and are taking aspirin concurrently. 65 In another report, a patient developed a thrombus in a mechanical valve after nearly a year of NK use without warfarin and underwent a successful repeat valve replacement. 66 Use of multiple drugs concurrently is common in patients with CVD. Therefore, it is particularly important to elucidate the potential possible interactions of NK with other drugs acting on the cardiovascular system.
Given the fact that serine proteases, currently in clinical use for stroke patients, have been implicated in tumour growth and metastasis, 67,68 it is necessary to design studies to evaluate the potential carcinogenic effects, and effects on cell invasiveness, of NK use in in vitro and in vivo model systems, Chen et al 7 and to ascertain the presence of any cardiotoxicity associated with this compound.

Concluding Remarks
In recent years, NK has gained popularity due to its beneficial cardiovascular health effects. Research defining the pharmacologic effects of NK and refining the production and purification of NK is increasing. It is particularly exciting to observe that every year contemporary research on NK continues to find more positive and favourable clinical benefits resulting from NK use. For example, NK was recently demonstrated to effectively reduce nasal polyp tissue through fibrin degradation, 69 indicating its potential in the treatment of chronic rhinosinusitis. The emerging picture of NK is one of the unique natural compounds with many favourable health effects, especially for patients with CVD. Importantly, NK simultaneously effects several key favourable benefits for thrombosis, hypertension, atherosclerosis, hyperlipidaemia, platelet aggregation, and neuroprotection in patients with CVD. These multiple benefits to patients make the role of NK very unique in the prevention and treatment of CVD with reduced side effects commonly associated with conventional CVD drugs. In addition to the unique multiple functions of NK, other advantages to the patient include ease of oral administration, a robust safety profile and a long in vivo half-life.
In summary, compared with traditional antithrombotic and antihypertensive drugs, NK is characterised by high safety, low cost, simple production process, oral availability, and long in vivo half-life. As such, it is expected to become a new-generation drug for thrombotic disorders or CVDs. Although human trials and clinical studies demonstrate the clinical benefits of NK in the clinical settings, there remain a number of limitations. However, all the available data are encouraging and promising and further clinical trials are needed to fully examine the prospect of NK as an alternative medication to tPA, aspirin, warfarin, or newer anticoagulants in the management of CVD. In the near future, it is possible that patients with CVD may need only a single NK pill to replace multiple drugs administered for the prevention and management of CVD, including tPA, antihypertensives, statins, aspirin, and warfarin.
There are a few obstacles that remain to be overcome for NK to become more widely accepted as a drug, or drug candidate, for CVDs and these include establishment of pharmacokinetic evidence of the absorption and metabolism of NK in humans and the requirement for extensive studies to elucidate potential drug interactions between NK and other cardiovascular drugs, which are commonly, and concurrently, used in the prevention, treatment, and management of CVD in affected patients.

Author Contributions
YL conceived and designed the content and structure of the review; HC, EM, and YL drafted, revised, and edited the review; NR, SL, NN, FS, and XQ revised and edited the manuscript.

Disclosures and Ethics
The authors have read and confirmed their agreement with the ICMJE authorship and conflict of interest criteria. The authors have also confirmed that this article is unique and not under consideration or published in any other publication, and that they have permission from rights holders to reproduce any copyrighted material.