Systems analysis-based assessment of post-treatment adverse events in lymphatic filariasis

Background Lymphatic filariasis (LF) is a neglected tropical disease, and the Global Program to Eliminate LF delivers mass drug administration (MDA) to 500 million people every year. Adverse events (AEs) are common after LF treatment. Methodology/Principal findings To better understand the pathogenesis of AEs, we studied LF-patients from a treatment trial. Plasma levels of many filarial antigens increased post-treatment in individuals with AEs, and this is consistent with parasite death. Circulating immune complexes were not elevated in these participants, and the classical complement cascade was not activated. Multiple cytokines increased after treatment in persons with AEs. A transcriptomic analysis was performed for nine individuals with moderate systemic AEs and nine matched controls. Differential gene expression analysis identified a significant transcriptional signature associated with post-treatment AEs; 744 genes were upregulated. The transcriptional signature was enriched for TLR and NF-κB signaling. Increased expression of seven out of the top eight genes upregulated in persons with AEs were validated by qRT-PCR, including TLR2. Conclusions/Significance This is the first global study of changes in gene expression associated with AEs after treatment of lymphatic filariasis. Changes in cytokines were consistent with prior studies and with the RNAseq data. These results suggest that Wolbachia lipoprotein is involved in AE development, because it activates TLR2-TLR6 and downstream NF-κB. Additionally, LPS Binding Protein (LBP, which shuttles lipoproteins to TLR2) increased post-treatment in individuals with AEs. Improved understanding of the pathogenesis of AEs may lead to improved management, increased MDA compliance, and accelerated LF elimination.


Introduction
Lymphatic filariasis (LF) is a disabling neglected tropical disease that is caused by the mosquito-borne filarial parasites Wuchereria bancrofti, Brugia malayi and B. timori. Adult worms live in the human host's lymphatic system and release larval parasites (microfilariae or Mf) that circulate in the blood. Infection and host inflammatory responses to the parasite can lead to severe morbidity including lymphedema, hydrocele and elephantiasis [1]. To combat this disease the WHO launched the Global Program to Eliminate Lymphatic Filariasis (GPELF) in the year 2000 with the goal of eliminating LF as a public health problem by 2020. The program uses mass drug administration (MDA), to cure infections, prevent disease, and reduce transmission of new infections. As of 2016 a total of 6.7 billion treatments had been delivered to more than 850 million individuals [2], making GPELF the largest public health intervention for an infectious disease to date based on MDA. Drugs used for LF MDA include albendazole (ALB), ivermectin (IVM) and diethylcarbamazine (DEC). MDA with two-drug combinations is usually provided annually for 4-6 years. The combinations used are ALB with IVM in sub-Saharan Africa and ALB with DEC in other regions [1]. New studies have shown that combining all three drugs increases the anti-filarial effect and potentially decreases the number of required treatment rounds [3][4][5][6][7]. This new triple therapy (IDA) was recently recommended by the WHO as the preferred regimen for LF elimination in some settings [8].
Although LF treatment is safe, transient mild to moderate systemic adverse events (AEs) are common following treatment, and these are especially common in individuals with Written informed consent was obtained from all participants. Adults with W. bancrofti microfilaremia were randomly assigned to one of four treatment arms (all oral medications): the standard LF treatment regimen for Côte d'Ivoire (200μg/kg IVM plus 400mg ALB), IDA: 200μg/kg IVM plus 6mg/kg DEC and 400mg ALB, a single 400 mg dose of ALB, or a single 800 mg dose of ALB. A subset of ninety-five individuals treated with either IVM/ALB, IDA or 400mg ALB had samples processed for use in the AE study described in this paper. We selected these individuals based on the availability of pre-and post-treatment samples and clinical AE data. Metadata of these 95 individuals is shown in S1 Table. A physical examination was performed shortly before treatment, and vital signs were recorded. A review of systems (ROS) questionnaire was also completed to assess subjective symptoms prior to treatment. Venous blood (3 to 4 mL in EDTA) was collected immediately before participants received treatment. Participants were interviewed and examined the next day to assess AEs, and venous blood was collected approximately 24 hours after treatment. Blood samples were centrifuged within an hour of collection, and plasma was removed. The buffy coat (approximately 500μL) was carefully aspirated with a pipette and added to 1.8mL of RNAlater (Ambion, Foster City, CA). The plasma samples and buffy coat/RNAlater samples were stored at the study site at -20˚C, shipped frozen, and later stored at -80˚C.

Adverse event classification
AEs were categorized as mild, or moderate. Those with moderate AEs (n = 9) had at least two new or worsening subjective symptoms plus one objectively measured change in their vital signs (an increase in axillary temperature of � 0.8˚C to at least 37.4˚C post-treatment and/or a decrease in sitting systolic blood pressure of at least 20 mm Hg). Individuals with subjective or objective AEs that did not fulfill the criteria for moderate AEs were considered to have had mild AEs (n = 24). Individuals with no new objective or subjective symptoms after treatment were considered to have no AEs (n = 62).

