Sixteen individuals, enrolled in a protocol approved by the Institutional Review Boards of the University of Florida and the University of South Florida, included four groups (n = 4 per group): healthy controls (HC), subjects with SLE, and HIV-1 infected individuals either therapy-naïve (HIV) or receiving combination antiretroviral therapy (cART) (HIVTx) (Table 1). Groups within the cohort were age-balanced with median (25–75% quartile range) age of 20 (16–22) years. SLE subjects were untreated and diagnosed by clinical and laboratory criteria defined by the American College of Rheumatology (42). Therapy-naïve subjects were HIV-1 infected by sexual transmission for at least 6 months (Table 1). HIVTx individuals were infected either through maternal to child transmission or by contaminated blood products for 16.0 (12.3–21.3) years, treated for 8.0 (1.5–17.5) years, and achieved viral suppression (log₁₀ HIV-1 RNA copies/ml < 1.7), and CD4 > 25% for 1.9 (0.5–6.0) years at the time of study (Table 1). No subject received any vaccination 30 days prior to study entry, and/or had acute infection to reduce the chance of plasma cell circulating in peripheral blood. Informed consent was obtained from all subjects.

Immunofluorescence staining was performed using the whole blood lysis method (43) and analysis by LSR2 flow cytometer (BD Biosciences, Franklin Lakes, NJ, USA) with Diva software (BD Biosciences, San Jose, CA, USA). Monoclonal antibody panel included pan B cell marker PE Cy7-conjugated anti-CD19 (BD Biosciences), memory B cell marker Qdot 655-conjugated anti-CD27 (Invitrogen, Carlsbad, CA, USA), and surface IgM with APC-conjugated anti-IgM (BD Biosciences). B cell percentages ranged from 6.6 to 16.1% of peripheral blood mononuclear cells (PBMC) across the groups, CD19⁺ IgM B cells ranged from 70 to 90% of total B cells (Table 1).

Analysis of the IGHM CDRH3 transcriptome was performed using total peripheral blood mRNA without isolation of IgM B cells. Specificity for the IgM transcriptome was achieved by reverse amplification primers Cμ15, Cμ5, and Cμ2 homologous to the IgM constant region (44) (Figure 1A). An amplicon library was constructed for each subject by RT-PCR using SuperScript™ One-Step RT-PCR with Platinum Taq (Invitrogen) and GoTaq colorless Master Mix (Promega, Madison, WI, USA) from ∼200 ng mRNA (equivalent to about 100,000 B cells or 80,000 IgM B cells) extracted from total PBMC by MicroPoly[A]Purist™ kit (Ambion, Austin, TX, USA).

Forward inner primer pan V_H-FR3 is located in framework 3 and does not overlap CDRH3 (44). Using the software Oligo (Molecular Biology Insights, Inc., Cascade, CO, USA), pan V_H-FR3 primer was predicted to bind with similar capacity to a unique sequence in each IGHV family when evaluated against IGHV family specific reference sequences. In contrast, the pan V_H-FR1 primer (44) displayed different binding capacity to different IGHV families. Specifically, pan V_H-FR1 primer bound preferentially to IGHV3; displayed low binding capacity with IGHV4, IGHV5, and IGHV7; failed to bind to IGHV2 or IGHV6; and bound to false priming sites in IGHV1, IGHV4, and IGHV5. To overcome the different binding capacity of the published primer sequence, we designed a modified panV_H-FR1* primer (CAGGTGCAGCTGGAGCAGTCTGG_) that was one nucleotide shorter and substituted A and C for T and G in the published primer sequence, respectively. When evaluated by Oligo, the modified panV_H-FR1* primer displayed similar binding capacity to each IGHV family reference sequence and no binding to false priming sites. Forward and reverse primers, pan V_H-FR3, and Cμ2, in the final amplification were conjugated at 5′ ends with A (GCCTCCCTCGCGCCATCAG) or B (GCCTTGCCAGCCCGCTCAG) adaptor, respectively, to generate amplicons ranging from 150 to 200 nucleotides. Gel-purified amplicons were submitted to the Interdisciplinary Center for Biotechnology Research (University of Florida) for 454-pyrosequencing using a Genome Sequencer FLX (454 Life Sciences) according to the manufacturer’s protocol.

