|Home | About | Journals | Submit | Contact Us | Français|
Phenotypic manifestations of infectious diseases are closely related to individual immune responses. Methods to extract information from patients’ own immune reactions would be of great use for both diagnosis and treatment. Dengue fever is one of the diseases that clinical aggravations could occur paradoxically after humoral immunity appears. This property makes dengue fever an excellent disease model to explore. A principal component analyses (PCAs)-based framework derived from a prior vaccination study was developed. The framework was verified by successful demonstrations of known IgG signatures from a Mexico Dengue data set. Afterward the pipeline was tested upon de novo IgG and IgA libraries of Dengue patients from southern Taiwan. We discovered four infection signatures within IgG repertoires, two of which were identical to previous reports. However, it was IgA but not IgG that could differentiate hemorrhagic from non-hemorrhagic patients. IgA repertoires were found more diversified among bleeders, from whom seven signature clusters were characterized. The expressions of transforming growth factor beta 1 (TGFβ1) and accordingly mediated class-switch activity of IgA were distinct only among the PCA-segregated bleeding group. In sum, intercontinental sharing of IgG signatures in dengue fever was demonstrated via a unified working flow. Differential regulation of IgA class-switch with associated diversity expansion plus existences of hemorrhage-restricted clusters were shown. The ability of the framework to find common IgG signatures would implicate applications to infections even from unknown pathogens. The clusters within IgA repertoires could offer perspectives to other IgA-related bleeding disorders such as Henoch-Schönlein purpura or IgA nephropathy. Substantiated grounds for IgA-specific effector function via TGFβ1-mediated class-switch would be a new factor to consider for infectious diseases.
Humoral immunity could be both friends and foes in human diseases. For example, IgG and IgE antibodies specific for double-stranded DNA could differentially induce pathogenic inflammation in systemic lupus erythematosus (1). Reactions elicited by past viral infections could have effects beyond the same species of viruses (2). Therefore, an effective and versatile pipeline to distil signals out of humoral immunity could be of great value for both basic research and clinical interpretation of human diseases. Dengue fever is an acute febrile illness caused by four groups of dengue viruses (3). Infected mosquitoes of either Aedes aegypti or Aedes albopictus are the principal vectors in transmitting the pathogens to 4 millions of people in the tropical and sub-tropical areas every year (4). Clinical courses of dengue fever can be divided into three phases: febrile, critical, and convalescent (4). Most patients recover spontaneously, but a few suffer from hemorrhages, plasma leakage, or even circulatory collapses at the critical stage. These severe forms of dengue fever usually occur after the febrile stage, when viral loads are actually very low (5). The most well-recognized factor that significantly increases the likelihood of these serious consequences is the secondary heterotypic infection. Aberrant adaptive immunity might play roles in these scenarios. For example, cross-reactive antibodies against dengue NS1 protein can also induce apoptosis of endothelial cells (6) or enhance activation of plasminogen (7). In addition, the so-called antibody-dependent enhancement (8, 9) hypothesis has been proposed to explain overactivation of myeloid cells after dengue viral infection. Therefore, dengue fever is an excellent disease model for analyses of humoral immunity or antibody repertoires in an infection- and phenotype-specific context.
Recent innovations of next-generation sequencing (NGS) have made possible clonal examinations of adaptive immune responses in dengue fever. Because nucleotides in the complementarity-determining region 3 of the heavy chain (CDR-H3) on most antibodies are sufficient to determine specificities (10), sequence repertoires of this region can effectively serve as clone proxies of humoral immunity. Recently, Parameswaran et al. described convergent IgG signatures among dengue-recovered patients in Nicaragua (11), and Godoy-Lozano et al. found fewer somatic hypermutations among IgG immune repertoires in Mexico (12). It is not clear, although, if there exist specific repertoire signatures that are linked to any of the severe forms of dengue fever.
In this study, we adapted a selection heuristic of antibody repertoires that was developed to characterize carrier children of chronic hepatitis B and vaccination responses of healthy siblings on the platform of NGS (13). After reproducing the same result of the infection signature in the hepatitis B data set by the new pipeline, we applied the renewed scheme to a public data set of dengue fever from Mexico. The infection signatures as reported previously (11, 12) were successfully identified. We then collected blood samples from Taiwan’s Dengue patients. The heuristic indeed revealed four infection signatures, two of which were identical to prior discoveries in Nicaragua and Mexico (11, 12). We further compared IgG and IgA immune repertoires between patients with or without hemorrhages. It was noticed that alterations of diversity profiles in IgA repertoires were more prominent among patients with bleeding phenotypes. Instead IgG repertoires were unable to tell the differences. The identical cluster-seeking heuristic then spotted seven IgA clusters that were closely linked to the clinical bleedings. Furthermore, we found transforming growth factor beta 1 (TGFβ1)-mediated class-switch activity of IgA was differentially regulated in the principal component analysis (PCA)-segregated hemorrhagic group only. In conclusion, our results demonstrated an efficient pipeline to characterize the stochastic IgG signatures of viral infection and revealed the hidden relationships between IgA and hemorrhages in dengue fever. Host variations in TGFβ1-mediated switch to a more diversified profile of IgA effectors could be a novel contributing factor to dengue hemorrhages.
The protocol in recruiting dengue subjects (Table (Table1;1; ND5, SD1-7, and HD1-7) was approved by the Institutional Review Board of Kaohsiung Medical University Hospital. Some of non-dengue controls (Table (Table1;1; ND1-4) were from a previous study (13), the protocol of which was approved by both the Institutional Review Board of National Health Research Institutes and the Institutional Review Board of En Chu Kong Hospital. Study methods were carried out in accordance with the relevant guidelines. Informed consents were obtained from all attendants.
Blood samples were collected, handled, and discarded in accord with clinical standards. NS1 protein of dengue viruses was tested with the enzyme-linked immunosorbent assay. The viral RNA was extracted from serum with the QIAamp viral RNA kit (Qiagen) according to the manufacturer’s instructions. After reverse transcription, a TagMan real-time PCR was carried out on an ABI thermocycler with a published protocol (14) for serotyping.
Total RNA from 2.5-ml whole blood in the PAXgene tube was extracted with PAXgene Blood miRNA Kit (QIAGEN). 1.2µg was reverse transcribed by IgG or IgA constant region-specific primers (Table S4 in Data Sheet S2 in Supplementary Material) with SuperScript III Reverse Transcriptase and RNaseOUT (Life Technologies), going through 65°C 5min, on ice at least 1min, 55°C 60min, and final 70°C 15min.
Normalized cDNA were multiplex amplified with a set of forward V-primers from known IgG alleles (15) that were Primer 3 (16) optimized by the same parameters and two reversed primers (Table S4 in Data Sheet S2 in Supplementary Material). All forward primers shared similar amplification efficiencies with linear correlations to cDNA concentrations across 5 logs and had low backgrounds as well as specific products in pilot optimizations. The reaction condition was 98°C for 2min, 15 cycles of 98°C for 80s+60°C for 60s+65°C for 30s, 10 cycles of 94°C for 30s+65°C for 90sec, and final 65°C for 10min with AccuPrime Taq DNA Polymerase High Fidelity (Life Technologies).
The product was amplified in a second PCR by the same polymerase with double-indexed P5 and P7 primers (Table S4 in Data Sheet S2 in Supplementary Material) under the condition 94°C for 2min, 10 cycles of 94°C for 30s+60°C for 30s+68°C for 40s, 10 cycles of 94°C for 30s+72°C for 90s, and final 72°C for 5min. Products were sieved by 2% agarose gel under 30V × 8h or 60V × 4h, and the target bands around 300bp were eluted with MinElute Gel Extraction Kit (QIAGEN) before 150bp paired-end sequencing (NextSeq, Illumina, USA).
Demultiplexed raw sequences were processed as reported before (13). Briefly, sequences were paired with PEAR (17); poor quality reads were removed at this step with default parameters, including a Phred filter score at 33. Reads must have no ambiguous “N” nucleotides as well as matched end sequences on both 5′ and 3′ terminals to the amplifying primers. The expected CDR-H3 region as implicated by matched V-primers had to contain multiplicities of three nucleotides without stop codons. The N-end of translated amino acids from position 100 to C104 had to align well to anticipated sequences as suggested by matched V-primers and the C-end had to have the W118 as well as the following GXG signatures, where X denotes any amino acids (18). At last, CDR-H3 without the minimal length of two amino acids was discarded.
For CDR-H3 libraries from the Mexico data set, the raw sequencing files were downloaded and subjected to MIXCR alignment (19). Identified CDR-H3 by MIXCR contained the heading C104, trailing W118, and occasional non-amino acid letters like “*” or “_” in between. We preprocessed these CDR-H3 by trimming C104/W118 and discarding “*” or “_” containing reads before bridging to the downstream analyses.
cDNA was synthesized with random hexamers primed to the same RNA used in preparations of CDR-H3 libraries. TGFβ1 was qPCR-assayed (KAPA SYBR FAST Master Mix, Kapa Biosystems) with a published primer pair (20), 5′-CCCAGCATCTGCAAAGCTC-3′ and 5′-GTCAATGTACAGCTGCCGCA-3′, under the condition 95°C for 3min, 40 cycles of 95°C for 1s+60°C for 25s, and a final dissociation stage on 7900HT Fast Real-Time PCR System (Applied Biosystems). IgA GLT was first PCR amplified from RNase-H cleaned cDNA with primers 5′-CAGCAGCCCTCTTGGCAGGCAGCCAG-3′ and 5′-TTTCGCTCCAGGTCACACTGA-3′ (21) under the condition 94°C for 2min, 14 cycles of 94°C for 30s+60°C for 30s+68°C for 30s, and 68°C for 5min with AccuPrime Taq DNA Polymerase High Fidelity (Thermo Fisher Scientific). The product was then nested qPCR-assayed (KAPA SYBR FAST Master Mix, Kapa Biosystems) with primers 5′-TTGGCAGGCAGCCAGACG-3′ and 5′-TGGGGCTGGTCGGGGATG-3′ under the condition 95°C for 3min, 40 cycles of 95°C for 1s+60°C for 20s, and a final dissociation stage with the same qPCR reagent and instrument as used for TGFβ1. The nested primer pair was based on Sanger sequences obtained from outer PCR products, and the final nested products were confirmed by Sanger sequencing as well. Cycles of thresholds (Cts) for both transcripts were normalized to Ct of Actin. Expressions of TGFβ1 and GLT were compared with Mann–Whitney U and Student’s t-tests, respectively. Of note, GLT comparisons excluded patient SD6 to control the time spans between first and second samples within 5days (Table (Table1).1). The significance of the p value of 0.48 for the hemorrhage-absent central group was found <0.001 by fitting against the p value distribution exhaustively built from all p values at n=5 from the hemorrhage-prevalent peripheral group.
IgA-specific first-strand cDNA from patient HD3 as used for preparations of CDR-H3 libraries was assayed with the following primers to quantify total IgA1 vs. IgA2 and cluster-specific IgA1 vs. IgA2. For the former, a common forward primer 5′-ACCAGCCCCAAGGTCTTCC-3′ was paired with 5′-GATGACCACGTTCCCATCTG-3′ and 5′-GACGACCACGTTCCCATCTT-3′ reverse primers to detect total IgA1 and IgA2, respectively. A cluster-specific forward primer 5′-TGCGACGGTCTTCACTACAG-3′ was chosen among the most abundant CDR-H3 sequences of cluster ATVFTTVHY to examine cluster-specific IgA1 and IgA2 with the same reverse primers above. The reactions were performed in triplicates with KAPA SYBR FAST qPCR Master Mix (Kapa Biosystems) on 7900HT Fast Real-Time PCR System (Applied Biosystems) under the condition 95°C for 3min, 40 cycles of 95°C for 1s plus 60°C for 10s, and final addition of a dissociation stage to confirm product specificity.
3′ extension sequences after the two reverse primers paired with V-primers (Table S4 in Data Sheet S2 in Supplementary Material) were used to classify NGS reads into IgG subclasses. 5′-GCCCTTGGTGGAG-3′ would categorize reads into the composite group of IgG1+IgG2, while 5′-GCCCTTGGTGGAA-3′ would segregate reads into the other group of IgG3+IgG4. Relative percentages of cluster-specific reads among both groups were calculated for all dengue samples. The results of two series were subjected to the Mann–Whitney U test, which yielded p=0.28.
The normalized clone frequencies for each sample were Hellinger transformed (22) before PCA (23) in SciPy (24), where “svd_solver” was set to “full.” Clones ranked beyond 2 SDs in absolute values of loadings were used in clustering as specified in the text and below.
Samples with more reads were randomly resampled with Python to match the repertoire with minimal read counts in the indicated data set. Rarefactions were performed 10 times and saved for downstream analyses. Clusters were identified in two steps. Indel-free Hamming distances (25) between clone pairs were calculated into an adjacency matrix in Python, which was used to initiate a graph in igraph (26) discarding edges weighted higher than distance 1. Cluster-associated functions in igraph were applied to discover independent clusters and to calculate associated PageRank scores (27).
For repertoire A and B with clone frequencies FA and FB, the Morisita dissimilarity index (28) is a measure to quantify the distance between two sets of clone sequences. This measure is affected mainly by abundant clones; relatively rare clones have little effect, even if there are many of them. The measure (M) was calculated as follows:
In each rarefaction data set, samples with the same hemorrhagic phenotype and immunoglobulin class were pooled. Within each pool, clusters were defined by an adjacency matrix as described above. We then used Hill numbers, i.e., effective number of clusters in this case, to quantify cluster diversity. Diversity profiles in Hill numbers (D) (29) were calculated following the formula with frequencies normalized in each pool:
The parameter q determines each measure’s sensitivity to normalized frequencies. The measure of q=0 (the total diversity) counts clusters equally without regards to their normalized frequencies. The measure of q=1 [Shannon diversity, the exponential of Shannon entropy (30)] counts clusters in proportional to their normalized frequencies and thus can be interpreted as the number of common clusters in the data. The measure of q=2 discounts all but the dominant clusters and can be interpreted as the number of dominant clusters in the data. The plot of Dq with respect to the parameter q is referred to as a “diversity profile” in ecological science. The profile is generally a decreasing function. The slope of the curve reflects the unevenness of cluster-normalized frequencies. The more uneven the distribution of normalized frequencies, the more steeply the curve declines. Graphs with q=(0, 5] were generated with igraph (26).
In this study, significantly contributing clones in PCA were selected as building blocks for clusters. In the previous work of chronic hepatitis B (13), carrier and non-carrier clone pools were instead made of 2-occurrence CDR-H3 sequences among the four carriers and the four non-carriers, respectively, regardless of vaccination histories. Clusters were then defined among the indicated clone pools with adjacency matrices before downstream feature selections.
Linear support vector classification (SVC) and logistic regression (LR) from Scikit-learn (31) were used to select clusters with “l1” penalties. Hellinger transformed frequencies served as independent variables for classification models. The top 1 or 0.5% of clusters with most abundant members were candidate clusters subjected to selections. In choosing clusters marking Dengue or hepatitis B infections (Figures (Figures1C1C and and2B;2B; Figure S1B in Data Sheet S1 in Supplementary Material), repertoires from confirmed patients were labeled positive. For hemorrhage-associated clusters, repertoires from bleeding patients were marked positive. Regulatory parameters were optimized by leave-one-out cross validations. Clusters were selected if both SVC and LR models gave positive supports. In rarefied experiments, clusters that manifested only in one data set were discarded.
The data set of Mexico Dengue samples is available at NCBI-SRA repository, BioProject ID PRJNA302665. The data sets of Taiwan samples for healthy controls and dengue samples are available at European Nucleotide Archive, PRJEB9332 and PRJEB13768, respectively.
A previous study in characterizing antibody repertoires among carrier children of chronic hepatitis B found that PCA (23) could readily separate carrier from non-carrier repertoires in a non-supervised manner (13). We hypothesized that those CDR-H3 sequence clones (briefed as clones in the following text) with higher absolute loading values in PCA would be sufficient to spot the infection signature (13) (Figure S1A in Data Sheet S1 in Supplementary Material). We identified 20,429 and 41,106 clones from carriers and non-carrier repertoires, respectively, that had absolute loading values beyond 2 SDs in the PCA. These clones were each categorized by one amino acid difference criterion, i.e., Hamming distance 1 without indels (25), into 17,392 and 32,083 clusters, respectively. The top 0.5% of clusters with most abundant members in either category were subjected to feature selections by both SVC and LR. Among a total of 246 clusters, only 1 (Figure S1B in Data Sheet S1 in Supplementary Material) gained very strong supports from both models (Table S1 in Data Sheet S2 in Supplementary Material), and it happened to be the previously published infection signature (13) comprising 28 unique clones (Figure S1C in Data Sheet S1 in Supplementary Material). We expected that the same selection scheme could be generalized to other diseases as well. The idea was tested with a public data set of dengue fever from Mexico (12). We found that repertoires from acute dengue and postconvalescent samples overlapped to a large extent on the PCA plot (Figure (Figure1A).1A). To improve the contrast, we mixed IgG repertoires from 4 healthy adults and 4 healthy children in the vaccination study (13) and 19 acute dengue samples from the Mexico data set. We further conducted 10 runs of rarefactions to match the lower sequencing depths in the dengue data set. Effectively the acute dengue repertoires separated fairly from those of adults but not from children (Figure S2 in Data Sheet S1 in Supplementary Material). A typical PCA plot from one of the rarefied data set is shown (Figure (Figure1B).1B). Candidates were set to the top 1% of member-rich clusters as constructed from those clones contributing variations most on PCA plots. Intersected results from both SVC and LR models (Table S2 in Data Sheet S2 in Supplementary Material) yielded nine signature clusters (Figure (Figure1C),1C), which appeared at least twice among 10 rarefactions (Table S3 in Data Sheet S2 in Supplementary Material). Each cluster was designated with the clone sequence that had the highest PageRank score (27). There were 143, 9, 15, 9, 7, 140, 16, 7, and 9 members for cluster ARQRGNWFDS, VGGGHYGP, ARGLNFLDSSGYHQSNWFAP, ARVSFNWNDVNYYYYMDV, AKDCGLGGDRDY, ARDLRYCSGGSCYDGAFDI, ARRRPLWPSYYFDY, ARLDYYYYYGMDV, and ARGGGFAWYNFDN, respectively (Table S3 in Data Sheet S2 in Supplementary Material). The most prevalent cluster, ARQRGNWFDS, occurred in all rarefactions (Table S3 in Data Sheet S2 in Supplementary Material) and was the same as reported by the original dataset contributors (12). One of the clusters, ARLDYYYYYGMDV, was the convergent dengue signature reported by Parameswaran et al. (11). On the basis of these establishments, we set out to investigate antibody repertoires of dengue fever in Taiwan.
15 feverish patients from a single medical center in southern Taiwan were recruited to the study (Table (Table1;1; ND5, SD1-7, and HD1-7). Of the 15 patients, 14 had dengue NS1 protein detected in blood, but the other did not. Half of the dengue-confirmed cases had clinical evidences of hemorrhages plus other signs like thrombocytopenia or elevated hepatic transaminases (Table (Table1;1; HD1-7), meeting at least the criterion of dengue fever with warning signs (4) in terms of severity. However, none of the recruited subjects suffered from plasma leakage or circulatory collapses. All of the 15 patients contributed blood samples at two different time points (Table (Table1;1; referenced from the disease onset). There were no significant differences of age, sex, first sampling day, second sampling day, lowest platelet count, and highest AST or ALT values between both groups (ANOVA, minimal p=0.13).
We prepared CDR-H3-based immune repertoires with optimized PCR primers (Table S4 in Data Sheet S2 in Supplementary Material) for parallel sequencing. On the basis of total RNA without enrichments, we got average reads of 792,819±475,770 (SD) for IgG libraries and 847,278±442,526 (SD) for IgA libraries; unique sequence clones amounted to 45,498±20,894 (SD) for IgG repertoires and 52,033±22,358 (SD) for IgA repertoires, respectively (Table S5 in Data Sheet S2 in Supplementary Material). All comparisons between hemorrhagic and non-hemorrhagic patients in terms of counts of IgG reads, IgG clones, IgA reads, and IgA clones at both time points yielded insignificant results except the counts of IgG clones at the second time point (Student’s t-test, p=0.02; minimal p=0.29 for the other comparisons).
To identify infection signatures of dengue fever in Taiwan, the four IgG adult control repertoires as used above with the Mexico data set were again combined with eight dengue IgG repertoires in a strict chronological order of case recruitment, without exclusions, for PCAs. The repertoires at the second time point for each dengue patient were adopted because longer infection exposures presumably would better mark antibody repertoires. We found that dengue repertoires could be fairly separated from healthy controls (Figure (Figure2A).2A). The two components explained 10.1 and 9.4% variances, respectively. A set of 13,548 unique clones whose absolute loading values passed beyond 2 SDs in PCA was grouped into 10,853 clusters with the Hamming distance 1 criterion without indels. The top 1% of clusters with most abundant members were subjected to feature selections by both SVC and LR. Clusters were selected if both SVC and LR models gave positive supports. We found four clusters satisfied the above filtering conditions (Table S6 in Data Sheet S2 in Supplementary Material). Each cluster was designated with the clone sequence that had the highest PageRank score (27). There were 321, 10, 6, and 14 members for cluster ARLDYYYYYGMDV, ATAFTTVDY, ARQYGNYFDY, and ARANVRNHIYSSSWAYFDY, respectively (Table S7 in Data Sheet S2 in Supplementary Material). Cluster ARLDYYYYYGMDV and ARQYGNYFDY were reported before (11, 12), but the others were not. Except cluster ARLDYYYYYGMDV, topographies of all other three clusters were relatively simple (Figure (Figure2B).2B). We calculated the percentages of NGS reads for each cluster in the data set (Figure (Figure2C).2C). Although the four clusters were selected from a half of dengue patients who were recruited earlier, the other half of patients also had significant percentages of readings for cluster ARLDYYYYYGMDV (Figure (Figure2C).2C). Many of the patients had that cluster in their IgG repertoires at the first time point already (Figure (Figure2C).2C). Cluster ATAFTTVDY, ARQYGNYFDY, and ARANVRNHIYSSSWAYFDY were not shared among subjects, but the latter two had been detectable in the first IgG repertoires for indicated patients (Figure (Figure2C).2C). Non-dengue controls had no discernible presence of any of these four clusters in their IgG repertoires (Figure (Figure2C).2C). We further segregated members of cluster ARLDYYYYYGMDV into two composite groups of either IgG1+IgG2 or IgG3+IgG4. The same categorization was also performed on all IgG reads at both time points. Relative percentages for cluster ARLDYYYYYGMDV among both groups were illustrated (Figure (Figure2D).2D). There were no subclass preferences with p-value at 0.28 by Mann–Whiney U test.
IgA levels in body fluids are known to fluctuate along the course of dengue fever (32, 33). Because cross-reactive IgG antibodies have been reported to cause undesired consequences (6, 7), the pathogenic roles that might be attributable to IgA immunoglobulins are worth further investigations. First, we used Morisita index (28) to compare dissimilarities between the same isotypes for each dengue patient at both time points (Figure (Figure3A).3A). We did not see significant differences for either type of immunoglobulins, although hemorrhagic patients seemed to have higher dissimilarities for IgA (Figure (Figure3A).3A). We then compared dissimilarities between IgG and IgA antibodies (Figure (Figure3B).3B). The results showed that IgG and IgA were more dissimilar for those patients with hemorrhages (Figure (Figure3B).3B). The average Morisita indices for bleeders and non-bleeders were 0.92±0.02 (SE) and 0.78±0.07 (SE), respectively. Mann–Whitney U test exhibited a significant difference (p=0.028). A numerical assay to absolutely evaluate the status shifts of IgA immune repertoires can provide more insights than relative comparisons between IgG and IgA alone. Diversity profiles in terms of Hill numbers that are widely used in numerical ecology (29) could be useful in this context (13, 34).
10 runs of random rarefactions were conducted to normalize NGS reads among different samples. In each rarefied data set, samples were pooled together by both the clinical phenotypes of hemorrhages vs. non-hemorrhages and immunoglobulin classes of IgG vs. IgA. Diversity profiles in Hill numbers were plotted for each data set (Figure S3 in Data Sheet S1 in Supplementary Material), where clusters as defined by the indel-free Hamming distance 1 criterion were compared to ecological species in the calculation (29). All of the plots were nearly identical (Figure S3 in Data Sheet S1 in Supplementary Material), and the one from rarefied data set 1 was illustrated (Figure (Figure3C).3C). Upward shifts of the curves from hemorrhagic patients were noted for both IgG and IgA, but the latter apparently had a stronger momentum (Figure (Figure3C).3C). The trend was especially prominent for more abundant clusters as revealed by enlarging gaps along the increasing parameter of Hill numbers (Figure (Figure3C).3C). Therefore, the bleeding phenotypes in dengue patients were associated with a diversity expansion of immunoglobulins, of which IgA seemed to be more dominant.
With differential diversity shifts of IgG and IgA in hemorrhagic patients (Figure (Figure3C),3C), we hypothesized that these two classes of immunoglobulins might relate differentially to the bleeding phenotypes. We used PCA to summarize IgG and IgA variations among all 15 patients with samples at both time points (Table (Table1;1; ND5, SD1-7, and HD1-7). We found IgG repertoires at either time point were unable to classify patients in accord with their hemorrhagic phenotypes (Figure (Figure4A).4A). IgA repertoires instead were distinctly clustered into several groups; most non-hemorrhagic patients were clearly aggregated to the central area on PCA plots (Figure (Figure4B),4B), including the NS1-negative patient ND5. Clustering in PCA coordinates were significant at p=0.001 by the PERMANOVA test (35).
To stress the clustering power, we did IgA PCA with the 10 rarefied data sets as prepared above (Figure S4 in Data Sheet S1 in Supplementary Material). Although the repertoires from the first time point became suboptimal in resolving patient groups, the repertoires from the second time point were still as robust as those from the full data set (Figure (Figure4B).4B). Therefore, the clones with high loading values in PCAs as constructed from the second IGA repertoires could be closely associated with the hemorrhagic phenotypes.
We took a similar strategy used in identifying infection signatures of IgG repertoires (Figures (Figures11 and and2)2) to analyze those clones with highest absolute loading values in PCAs of the second IgA repertoires. For each data set, these clones were clustered first with the indel-free Hamming distance 1 criterion (29). The top 0.5% clusters with most abundant clone members were subjected to both SVR and LR feature selections. We found seven clusters that gained strong supports from both SVR and LR models in all 10 rarefactions (Table S8 in Data Sheet S2 in Supplementary Material). Each cluster was designated with the clone sequence carrying the highest PageRank score (27) among pooled CDR-H3 sequences of the indicated cluster from all 10 data sets (Figure (Figure5A).5A). The member counts for cluster ATVFTTVHY, ARPRTTRGGAPDV, ARDPVRWWSPSRGLYYYYMDV, VRGPDGQYGMDV, ARDCSGGNCYGSSYYGMDV, AKEAHTWNDVAGLDV, and AKFASSGSYADI were 106, 66, 53, 47, 38, and 50, respectively (Table S9 in Data Sheet S2 in Supplementary Material). Only among hemorrhagic patients did clusters have high read percentages in the second IgA repertoires (Figure (Figure5B).5B). Cluster ARPRTTRGGAPDV and AKEAHTWNDVAGLDV were already detectable in the first IgA repertoires for indicated patients. There was no sharing of clusters among hemorrhagic patients. Of note, IgA cluster ATVFTTVHY (Figure (Figure5A)5A) and IgG cluster ATAFTTVDY (Figure (Figure2B)2B) were essentially related (Tables S7 and S9 in Data Sheet S2 in Supplementary Material), but the topology in the latter (Figure (Figure2B)2B) was much simpler. We suspected that differential class-switch activity of IgA might be the underlying driving force for the appearance of the more complex IgA variant along with more diversified IgA profile (Figure (Figure3C)3C) among the bleeding patients.
We first used real-time PCR to quantify IgA subclasses. Cluster ATVFTTVHY-specific IgA1 vs. IgA2 was compared with total IgA1 vs. IgA2 on the second sample of patient HD3 (Figure (Figure5C).5C). Results showed that IgA1 was the more dominant subclass. The difference of Ct for total IgA was 3.1, but the Ct gap for cluster ATVFTTVHY got higher to 6.0 (Figure (Figure5C).5C). We then set out to assay global switch activity of IgA as induced by TGFβ (36). Dengue patients with hemorrhages are known to have higher expressions of TGFβ (37). We used qPCR to compare blood TGFβ1 levels between the hemorrhage-prevalent peripheral group on the IgA PCA plot (Figure (Figure4B4B and Table Table1;1; SD1,6,7, and HD1-7) and the hemorrhage-absent group at the center (Figure (Figure4B4B and Table Table1;1; SD2-5 and ND5). Indeed shortly after dengue onset the peripheral group had much higher TGFβ1 levels than the central group (Figure (Figure6A;6A; p<0.029, U test). However, the difference was no longer discernible in a few days (Figure (Figure6A;6A; p=0.38, U test). Because TGFβ1 is an IgA-specific class-switch factor (36) through inductions of IgA GLTs (38), we tested the class-switch activity by quantifying the blood GLT levels. An increase of switch activity along the disease course was clearly demonstrated in the peripheral group but not in the central group (Figure (Figure6B;6B; p=0.013 vs. p=0.48, t-test, more details in Section “Materials and Methods”). The activity boost was compatible with the higher TGFβ1 levels for the peripheral group at dengue onset (Figure (Figure66A).
In this study, we applied a unified selection scheme among several data sets to find signatures of antibody repertoires. We took advantage of PCA to focus on a smaller set of contributing CDR-H3 clones and used SVC plus LR models to detect signature clusters out of those contributing clones. This pipeline successfully reproduced the infection signature of chronic hepatitis B carrier children (Figure S1B in Data Sheet S1 in Supplementary Material) (13) and identified previously reported IgG signatures of dengue fever in South America (Figures (Figures11 and and2)2) (11, 12). Although cluster ARQRGNWFDS was more prevalent in Mexico than in Nicaragua or Taiwan, cluster ARLDYYYYYGMDV was detectable in all three areas (Figures (Figures11 and and2).2). It is not clear why cluster distributions differ with geographic locations. Potential biological factors could include ethnic groups, viral strains, or even transmission vectors. However, the possible technical biases from processing pipelines could not be ruled out, either. The intrinsic stochastic nature of human antibody responses to viral infections could also play an important role. More signatures from a wider coverage of dengue-endemic areas might provide a comprehensible pattern to address the question.
IgA has been proposed as a diagnostic approach to dengue fever (39, 40). Recently, the relationships between IgA and hemorrhages have also been reported in several studies (41–44). However, it is still not clear what mechanism causes the bleeding phenotype. In this study, we found seven IgA clusters that were closely associated with hemorrhages (Figure (Figure5).5). Antigenic stimulations would likely account for the appearances of these CDR-H3 clusters. With recent advances in single-cell sequencing, it would be possible to get complete sequences of relevant monoclonal antibodies carrying these signatures, including the identities of light chains. Monoclonal antibodies accordingly reproduced would be of great use in addressing the underlying pathophysiology of dengue hemorrhages. A similar approach could be taken for other IgA-related bleeding disorders such as Henoch-Schönlein purpura or IgA nephropathy.
Innate IgA reactivity is an important contributing factor to shape IgA repertoires in mice (45). If the same is true for humans, the IgA specificity within human milk could also play a similar role. The corollary to the assumption would imply the compositions of IgA repertoires among dengue-infected infants would be affected by the prior IgA specificity within ingested milk. Thus, dengue-experienced mothers would influence their babies’ reactions to dengue viruses in an IgA-dependent manner. Existence of dengue-primed IgA in maternal milk could make these young kids prone to IgA-related bleeding disorders by differentially shaping their IgA repertoires. Such an effect would be absent from dengue-naive mothers. It would be intriguing to see if this hypothesis could address the unexplained dengue hemorrhages in young infants born to dengue-inflicted women. An easy prevention might be achievable by strict formula feeding for these young children.
In this study, we took healthy controls from another database to contrast immune repertoires from dengue patients. With this approach, we did reveal common infection signatures among dengue patients in both North or Central America and Taiwan (Figures (Figures11 and and2).2). It would be useful if we could expand the collections of immune repertoires from even more clinical scenarios. With more sequences available, we would expect to discover relevant information from disease in other categories as well. The identified sequences might play essentials roles under those circumstances. The reusability of these repertoire data could de facto expand the possibility of sequence-based medicine to a new horizon.
In sum, we developed a unified pipeline to extract signatures of antibody repertoires and used dengue fever as a model disease to gain insights in infection- and phenotype-specific contexts. The results well exemplified the stochastic nature of reacting IgG repertoires and uncovered cluster-yielding characters of IgA repertoires in dengue hemorrhages. Differentials of class-switch activity and isotype effectors in humoral immunity could be new factors to consider in viral infections. For inflammatory diseases without animal models or known pathogens, the framework could be valuable in offering perspectives not derivable elsewhere.
The protocol in recruiting dengue subjects was approved by the Institutional Review Board of Kaohsiung Medical University Hospital. Some of non-dengue controls were from a previous study, the protocol of which was approved by both the Institutional Review Board of National Health Research Institutes and the Institutional Review Board of En Chu Kong Hospital. Study methods were carried out in accordance with the relevant guidelines. Informed consents were obtained from all attendants.
C-HH, C-YL, H-CK, Y-JH, Y-WW, and J-YC helped in recruiting subjects. Y-HC and C-HY prepared sequencing libraries. W-HW did virus serotyping. S-FT, Y-HC, and H-HL designed the study. H-HL analyzed the data, wrote the manuscript and programmed the Python scripts.
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
We thank Dr. Anne Chao at Tsing Hua University for statistical consultations.
Funding. This work was supported by Ministry of Health and Welfare (106-0324-01-10-07) and Ministry of Science and Technology (104-2320-B-037-004).
The Supplementary Material for this article can be found online at http://www.frontiersin.org/article/10.3389/fimmu.2017.01726/full#supplementary-material.