Circulating filarial antigen assay
A direct sandwich enzyme immunoassay (EIA) was performed as previously described [13]. This assay uses the monoclonal antibody AD12 that binds to a carbohydrate epitope on circulating filarial antigen (CFA). It is important to note that the carbohydrate epitope recognized by AD12 is present in many filarial glycoproteins [26]. However, the high molecular weight CFA is the only filarial antigen that is frequently detected in the blood of W. bancrofti-infected individuals. Pre-and post-treatment plasma samples from 95 individuals were tested in duplicate. The detection range of the CFA EIA was 6.3 to 400 ng/mL. CFA was detected in all samples, but two individuals had extremely high CFA levels that were above the upper detection limit of the assay. These samples were retested after dilution to obtain baseline CFA concentrations. The percent change in CFA relative to baseline following treatment was calculated for each participant. Sample pairs with pre-treatment values less than 20 ng/mL (7 individuals) were excluded from the percent calculations, because they were near the lower detection limit of the assay. Kruskal-Wallis H tests were used to compare percent change values and absolute values between the three AE groups and the three treatment arms. Wilcoxon signed-rank tests were used to compare pre-and post-treatment CFA levels within the three AE groups.
as AD12 was directly conjugated to 2mL of agarose Affigel 10 beads (Bio-Rad, Hercules, CA) according to the manufacturer's protocol. Conjugated beads were stored as a 50% solution in PBS. 40μL of conjugated beads were mixed with 50μl of human plasma and 300μl PBS and rocked overnight at 4˚C. The beads were washed four times with cold PBS and then boiled in 1X NuPAGE LDS sample buffer (Invitrogen, Carlsbad, CA) to release bound antigens. Proteins were resolved by SDS-PAGE using a 4-12% bis-tris NuPAGE gradient gel (Invitrogen) and transferred to 0.45μM nitrocellulose membrane (Amersham, Piscataway, NJ). Membranes were blocked with 5% milk in phosphate buffered saline with tween-20 (PBS-T) followed by incubation with a peroxidase-conjugated AD12 antibody (1:3000 dilution) for one hour at room temperature. After washing, membranes were incubated with Clarity Western ECL substrate (Bio-Rad). Chemiluminescence was detected by a ChemiDoc imager (Bio-Rad), and results were analyzed using Image Lab 5.2.1 software.

Immune complex assay
CIC were measured with a C1q ELISA. C1q was purchased (Sigma-Aldrich, St. Louis, MO), and a previously published protocol was followed [13]. Plasma samples were available from 41 individuals for this assay (8 with moderate AEs, and 33 with no AEs), and both pre-and posttreatment samples were tested in duplicate. Negative control samples (plasma samples from healthy North American subjects and deionized water) were tested on each plate. Values were expressed as μg/mL of AHG (aggregated human gamma globulin) (Invitrogen). The range of detection for the CIC ELISA was 0.0006 to 6 μg/mL of AHG, and all samples had detectable CIC. Mann-Whitney U tests were used to compare absolute CIC levels between the two AE groups pre-and post-treatment. The Wilcoxon signed-rank test was used to compare posttreatment CIC levels to baseline levels within AE groups. The Kruskal-Wallis H test was used to compare absolute CIC levels between the three treatment arms post-treatment.

Complement component assays
Nine individuals with moderate AEs were matched to individuals with no AEs following treatment. Matching was based on age, sex, baseline Mf count, and treatment arm (S2 Table). Complement component 3 (C3), complement component 4 (C4) and Factor B (FB) were measured in the 36 samples (18 matched case-control subjects pre-and post-treatment) with ELISA kits (AssayPro, St. Charles, MO). The C3 and C4 assays were competitive enzyme immunoassays, and the FB assay was a sandwich ELISA. Each sample was tested in duplicate and manufacturer's protocol was followed. Paired t-tests were used to compare pre-and post-treatment complement component levels by AE group.

LPS binding protein assay
LBP was measured with a sandwich ELISA kit (Abnova, Taipei, Taiwan). Plasma samples from the same 18 matched case-control subjects were included. Each sample was tested in duplicate and manufacturer's protocol was followed. Paired t-tests were used to compare pre-and posttreatment levels within both AE groups. The range of detection for the LBP ELISA was 5 to 50 ng/mL.

Cytokine assays
Twenty-seven cytokines (IL-1β, IL-1RA, IL-2, IL-4, IL-5, IL-6, IL-7, IL-8, IL-9, IL-10, IL-12  (p70), IL-13, IL-15, IL-17, basic FGF, eotaxin-1, G-CSF, GM-CSF, IFN-γ, IP-10, MCP-1, MIP-1α, MIP-1β, PDGF-BB, RANTES, TNF-α, and VEGF) were measured with the MAGPIX system with the Bio-Plex Human 27-Plex Cytokine Panel and Bio-Plex Cytokine Reagent Kit (Bio-Rad). Plasma samples from the same 18 matched-control subjects were included. A previous paper includes the detailed protocol [13]. Briefly, all samples were tested in duplicate, standard curves were calculated using the manufacturer's software, and the analysis considered mean concentrations (pg/mL) from two duplicate wells. Wilcoxon signed-rank tests were used to determine whether cytokine levels changed after treatment in either of the two AE groups. Mann-Whitney U tests were used to determine whether pre-or post-treatment cytokine levels were different between the two AE groups. For graphing, fold changes were calculated for each cytokine by AE groups, and samples with cytokines below the detection limit were assigned a value equal to half the value of the lowest pre-treatment sample concentration measured for that cytokine.

RNA preparation
RNA was extracted from pre-and post-treatment buffy coat samples from the same 18 matched case-control subjects. Total RNA was extracted using Qiagen RNeasy kits (Qiagen, Hilden, Germany) according to the manufacturer's protocol with an added homogenization step and on-column DNase digestion as follows. For each sample 200μL of the buffy coat/ RNAlater mixture was added to 700μL of the kit's RLT buffer and vortexed. The mixture was then added to a QIAShredder column (Qiagen) and centrifuged for 2 min at 16,000g. The flow-through was added to 700μL of 70% ethanol, and this mixture was added to a RNeasy column. Bound RNA was eluted in 30μL RNase-free water and stored at -80˚C. The quality and quantity of RNA was verified with a Bioanalyzer 2100 (Agilent Technologies, Cedar Creek, Texas). Samples were processed with the TruSeq Stranded Total RNA LT Sample Prep Kit with Ribo-Depletion using the manufacturer's protocol (Illumina, San Diego, CA). The RNA was high quality (average RIN value 9.3, range 8.5-10).