A bioinformatics pipeline was developed to facilitate analysis of large numbers of relatively short IGHM CDRH3 sequences that could not be processed by conventional IMGT/V-QUEST analysis (Figure 1B) (45). Raw reads ranged from 4,500 to 21,000 pyrosequences per subject. A quality control step filtered 9–21% low quality reads with ambiguous nucleotides, more than one error in either primer tag, or failure to align to reference sequences in the germ line IGH V-D-J recombination amino acid reference sequence database (see below), leaving 4,000–17,000 quality sequences per sample with no significant difference in number of sequences among the groups (Table 2).

To overcome the limitation of IMGT/V-QUEST for high throughput classification of relatively short sequences, a novel custom reference sequence database containing over 700,000 germ line IGH V-D-J amino acid (aa) sequences was established by generating in silico all possible combinations of germ line IGHV (285 unique aa sequences), IGHD [44 unique nucleotide sequences translated in 6 reading frames (RF) producing 222 aa sequences], and IGHJ (13 unique aa sequences) alleles downloaded from IMGT (45). Each reference sequence was annotated for IGHV, IGHD, and IGHJ allele use, IGHD RF in the header with C104 and W118, the 5′ and 3′ boarders of CDRH3, marked by aligning to germline IGHV and germline IGHJ reference sequences using FASTA (http://fasta.bioch.virginia.edu/fasta_www2/fasta_down.shtml) (46). After removing adapter sequences, query sequences were then aligned to the reference sequences using FASTY. Use of IGHV, IGHD, and IGHJ alleles, and IGHD RF of each query 454-sequence was extrapolated from the header of the best-matched reference sequence, the C104 and W118 positions of each query sequence were identified, the part of the query sequence matched to the C104-W118 region of the reference sequence was extracted from the entire sequence, and CDRH3 charge and hydropathy index calculated using an in-house code (47, 48). IMGT/Junction Analysis was then performed allowing 3D-GENEs (49). Using this strategy, the frequency distribution of IGHD and IGHJ alleles, with IGHD2, IGHD3 and IGHD6, and IGHJ4 as the dominant genes, was similar to previous reports (50, 51). However, due to the relatively short IGHV sequence, we observed some ambiguity in IGHV4 allele assignment; consequently, no comparisons of IGHV alleles without or with SHM were made.

The total number of mutations including silent- and non-silent-mutations in CDRH3 between C104 and W118 (40) was calculated to summarize extent of SHM as: (ΣSHM within CDRH3 in each SHM⁺ sequence ÷ ΣCDRH3 nucleotide length in each SHM⁺ sequence) × 100 nucleotides. The effect of N and P nucleotides to the SHM was similar among the study groups because the number of N and P nucleotides was minimized when the setting for D-GENES was 3, and the frequency of sequences with one D, D-D fusion and D-D-D fusion was similar among study groups (data not shown).

To avoid potential ambiguity in scoring SHM, an average of 14.3 ± 0.6% of reads among groups with a single nucleotide difference from aligned reference sequences was removed from analysis (Table 2) (40). IGHM sequences with two or more mutations in CDRH3 were classified as SHM⁺ representing IgM memory B cells, while SHM⁻ sequences represent naïve B cells. After correction of deletion and insertion, the PCR and pyrosequencing-induced rate of misincorporation tested in our control analysis of clones of sequences was 0.18 errors per 200 nucleotides, the longest sequence length, similar to other reports (52–54) and well below the level of SHM identified in the IgM populations.

ESPRIT was applied to study the biodiversity (genetic complexity) of nucleotide sequences as well as V-D-J combinations of the IGHM CDRH3 transcriptome repertoire by clustering the sequences at 0% genetic distance (55). Rarefaction analysis measures increase in biodiversity along the depth of sequencing (number of sequences). The deeper the initial slope is, the higher the biodiversity. Left shift or right shift of the curve indicates an increase or decrease of biodiversity. Chao1 analysis inferred maximum biodiversity within the input templates (55, 56). Biodiversity is influenced by input cell number, the coverage (sequence number/input cell number), and clonality (frequency of clusters with more than 10 repeated sequences). Coverage was ∼10% in each individual to minimize the influence of preferential amplification of replicate templates. Biodiversity was weighted by the absolute number of input IgM B cells (calculated from the absolute lymphocyte count multiplied by the percentage of CD19⁺ IgM B cells) to make the data comparable among study groups. To avoid a potential bias from sequence number, clonality was averaged on the same number of sequences randomly drawn 1,000 times from each sequence pool.

To rule out an influence by differences in sequence numbers among individuals, a bootstrap resampling method was used to generate 1,000 rarefaction curves and Chao1 values from random data sets with the same number of sequences for each individual. Calculated and estimated biodiversity derived from averages of 1,000 rarefaction curves produced similar findings (result not shown).

Comparisons among study groups were performed by ANOVA. Kolmogorov–Smirnov test was used to test Gaussian distribution. Overall difference of distribution of CDRH3 length frequencies, and charge and hydropathy along CDRH3 lengths among study groups were evaluated as described (57). Comparison of frequency of sequences in each gene allele, or each CDRH3 length between SHM⁺ and SHM⁻ sequences was performed by paired t-test. Pearson correlation evaluated the relationship between frequency of SHM⁺ sequences and frequency of peripheral CD27⁺ IgM B cells. A primary purpose of this pilot study with relatively small sample size was to detect meaningful trends, thus parametric methods instead of conservative non-parametric methods were performed. P values were unadjusted for multiple tests to increase sensitivity to identify differences. Statistical analyses were performed using SAS version 9.1 (SAS Institute, Cary, NC, USA) with p < 0.05 (two sided) defined as significant.

IGHM CDRH3 sequences were initially classified by absence or presence of SHM as a molecular means to distinguish naïve B cells from IgM memory B cells. Frequency of sequences with SHM in HIV or HIVTx groups was similar to HC group, but significantly reduced in individuals with SLE (p = 0.02) (Figure 2). In contrast, extent of SHM (mutations/100 nucleotides) in IGHM CDRH3 among groups was similar. While the proportion of B cells expressing CD27⁺ IgM memory phenotype was similar among groups, percentage of phenotypically mature CD27⁺ IgM B cells was significantly less than the frequency of IgM sequences with SHM within each group (Figure 2). No correlation occurred between the proportion of CD27⁺ B cells and SHM (r² = 0.06, p = 0.35). Subsequent analysis of the IgM transcriptome was based on molecular assessments of sequences without or with SHM in IGHM CDRH3.

Deep sequencing data sets support application of novel rarefaction analysis to Ig biodiversity that cannot be inferred from analysis of limited numbers of sequences. Rarefaction curves provide comparison of biodiversity along with the depth of sequencing among study groups; deeper slopes indicate greater biodiversity. IGHM CDRH3 among healthy young adults displayed a range of biodiversity that overall was greater than biodiversity in SLE or HIV-1 infection (Figures 3A,B). Differences in biodiversity between health and disease were apparent in populations of IgM sequences without SHM, and to a greater extent among sequences with SHM due to contribution of SHM. At the depth of sequencing, none of the rarefaction curves reached saturation, indicating a preponderance of unique sequences and an extent of IGHM CDRH3 diversity that exceeded the depth of sequencing.

When maximum IGHM CDRH3 biodiversity was inferred, the SHM-negative repertoire within healthy individuals displayed significantly greater maximum biodiversity than SLE (p = 0.031), HIV (p = 0.006), or HIVTx (p = 0.020) (Figure 3C). Likewise, the repertoire with SHM among healthy individuals also displayed significantly greater biodiversity than SLE (p = 0.021), HIV (p = 0.004), or HIVTx (p = 0.010) (Figure 3D). ART failed to normalize biodiversity in IGHM CDRH3 repertoire either without or with SHM in HIV-infected subjects (Figures 3C,D).

Multiple factors contribute to biodiversity. SHM is one factor, but extent of SHM was similar among study groups. Likewise, extent and distribution of clonality were similar among sequences without or with SHM in each study group with a frequency of ∼99% for unique clusters about 0.2% of clusters with more than 10 repeated sequences. Consequently, to identify molecular differences that might contribute to reduced biodiversity without or with SHM among SLE or HIV-infected groups, CDRH3 length variation, charge, and hydropathy distribution, as well as IGHD and IGHJ allelic use, and diversity of allele combinations were investigated.

In all groups, CDRH3 length variation displayed Gaussian distributions (Figure 4). CDRH3 regions without SHM ranged from 4 to 31 amino acid residues with a peak frequency of 15 amino acids (Figure 4A). SHM⁺ CDRH3 regions ranged from 5 to 34 amino acids, with lengths of 14 amino acids occurring most frequently (Figure 4B). Overall length distribution differed significantly between sequences without or with SHM (p = 0.004) (Figure 4C). In general, CDRH3 lengths of 10–17 amino acids occurred with greater frequency among sequences with SHM than in sequences without SHM (p < 0.0001). The frequency of long CDRH3 regions (27–34 amino acids) ranged from 0.03 (± 0.06%) to 0.13% (± 0.12%) in sequences without SHM and 0.13 (± 0.04%) to 0.16% (± 0.10%) in sequences with SHM with no significant differences among study groups (Figure 5). Charge and hydropathy distribution across CDRH3 amino acid length range were similar among all groups indicating that junctional modifications in sequences without or with SHM were unchanged (Figure 6).

Nearly 90% of the 44 IGHD alleles were detected across the study groups (Figures 7A,B). Sequences with SHM showed significant increases in use of multiple IGHD alleles compared with sequences lacking SHM. Similarly, over 75% of 13 IGHJ alleles, predominantly in IGHJ3 and IGHJ4 families, were identified independent of SHM (Figures 7C,D). A significant increase in three IGHJ alleles was detected in sequences with SHM (Figure 7D). While all study groups had similar IGHD and IGHJ allele usage, the major difference in allele usage was between sequences with or without SHM.

Biodiversity resulting from different combinations of alleles comprising the IGHM CDRH3 regions was assessed in health and disease groups. Every sequence was composed of V-D-J gene segments suggesting intact recombination machinery in all individuals independent of health or disease. When comparing numbers of different combinations along the depth of sequencing using rarefaction curves, a left shift was observed in healthy individuals, indicating greater diversification of combinations whether or not sequences contained SHM (Figures 8A,B). The difference between health and disease was more pronounced when maximum biodiversity was inferred by Chao1 algorithm based on both the richness and evenness of different V-D-J combinations (Figures 8C,D). ART failed to restore biodiversity in either group of IgM transcripts. Overall, reduced diversity of V-D-J combinations was a major contributor to the difference in biodiversity within IGHM CDRH3 sequences between groups.

Combinations of alleles within IgM sequences that might be similar to IgG bn-HIV-Abs were evaluated, but undetected among any individuals (Table 3) (23, 26, 27, 58–68). In contrast, combinations of families found in many bn-HIV-Abs were identified among IgM sequences across the four study groups (Table 4). While few, if any, family combinations similar to M66.6/MPER or CD4 binding site antibodies b12 and VRC-PG04 were detected, combinations similar to other bn-HIV-Abs were identified among IgM sequences in a majority of subjects. In general, family combinations were more frequently represented among IgM sequences with SHM, although frequencies greater than 0.2% appeared in sequences with or without SHM. Frequency of family combinations was not a function of the targets in Envelope, as combinations directed toward CD4 binding site, N-linked glycans, V2/V3, or MPER were detected in almost all individuals.