RNA sequencing and mapping
The 36 samples were sequenced in two batches. The first 14 samples were sequenced with the HiSeq2000 (2x 100 PE run and Illumina TruSeq Stranded Total RNA) platform, and the remaining 22 samples were sequenced with HiSeq4000 (2x 150 PE run and Illumina TruSeq Stranded Total RNA). Between 28-41 million read fragments per sample were mapped to 19,864 protein-coding genes. Raw reads were mapped to protein coding genes using HISAT2 (version 2.0.5) [27], and the human reference genome GRCh38.84. FeatureCounts [28] was used to count reads per gene.

Differential gene expression and overall expression patters
DESeq2 [29] was used to generate normalized read counts and to identify differentially expressed genes between the different comparator groups (namely AEs vs. no AEs and pre-vs. post-treatment). The program "R" with the biocLite package "DESeq2" was used. Gene expression results from individuals before and after treatment were considered to be repeated measures for the analysis. Principal component analysis (PCA) was performed for 500 genes with the greatest variability in expression (based on DESeq2 output, default settings), and distance metrics statistics [30] were used to determine whether groupings affected overall expression patterns. A clustering dendrogram (Euclidean distance, complete linkage) was also used to illustrate overall expression patterns, and this method considered all genes. A two-tailed binomial distribution with unequal variance (for categorical data), and Mann-Whitney U tests (for continuous variables) were used to identify over-represented metadata variables in the different clustering groups such as baseline Mf/mL and treatment group.

Functional enrichment in the post-treatment AE group
The online tool WebGestalt [31] was used to identify enriched KEGG pathways within genes that were upregulated post-treatment during AEs. The reference set was a list of all 19,864 genes with expression signals in the RNAseq data, and the default values were used except the significance level (FDR < 0.05). The program i-cisTarget [32] was used to identify enriched transcription factor binding sites in the upregulated gene set using default settings and database version 4.0.

Identification of similar expression profiles
GeneQuery is an online tool that can search the PubMed GEO database and compare transcriptional signatures to published gene expression profiles [33]. The input for the post-treatment AE profile were 744 genes that were identified by differential gene expression (DESeq2) to be upregulated post-treatment in individuals with moderate AEs.

Changes in peripheral blood leukocyte populations after treatment
CIBERSORT [34] is an analytical tool that can estimate the abundances of 22 leucocyte subtypes based on RNA-seq data. Pre-and post-treatment DESeq2 normalized read counts were used as input for the program. The standard LM22 (22 immune cell types) was the signature gene file, and all default settings were used. Thirteen cell subtypes had very low representation in this dataset (totaling less than 4% in all 36 samples), so the analysis was limited to the remaining subtypes. Percent change for each cell type post-treatment was calculated for the two AE groups, and Mann-Whitney U tests were used to assess the significance of differences by AE group.

Prioritization of the genes with altered regulation post-treatment in individuals with AEs
Random forest (RF) analysis was used to prioritize the 678 genes that were differentially expressed in individuals with moderate post-treatment AEs. RF was performed using the "R" package "randomForest" with 1000 trees and default values to analyse DESeq2 normalized read counts. Differentially expressed genes were ranked based on decreasing Mean Decrease in Accuracy values for 10 separate RF models. The Mean Decrease in Accuracy is the decrease in model accuracy from permuting the values in each feature. This metric is used to compare the impact of the variables in the model, and a large positive value indicates that a variable was closely linked to AE group across the dataset.

Preparation of cDNA for validation of selected differentially expressed genes
Additional RNA was extracted from residual buffy coat samples that were available for 34 of the 36 samples that were subjected to expression profiling (17 individuals pre-and post-treatment) as described above. Extracted RNA was treated with DNase I (Invitrogen), and RNA was measured with a NanoDrop 1000 Spectrophotometer (Thermo Scientific, Waltham, MA). Each sample was diluted to approximately 0.5ng/μL RNA with RNase-free water. cDNA was prepared with SuperScript II Reverse Transcriptase (Invitrogen) and with Oligo(dT) [12][13][14][15][16][17][18] according to the manufacturer's protocol .

Validation of the top differentially expressed genes by quantitative reversetranscription PCR (qRT-PCR)
SYBR Green based assays were performed for the top eight genes based on the RF analysis (DIP2B, ZCCHC6, RBPJ, PELI1, FNDC3B, TLR2, LTBR, NT5C2) that were upregulated in peripheral blood leukocytes (PBL) after treatment in participants who experienced moderate AEs. Four housekeeping genes (SDHA, ACTB, HPRT1 and YWHAZ) were used as controls for these experiments. We chose these based on prior validation as housekeeping genes by others [35] and because our results confirmed their stable gene expression before and after treatment. Pre-validated primer sets for the eight target genes were purchased from KiCqStart SYBR Green Primers (Sigma-Aldrich), and primers for the four housekeeping genes were made using previously published sequences (IDT, Coralville, IA) (S3 Table). Real-time PCR reactions were performed with 10μL of SYBR Green Master Mix (Applied Biosystems, Foster City, CA), 450 nmol/L of each primer, and 2μL cDNA (approx. 1ng RNA) with a final volume of 20μL. Thermal cycling was performed for 40 cycles with a QuantStudio 7-Plex Real-Time PCR System (Applied Biosystems), and cycle threshold (Ct) values were determined using the manufacturer's software. All samples were tested in duplicate, and each plate included a negative water control and a RNA sample that had not been treated with reverse transcriptase. Delta delta Ct values were calculated [36], using the geometric mean Ct value of three housekeeping genes (SDHA, ACTB and YWHAZ) as a normalization factor [35]. Student's t-tests were performed to compare baseline and post-treatment delta Ct values by AE group.

Statistical methods
All statistical analyses were performed with IBM SPSS (version 23). Shapiro-Wilk tests were used to test for normality in each sample set, and additional tests were performed as described in each section above. Logistic regression analysis was performed with the binary dependent variable AEs (moderate AEs vs. no AEs). The independent variables considered included age, sex, treatment arm, baseline Mf/mL, and baseline CFA level.
Separate RF analyses of gene expression and plasma biomarker data were performed 10 times using 1000 trees. The output was the average Mean Decrease in Accuracy over the 10 runs for each variable.

Ethical review
Institutional review boards in Cleveland, USA (University Hospitals Cleveland Medical Center IRB #08-14-13) and in Côte d'Ivoire (Comité National d'Ethique et de la Recherche, CNER, N: 008/MSLS/CNER/-kp) approved the clinical trial study protocol. Written informed consent was received from all participants prior to inclusion in the study.

Study population
This study of the pathogenesis of AEs that occur after LF treatment used human samples that were obtained as part of a clinical trial for LF that was conducted in Côte d'Ivoire. Full results from that study have not yet been published, but early results from the study have been reported in published abstracts [4,5]. Briefly, 189 W. bancrofti-infected adults were randomly assigned to one of four treatment arms as described above, and all participants had AE assessments performed 24 hr after treatment. The AE study enrolled a subset of 95 treated participants. S1 Table summarizes the specific analyses that were performed on samples from each of the 95 individuals. Nine of these participants experienced moderate AEs (Table 1), 24 had mild AEs, and 62 had no AEs. There was no difference in age or sex distribution between the three AE groups (S4 Table).

Multiple additional filarial antigens were detected in post-treatment plasma samples
Baseline CFA levels were positively correlated with baseline Mf counts (Spearman's rho: 0.51, P < 0.001), and absolute CFA levels were significantly higher at baseline in individuals who developed moderate AEs compared to individuals who developed mild or no AEs after treatment (P = 0.012 by Kruskal-Wallis H test). Plasma CFA levels increased post-treatment in all three AE groups, but the increases were greater in persons with moderate AEs (P < 0.05 by Kruskal-Wallis H test) (Fig 1). Percent changes in CFA levels post-treatment were significantly lower in the individuals treated with only ALB compared to those in individuals treated with IVM/ALB or IDA (P < 0.0001 by Kruskal-Wallis H test).
Western blot analysis was performed for nine pre-and post-treatment plasma pairs to compare CFA patterns detected in plasma before and after treatment. All nine pre-treatment plasma samples contained only a single high molecular weight parasite antigen as expected, and this antigen was also present in post-treatment plasma samples. However, four of the posttreatment plasma samples contained many parasite antigens that were not present before treatment. Two examples are shown in S1 Fig (P1 had moderate AEs, and P2 had mild AEs), and this pattern was also observed in plasma from two other participants who experienced moderate AEs following treatment. Western blot results obtained with five other post-treatment plasma samples tested (3 from persons with moderate AEs, and 2 from persons with no AEs) were no different from those observed in pre-treatment samples.

CIC did not increase post-treatment and the classical complement pathway was not activated in individuals with moderate AEs
All pre-and post-treatment samples contained CIC. However, there was no difference in CIC levels between the two AE groups before or after treatment, and CIC levels did not significantly change after treatment in either AE group (S2 Fig). There was also no difference in post-treatment CIC levels between the treatment arms.
C3 levels significantly decreased post-treatment in individual with moderate AEs, but this change was not observed in individuals with no AEs (Fig 2A). C4 levels did not change in either AE group (Fig 2B). Factor B levels decreased post-treatment in most individuals with moderate AEs, but three individuals had increases in FB levels, and the group differences were not significant (Fig 2C).

LPS binding protein levels increased post-treatment in individuals with moderate AEs
LBP was detected in all pre-and post-treatment samples. LBP levels increased post-treatment in individuals with moderate AEs (P = 0.0007 by paired t-test), but they did not increase in individuals with no AEs (Fig 2D).

Many plasma cytokines increased in plasma after treatment in persons who experienced moderate AEs
Plasma cytokine levels before and after treatment are shown by AE group in Fig 3. Seven cytokines (IL-8, MCP-1, VEGF, TNF-α, MIP-1β, G-CSF and IFN-γ) increased post-treatment only in individuals who experienced moderate AEs (P < 0.05 by Wilcoxon signed-rank test). Five cytokines (IL-6, IL-10, IL-1RA, IP-10 and MIP-1α) increased post-treatment in individuals with and without AEs (P < 0.05 by Wilcoxon signed-rank test), but three of these (IL-6, IL-10 and IL-1RA) had significantly higher levels post-treatment in individuals with moderate AEs compared to individuals with no AEs (P < 0.05 by Mann-Whitney U tests) (Fig 3). The remaining 15 cytokines (IL-1β, IL-2, IL-4, IL-5, IL-7, IL-9, IL-12 (p70), IL-13, IL-15, IL-17, basic FGF, eotaxin-1, GM-CSF, PDGF-BB and RANTES) did not change after treatment in either AE group. There was no difference in pre-treatment cytokine levels between individuals that would develop moderate AEs and individuals that would not develop AEs for any of the 27 cytokines.

Changes in gene expression associated with moderate post-treatment AEs
Raw and processed RNA-seq data are available to the public on NCBI's Gene Expression Omnibus (Accession number: GSE110146).
We analyzed changes in gene expression in PBL after treatment to further elucidate host responses associated with AEs. Post-treatment gene expression profiles from persons who developed moderate AEs clustered together using a clustering dendrogram based on gene expression profiles across all genes (Fig 4A). Post-treatment AE samples were significantly overrepresented in the fourth group (bolded in Fig 4A, P-value < 0.0001 for enrichment within the cluster, two-tailed binomial distribution with unequal variance). Higher levels of baseline Mf/mL were also observed in this group (P-value = 0.038, Mann-Whitney U test), but none of the other metadata categories (treatment arm or village) were over-represented. Age also did not affect clustering. A similar pattern was observed by principal components analysis (Fig 4B), where post-treatment moderate AE samples clustered together and were clearly separated from their pre-treatment controls (P-value = 0.005 by PERMANOVA [30]). No other differences were significant between the four groups by PCA.
We used differential gene expression analysis to identify the genes that were responsible for the clustering of the post-treatment moderate AE samples. At a very stringent significance threshold (P < 10 −5 according to DESeq2 output), 783 genes were identified to be upregulated  (2), fever (2), rash(1), joint pain (1) Post-treatment adverse events in lymphatic filariasis after treatment (n = 744) or before treatment (n = 39) in individuals who experienced moderate AEs (Fig 4C). No differences were observed pre-or post-treatment in individuals with no AEs when this stringent significance threshold was used. However, at a less stringent P-value of 0.05, there were 126 genes upregulated post-treatment and 19 genes upregulated pre-treatment in individuals without AEs (Fig 4D). There was only one overlapping gene in the genes upregulated pre-treatment in individuals with and without AEs, whereas the majority of the genes upregulated post-treatment in individuals with no AEs were also upregulated post-treatment in individuals with AEs ( Fig 4E).
We then assessed whether there was evidence for functional enrichment in the genes upregulated post-treatment in individuals with AEs. Among the 744 upregulated genes post-treatment in the AE samples a total of 35 enriched biological pathways (KEGG) were identified (S5 Table), and these included TLR signaling and downstream pathways such as NF-κB, TNF and Post-treatment adverse events in lymphatic filariasis Jak/STAT. Many individual genes in the TLR signaling pathway, including TLR2, TLR6, STAT1 and STAT2, were identified by DESeq2 to be significantly upregulated post-treatment in individuals with AEs. A separate analysis (i-cisTarget) predicted that six transcription factors were over-represented in the differentially expressed genes (S6 Table), and three of these, Post-treatment adverse events in lymphatic filariasis STAT1, STAT2 and IRF1, are downstream of TLR signaling. Convincingly, STAT1 and STAT2 were therefore identified by two independent analyses, signifying the importance of these two transcription factors in the development of AEs. IRF1 is activated by IFN-γ and is a major Post-treatment adverse events in lymphatic filariasis transcription factor of IL-8, and correspondingly both IFN-γ and IL-8 levels significantly increased post-treatment in individuals with AEs. The complete TLR signaling pathway highlighting the individual upregulated genes, KEGG pathways and transcription factors, was constructed with the use of the online database SPIKE [37] (S3 Fig). Finally we wanted to compare our newly identified LF AE transcriptional signature to published gene expression profiles. The post-treatment AE transcriptional signature of the 744 upregulated genes was very similar to multiple published endotoxin exposure gene expression profiles, in addition to many other profiles (S7 Table).
After successfully identifying a significant transcriptional signature of post-treatment AEs, we explored whether a pre-treatment transcriptional signature could predict what individuals would go on to develop AEs after treatment. There were no significant differentially expressed genes at baseline between individuals that would develop moderate AEs, and individuals that would not develop AEs (by DESeq2).

Neutrophils increased and lymphocytes decreased post-treatment
Changing cell populations can have a large effect on gene expression profiles. Differential cell counts were unavailable, so cell type proportions were estimated using the RNA-seq data. This analysis suggested that neutrophils increased more and lymphocytes decreased more posttreatment in individuals with AEs compared to individuals with no AEs (Fig 5A). These changes are consistent with stress-type immune responses. For simplicity, B cells (memory and naïve), T cells (CD8, CD4 naïve and memory resting) and NK cells were combined into one category (lymphocytes) in Fig 5A, and ungrouped data are presented in S4 Fig. Estimated leucocyte proportions at baseline were very similar between the two AE groups, but individuals who experienced post-treatment AEs had significantly fewer estimated memory B cells compared to individuals who did not develop AEs after treatment (Fig 5B).

Prioritization of genes upregulated post-treatment in individuals with AEs show that TLR2 is one of the most important genes for the development of AEs
The genes upregulated post-treatment in individuals with AEs were prioritized by importance for AE development using a machine-learning tool (RF analysis). This was done in order to identify genes with the strongest associations between expression levels and development of AEs and to identify genes of interest for PCR validation. Table 2 shows the top 15 genes that were upregulated in persons who developed moderate AEs after treatment. However, based on this analysis it was not possible to determine whether the gene expression changes were the cause or effect of the AEs that were experienced.

Orthogonal validation of expression levels for candidate genes confirmed RNA-seq data
qRT-PCR studies were performed to confirm whether expression of genes identified by DESeq2 and RF analyses was actually increased post-treatment in individuals with moderate AEs. Increased expression after treatment was confirmed for seven of the top eight genes (DIP2B, ZCCHC6, PELI1, FNDC3B, TLR2, LTBR and NT5C2) (Fig 6). Expression of the eighth gene (RBPJ) did not change after treatment in either AE group. Expression of the housekeeping gene HPRT1, did not change with treatment in either AE group (as expected).

Modeling identified high baseline CFA levels as predictor for development of post-treatment AEs and post-treatment increases in LBP levels as another risk factor
We were unable to identify any pre-treatment transcriptional signature that could predict moderate AEs. We therefore wanted to assess if any metadata or baseline infection parameter was associated with the development of moderate AEs. A logistical regression was performed to consider effects of age, sex, treatment arm, baseline Mf/mL and baseline CFA on the risk for development of post-treatment moderate AEs. A total of 71 individuals were included in the model (9 moderate AEs and 62 no AEs). The logistic regression model was statistically significant, X 2 (6) = 22.1, P = 0.0012. The model explained 50.2% (Nagelkerke R 2 ) of the variance in AE outcome, and correctly predicted 93% of outcomes. However the model was better at predicting people who did not develop AEs; it correctly predicted only 44.4% of the individuals  Post-treatment adverse events in lymphatic filariasis who developed moderate AEs. Increasing baseline CFA levels were associated with increased likelihood of developing AEs (P = 0.022), but the other independent variables did not significantly contribute to the model. RF analysis was performed on the same dataset (71 individuals), and this also identified the baseline CFA level as the best predictor for subsequent development of AEs. However, treatment arm was also a positive predictor in the RF model, and could therefore be related with the development of AEs (Table 3). It was surprising that baseline Mf count was not identified by the logistic regression model or RF to significantly contribute to correctly predicting AEs. However, baseline Mf counts were higher in individuals who developed moderate AEs (geometric mean 343 Mf/mL) compared to individuals with  Post-treatment adverse events in lymphatic filariasis no AEs (geometric mean 188 Mf/mL), and this difference was significant (P = 0.036 by Mann-Whitney U test). RF analysis was also used to identify the variable (CFA, CIC, C3, C4, FB or LBP) that was best at classifying AE outcome based on post-treatment fold change in the 18 matched case-control subjects. LBP changes after treatment was the only variable that was significantly associated with the development of AEs (S8 Table).

Discussion
This study looked at changes in proteins in plasma and changes in gene expression in PBL in persons who experienced moderate AEs following treatment of LF.

Changes in proteins in plasma in persons who experienced AEs after treatment
Filarial antigen levels increased in plasma after treatment in individuals with moderate AEs, and this agreed with our recently published results from a separate clinical trial [13]. Western blot results from this study showed that many new filarial antigens with the carbohydrate epitope detected by the monoclonal antibody AD12 appeared in the blood 24 hours after treatment in some individuals. In contrast, only a single high molecular weight antigen circulates in the blood of W. bancrofti-infected individuals without treatment [26]. We postulate that treatment kills or injures worms so that they release internal antigens that are normally concealed inside the parasite.
Results from this study also confirmed our previous finding that plasma CIC levels do not increase after treatment of LF in persons who develop moderate AEs [13]. This finding suggests that AEs after LF treatment are not caused by CIC. The complement cascade (classical pathway) modulates pro-inflammatory effects of CIC [38]. Activation of the classical complement cascade leads to decreased C3 and C4, whereas activation of the alternative pathway (AP) leads to decreased C3 and factor B (FB). Our results are most consistent with activation of the AP by parasite antigens, as C3 and FB decreased in individuals with moderate AEs while there was no change in C4 levels. The RNA-seq data also supports the AP hypothesis, because expression of CFP (complement factor properdin-a positive regulator and initiator of the AP) significantly increased post-treatment in individuals with moderate AEs (adjusted P-value 0.001, DESeq2). In contrast, expression of C4B (basic form of C4, part of the classical pathway) significantly decreased (adjusted P-value 0.04, DESeq2) post-treatment in individuals with moderate AEs. Additionally, IFN-γ and TNF-α are known to induce FB synthesis [39]. This could account for the inconsistent FB levels post-treatment in individuals with moderate AEs, because both IFN-γ and TNF-α increased in these individuals; the positive stimulus of these cytokines may have counteracted decreases in FB levels as it was used in the AP. In summary, many different filarial antigens were transiently released post-treatment, but they did not appear to form CIC or activate the classical complement cascade.
PBL appeared to respond to the release of filarial antigens by releasing cytokines, and the cytokine profiles of post-treatment AEs in LF infected individuals are complex. In our previous study of samples from a treatment trial in Papua New Guinea, we reported that 16 cytokines increased post-treatment in individuals with moderate AEs [13], and eight of these cytokines were also increased after treatment in this study (IL-1RA, IL-6, IL-10, G-CSF, MCP-1, MIP-1β, TNF-α, and VEGF). It is not surprising that more cytokines increased in the previous study, because more time-points were sampled. Also, participants in the Papua New Guinea study had higher blood Mf counts and higher rates and severity of AEs than participants in the present study, and this may account for their stronger cytokine responses. A new finding in this study was the increase in IFN-γ post-treatment in individuals with moderate AEs. This was supported by the fact that IRF1 (downstream of IFN-γ) was identified to be an important transcription factor for AE development. The increase in IL-8 levels paralleled the increase in neutrophils post-treatment in individuals with moderate AEs. Results from this study did not confirm the finding from the prior Papua New Guinea study that high levels of eotaxin-1 pretreatment are a risk factor for development of post-treatment AEs. Again this discrepancy may be related to differences in infection intensity between participants in the two studies. Levels of IL-6 and TNF-α have previously been shown to be positively correlate to levels of Wolbachia DNA in human plasma 48 hours after treatment of LF [40].

Changes in gene expression in PBL in persons who experienced moderate AEs after treatment
RNA-seq was performed to better understand changes in leukocyte gene expression that occur in persons who experienced moderate AEs after treatment. We identified a distinctive transcriptional signature associated post-treatment moderate AEs with 783 genes that were differentially expressed (at the P < 10 −5 level) in persons who experienced moderate AEs. In contrast, no gene was differentially expressed at that level of significance before or after treatment in individuals who did not experience AEs. 95% of the 783 genes associated with AEs exhibited increased expression post-treatment. Thus, moderate AEs were primarily associated with upregulation of gene expression and not with gene suppression. A total of 126 genes were upregulated post-treatment in individuals with no AEs at the low stringency P < 0.05 level, but 83% of these genes were also upregulated post-treatment in individuals with moderate AEs. Thus changes in gene expression after treatment did not always lead to clinically evident AEs.
The transcriptional signature results are consistent with the hypothesis that Wolbachia lipoprotein activates TLR2-TLR6 [18,41], as bacterial lipoproteins can induce pro-inflammatory responses through TLR2 signaling and NF-κB and STAT1 activation [42]. The finding that TLR2 was one of the genes most highly associated with the development of moderate AEs also supports this hypothesis. Furthermore, LBP was found to increase post-treatment in plasma from individuals with moderate AEs, and RF analysis identified LBP (fold change post-treatment) as the best variable for classifying AE outcome. LBP is an acute-phase protein that is mostly known for its function of shuttling LPS to TLR4 via CD14. However, it can also shuttle lipoproteins to TLR2 also via CD14 [43] as would be the case with PAL. CD14 expression was upregulated post-treatment in individuals with moderate AEs, and it had one of the most significant adjusted P-values (6.4e -26 ) in the dataset. CD36 is another accessory receptor for the TLR2-TLR6 heterodimer [44], and this gene was also upregulated post-treatment in individuals with moderate AEs (adjusted P-value 0.004). The LF AE transcriptional signature was similar to previously published endotoxin exposure gene expression profiles further supporting the Wolbachia lipoprotein hypothesis. Since multiple filarial antigens that are normally not accessible to the immune system are released post-treatment, it is possible that Wolbachia components are released into the host's circulation in a similar fashion. Thus we cannot be certain that Wolbachia lipoprotein is the only or prime trigger for AEs. It is possible that other filarial components can signal through TLRs and contribute to the development of AEs, and this might explain why expression of many other TLRs (TLR1, TLR4, TLR5, TLR6 and TLR8) were upregulated post-treatment in individuals with moderate AEs. Many different ligands can activate TLR2, including protozoan ligands such as GPI anchors from Trypanosoma cruzi [45] and Leishmania major [46]. Wolbachia-independent activation of the immune system causing severe AEs is seen in Loa loa (a filarial worm that lacks Wolbachia [1]) infected individuals post-treatment suggesting that Wolbachia is not the sole cause of AEs after anti-filarial treatment. TLR signaling is clearly associated with the development of AEs, but complement activation has similar downstream effects, and there is considerable crosstalk between these two pathways [47]. Thus both TLR signaling and the complement AP could be actively involved in the pathogenesis of AEs. Another possible mechanism for AEs following treatment of LF is that treatment abrogates the normally dominant Th2 immune responses stimulated by helminth infections that interfere with the expression and function of TLRs [48]. Increased TLR expression and signaling after treatment may then induce pro-inflammatory Th1 responses causing AEs. Indeed, classic Th1 cytokines TNF-α and IFN-γ were increased after treatment in people with moderate AEs.
We did not detect a pre-treatment transcriptional signature that was a significant risk factor for development of post-treatment AEs. Baseline CFA levels were the best predictor for moderate AEs in this study, and this was the only variable that was a significant predictor by logistic regression. CFA levels are correlated with Mf counts and likely with adult worm numbers, and this result suggests that CFA levels are related to the total parasite biomass that can potentially contribute to the development of AEs after treatment.

Genes upregulated post-treatment in individuals with AEs
An important finding from this study was that post-treatment AEs in LF-infected individuals are associated with upregulation of hundreds of genes. A prioritized list of the top 15 genes important for the development of AEs is listed above in Table 2. These genes and their associated pathways may provide insight into the pathogenesis of AEs. In addition to the TLR pathway, Notch, NF-κB and IL-1 signaling were common themes. Four of the top 15 genes (RBPJ, TLR2, ALDH1A2 and APLP2) are involved in TLR/Notch signaling. The Notch pathway is involved in development and is conserved from Drosophila to mammals. RBPJ is involved in the crosstalk between TLR and Notch signaling that is thought to help fine-tune the immune response through negative and/or positive feedback (S5 Fig) [49]. NF-κB is another downstream pathway of Notch signaling, and three of the top 15 genes (PELI1, TLR2 and LTBR) are involved in NF-κB signaling. TLR2 is a receptor for the canonical NF-κB pathway, and PELI1 is involved in intracellular downstream signaling. The canonical pathway results in the release of pro-inflammatory cytokines, such as IL-6 and TNF-α. Thus, activation of this pathway is consistent with the cytokine profiles individuals who experienced moderate AEs after treatment. Interestingly, lymphotoxin beta receptor (LTBR) stimulation has been shown to enhance the LPS-induced expression of IL-8 via the combined action of NF-κB and IRF1 [50], and this is consistent with our results. Finally, IL-1 is clearly associated with the development of AEs. IL-1RAP was one of the top 15 genes identified by RF, and this gene is one of two coreceptors for IL-1. Additionally, expression of the second co-receptor IL-1R1 was increased post-treatment in individuals with AEs. Both expression (this study) and plasma levels (our prior study [13]) of IL-1β increased after treatment in individuals with AEs. Inhibitors of the IL-1 pathway were also upregulated post-treatment in individuals with AEs. This included both increased expression and protein levels of IL-1RA and increased expression of the IL-1β decoy receptor IL-1R2. These results illustrate the importance of the balance between the proinflammatory effects of IL-1β and the anti-inflammatory effects of IL-1RA for AE development.

Limitations of the study
One limitation of this study was that the RNA-seq was performed on mixed PBL samples. This makes it difficult to separate the effects of altered gene expression from the effects of changing cell type proportions. Separating different types of leukocytes was not feasible, because the samples were collected in rural Côte d'Ivoire and processed in a simple field lab. On the other hand, this study provides insight into the pathogenesis of post-treatment AEs. Additionally, it was possible to estimate the different cell subtypes present in PBLs using the RNA-seq data, so we could associate changes in gene expression with altered cell types to decrease the chance that the former was a directly result of the latter. For example, if the post-treatment AE transcriptional signature had been similar to a neutrophil gene expression profile it could have been caused by the increasing proportion of neutrophils post-treatment and not due to specific neutrophil activation during AEs. Another limitation was that we did not study an untreated control group, because it would have been unethical to withhold treatment from infected individuals. For the cytokine analysis we did not correct for multiple comparisons, however, the results are generally consistent with our previous findings [13], and this increases our confidence in the results. Based on the standard significance level of 0.05, approximately 3 differences would be expected to be significant by chance for 54 tests (27 cytokines measured preand post-treatment in two AE groups), whereas 17 comparisons were significant at the 0.05 level in this study. The semi-quantitative Filariasis Test Strip (FTS) was used to assay CFA in the field, whereas a quantitative ELISA was used to measure CFA levels in the laboratory setting in this study. Newer studies have demonstrated that tests that detect LF CFA (including FTS) can cross-react with L. loa antigens that circulate in blood from a subset of individuals with heavy infections [51], and with biological samples from animals infected with L. loa and Onchocerca ochengi [52,53]. This cross-reactivity was not a concern for this study because L. loa is not endemic to Côte d'Ivoire, and the area of Côte d'Ivoire where the study was conducted is non-endemic for O. volvulus.
In future studies it would be interesting to measure Wolbachia DNA in pre and post-treatment samples from individuals with AEs after treatment of LF and onchocerciasis to try to correlate bacterial DNA release with host expression profiles post-treatment. Wolbachia DNA has been shown to increase post-treatment in individuals with AEs after treatment of LF [14]. Additionally, peak Wolbachia DNA levels have been shown to be correlated with AE reaction scores in individuals treated with DEC or IVM for onchocerciasis [54].

Conclusions
This study included first global RNA-seq analysis of PBLs from LF-infected individuals, and it has provided novel insights into the pathogenesis of a clinically relevant problem. The samples were ideal for studying AEs after LF treatment, because each post-treatment sample was paired with a pre-treatment sample from the same individual. This internal control improved our ability to study the AE phenotype in humans. AEs represent a significant challenge for the global program to eliminate LF, and the fear of AEs in communities receiving MDA is a main factor that reduces compliance [25]. Minimizing the impact of AEs has therefore been identified as a key component for successful MDA programs [25]. More than 850 million individuals have been treated as part of GPELF, and a significant percentage of these individuals experience AEs.
This study has also provided a framework for investigating the host responses associated with severe AEs that occur after treatment of other filarial worms such as O. volvulus and L. loa. Treatment of other, more familiar infections can also result in severe AEs that are caused by host responses to dying pathogens. This Jarisch-Herxheimer reaction occurs after antibiotic treatment of spirochetal infections such as syphilis, Lyme disease, leptospirosis, and relapsing fever, and it is also hypothesized to be caused by the release of bacterial lipoproteins that activate TLR2 [55]. The transcriptomic response during the Jarisch-Herxheimer reaction has not been studied, so the dataset from this study could provide a valuable starting point for research on this related clinical problem.
To recap our major findings, this study has provided new insights regarding the pathogenesis of post-treatment AEs in LF-infected individuals. Our results are consistent with the hypothesis that a Wolbachia lipoprotein triggers AEs by binding to TLR2-TLR6, but other uncharacterized filarial antigens might also play a role. Since TLR, NF-κB, and TNF pathways are involved, these pathways could potentially be targeted to prevent or treat AEs after LF treatment. We also found that high pre-treatment CFA levels were the best predictors of posttreatment AEs. This finding could be relevant for treatment-naïve areas with high LF infection prevalence and intensities. Individuals with high CFA levels pre-treatment (assessed with the FTS [56]) could be offered non-steroidal anti-inflammatory medications together with antifilarial medications for home management of moderate or severe AEs. However, a positive FTS from an individual who resides in or has traveled to an area that is also endemic for L. loa needs to be interpreted with caution due to the issues of cross-reactivity mentioned above. Information from this study should allow program managers and drug distributors to reassure populations and communicate to them that AEs experienced after LF treatment are transient and caused by host responses to dying or injured parasites.    Table. Upregulated KEGG pathways pre-treatment in individuals that did not develop adverse events (Gene Set Enrichment Analysis). A KEGG pathway. Size: Number of genes in the gene set after filtering out those genes not in the expression dataset. ES: enrichment score for the gene set; that is, the degree to which this gene set is overrepresented at the top or bottom of the ranked list of genes in the expression dataset. NES: normalized enrichment score; that is, the enrichment score for the gene set after it has been normalized across analyzed gene sets. NOM P-val: Nominal P-value; that is, the statistical significance of the enrichment score. The nominal P-value is not adjusted for gene set size or multiple hypothesis testing; therefore, it is of limited use in comparing gene sets. FDR q-val; false discovery rate; that is, the estimated probability that the normalized enrichment score represents a false positive finding. FWER Pval: Familywise-error rate; that is, a more conservatively estimated probability that the normalized enrichment score represents a false positive finding.