Clonality Analysis of Immunoglobulin Gene Rearrangement by Next-Generation Sequencing in Endemic Burkitt Lymphoma Suggests Antigen Drive Activation of BCR as Opposed to Sporadic Burkitt Lymphoma

Patients

A total of 37 BL samples preserved in RNA later (RNA stabilization Reagent; Qiagen, Valencia, CA) were collected from the Department of Human Pathology of the Lacor Hospital, Uganda; from the Department of Medical Biotechnologies, University of Siena, Italy; and from “L.A Seragnoli” Haematopathology Unit, Policlinico S. Orsola, Bologna, Italy. Formalin-fixed, paraffin-embedded samples were available for all. Overall, 26 BL samples were collected from endemic areas in Uganda, Africa (eBL), while the other 11 BL cases were collected in Italy (sBL). Clinical epidemiologic details of the 37 BL cases are summarized in Table 1.

All cases were previously analyzed by gene expression profile (GEP) analysis and showed a signature consistent with molecular BL.^13,14

Total RNA was extracted from BL samples with Trizol, according to manufacturer’s procedure (Invitrogen, Life Technologies, Carlsbad, CA). The study was approved by each institutional ethical committee, and written permission and informed consent have been obtained before sample collection in accordance with the Declaration of Helsinki.

TCF3 and ID3 Mutation Assessment by RNAseq

The paired-end libraries (2 × 75 base pairs [bp]) were prepared according to the TruSeq RNA sample preparation v2 protocol (Illumina, San Diego, CA). The resulting libraries were sequenced on a HiScan SQ (Illumina) following the manufacturer's instructions. Briefly, 2 μg Poly(A)⁺ RNA was purified from total RNA using poly-T oligo–attached magnetic beads and then used for fragmentation into 130- to 290-bp fragments. The first strand of complementary DNA (cDNA) synthesis was performed using reverse transcriptase enzyme (SuperScript II; Invitrogen, Life Technologies) and random hexamer primer, followed by generation of double-stranded cDNA. AmpureXP beads (Beckman Coulter, Brea, CA) were used to purify the cDNA, and an end-repair step was performed to convert the overhangs, resulting from fragmentation, into blunt ends by 3′ to 5′ exonuclease activity. A single “A” nucleotide was added to the 3′ ends of the blunt fragments to prevent them from ligating to one another during the adapter ligation reaction. This approach was adopted to ensure a low rate of chimera (concatenated template) formation. Subsequently, sequencing adapters were added to the ends of the double-stranded (ds)–cDNA fragment, and a polymerase chain reaction (PCR) reaction was used to selectively enrich those ds-cDNA fragments that had adapter molecules on both ends, amplifying the amount of ds-cDNA in the final libraries. Last, PCR libraries products were purified by AmpureXP beads, and quality control analysis was assessed using a DNA-1000 (Agilent, Waldbronn, Germany). The quantification was performed by the Quant-it PicoGreen dsDNA Assay Kit according to the manufacturer’s protocol (Invitrogen, Life Technologies). Sequence variants were obtained using the SAVI (Statistical Algorithm for Variant Identification) algorithm independently for each sample. Candidate somatic mutations were obtained by eliminating common germline variants and recurrence. Sanger sequencing was used for validation and confirmation of TCF3 and ID3 single-nucleotide variants in sBL.^9-11,15,16

PCR Amplification and High-Throughput Sequencing by Roche 454 GS Junior Instrument

In total, 1 µg total extracted RNA per case was transcribed into cDNA using the QuantiTect Reverse Transcription Kit (Qiagen) according to the manufacturer’s instructions. Resulting cDNA was subjected to PCR amplification of the IGVH locus by using modified BIOMED-2 framework region 1 (FR1) consensus primers in conjunction with a consensus J segment primer.¹⁷

In this assay, each of the seven BIOMED-2 primers at the 5′ end contained an additional sequence, Multiplex Identifiers, that barcodes the sample Table 2. Thermal cycling was performed on an Applied Biosystems 2720 Thermal Cycler using the BIOMED-2 protocols. Small DNA fragments were removed using the AMPure PCR purification system (Agencourt, Beverly, MA) following the manufacturer's protocol (http://454.com/downloads/my454/documentation/gs-junior/method-manuals/454SeqSys_AmpliconLibraryPrepMethodManual_Jun2013.pdf). Amplicons were subsequently quantified using the Quant-iT PicoGreen dsDNA reagent (Invitrogen Corporation, Life Technologies). All amplicons were then pooled at an equimolar ratio. Subsequently, the sample pool was diluted to a final concentration of 1 × 10⁷ PCR fragment molecules/μL.

Table 2.

Modified FR1 BIOMED-2 Primers^a

Open in a new tab

NGS, next-generation sequencing.

^a

Forward and reverse primers are in blue, key sequences are in orange and yellow, and Multiplex Identifier (MID) sequences in green.

The amplicon-PCR–derived fragments were annealed to carrier beads and clonally amplified by emulsion PCR (emPCR). The emPCR was performed according to the manufacturer's protocol (http://454.com/downloads/my454/documentation/gs-junior/method-manuals/GSJunioremPCRAmplificationMethodManualLib-A_March2012.pdf). The beads were isolated and compartmentalized into droplets of an aqueous PCR buffer in oil emulsion. Subsequently, the emulsions were broken by isopropanol to facilitate collection of the amplified fragments bound to their specific beads. The beads carrying single-stranded DNA templates were enriched, counted, and deposited into the picotiter plate for sequencing (http://454.com/downloads/my454/documentation/gs-junior/method-manuals/GSJuniorSequencingManual_Jan2013.pdf).

Data analysis was performed using the Roche (Basel, Switzerland) proprietary software package for the GS Junior system. Image acquisition, image processing, and signal processing were performed during the run.

An unmutated chronic lymphocytic leukemia (CLL) case, BL cell lines, and reactive lymph node sequencing were applied as control samples for monoclonality and polyclonality.

Bioinformatical Analysis

By using the GS Junior system, we obtained, for each sample, an average depth of 3,700 reads with a minimum depth of 772 reads and a maximum depth of 7,869 reads.

The pooled raw sequences obtained by GS Junior sequencing were split by using the 454 seqTool script and custom script in Perl language according to different sample-specific tag sequences.

To reduce the single-nucleotide variant complexity, each single subset of sequences was clustered by cap3 software (http://bioweb.pasteur.fr/seqanal/interfaces/cap3.html) with more than 98% overlap identity by obtaining a list of consensus sequencing, and therefore the sequencing reads that differed by up to two single-nucleotide substitutions or gaps were aggregated.¹⁸ By using these stringent criteria, artifactual clonotypes were minimized. For each consensus sequence corresponding to different intrasample clones, we calculated the number of raw sequences used to indicate the abundance of each different rearrangement present in the sample. Finally, the consensus sequencing was multialigned by the Clustalw (http://www.ebi.ac.uk/tools/clustalw2) program with default parameters.¹⁹

Characteristic 454 sequencing errors, such as insertions/deletions near homopolymer tracts, were also identified, and sequences containing such errors were also discarded. Low-frequency clonotypes that differed by up to 1% from another clonotype with significantly higher frequency were deemed to be the result of technical artifact (PCR or sequencing artifacts), and the sequence differences were corrected for clonotype aggregation.

Validation of NGS Data by Sanger Sequencing and Subcloning

To validate the reliability of NGS data, 26 eBL and 11 sBL samples were amplified in two separated reactions. Both reactions were carried using the FR1 VH-JH BIOMED-2 primers as previously described.¹⁵

Purified PCR amplicons were subjected to direct sequencing and cloning into the pCRII-TOPO vector (TOPO TA Cloning Kit; Invitrogen, Paisley, UK) and were used to transform TOP 10 cells (Invitrogen) according to the manufacturer’s instructions. At least 20 clones per case were sequenced in both directions with an ABI PRISM 310 Genetic Analyzer (Applied Biosystems, Weiterstadt, Germany) using Sp6 reverse and T7 forward primers, respectively.

The sequences obtained from direct sequencing and cloning were analyzed with the IMGT/V-QUEST (http://www.imgt.org/IMGT_vquest/vquest?livret=0&Option=humanIg; Marie-Paule Lefranc [Marie-Paule.Lefranc@igh.cnrs.fr], Université de Montpellier, CNRS, LIGM, IGH, SFR, Montpellier, France).^20,21

Sequence Data Analysis

All the obtained sequences were aligned with IMGT tools for nucleotide analysis of immunoglobulin (IG) V-(D)-J repertoires, polymorphisms, and IG mutations.²¹

In NGS data, the following two parameters have been analyzed for each V-D-J combination: (1) the number of B-cell clusters, defined as sets of sequences having the same unique V-D-J and heavy complementarity-determining region 3 (VH CDR3) rearrangement (number of clusters), and (2) the total number of sequences obtained in this V-D-J combination (number of reads). Clusters composed of more than 50 reads and showing the same VDJ and CDR3 but variation in the VH domain were considered subclusters, indicating intraclonal diversity.

The following additional information was also extracted: (1) IGHV gene usage and percentage of homology to germline considering the 2% cutoff to discriminate between mutated and unmutated immunoglobulin genes; (2) the number of nucleotide mutations associated with amino acid changes (replacement mutation) and the number of nucleotide mutations not associated with amino acid changes (silent mutation) of the rearranged immunoglobulin genes^22-24; and (3) the amino acid changes, classified as very similar (+ + +), similar (+ + –, + – +), dissimilar (– +, – + –, + –), and very dissimilar (– – –) according to hydropathy, volume, and IMGT physicochemical classes (http:/www.imgt.org/IMGT_vquest/vquest?livret=0&Option=humanIg).^21,25 These amino acid changes are useful to understand affinity maturation mechanisms, which result in extensive alteration of the structural protein conformation that is important in determining the binding to the antigens.

All “nonubiquitous” sequence changes from the germline were evaluated and further characterized as follows: (1) unconfirmed mutation (a mutation observed in only one subclusterd sequence from the same specimen [“unique”]) and (2) confirmed mutation (a mutation observed more than once among subcluster sequences from the same specimen [“partially shared”]). Where a clonal expansion was identified, the sequences that were more mutated were chosen as representative of the cluster to examine the characteristics of the mutations.

Statistical Analysis

Descriptive statistics for discrete parameters included counts and frequency distributions. The significance of the bivariate relationship was assessed by the Fisher exact test. For all comparisons, a significance level of P at .05 was set, and all statistical analyses were performed using GraphPad InStat version 3.05 (GraphPad Software, San Diego, CA).

TCF3 and ID3 Are Mutated in eBL

ID3 and TCF3 single-nucleotide variants (SNVs) were found in few eBL cases. Specifically, ID3 SNVs were found in six (23%) of 26 eBL cases, while TCF3 was found in three (11.5%) of 26 eBL cases. Since one sample (BL23) was found to carry SNVs in both genes, eight (30.7%) of 26 eBL cases had ID3 and/or TCF3 variants Table 3. Sanger resequencing of genomic DNA confirmed ID3 and TCF3 mutations in 30% and 64% of eBL and sBL samples, respectively. In three samples, matched DNA obtained from peripheral blood was available; therefore, we could assess the absence of the SNVs, indicating that the recognized variants corresponded to somatic mutations Figure 1.

Comparison of Sequencing Data Derived From the NGS Analysis With Sanger Sequencing

Multiplex FR1 immunoglobulin H PCR products were generated and subjected to clonality analysis via 454 NGS, followed by data processing, and compared with capillary electrophoresis.

A total of 130,000 sequences/reads were generated with an average depth of 3,700 reads, with a minimum depth of 772 reads and a maximum depth of 7,869 reads for each sample. All 37 cases were clonal on NGS using the criterion that a clonal cluster(s) must be at least fourfold more abundant than the next largest clonotype of the background.²⁶ In particular, the presumed monoclonal clusters ranged from 554 to 6,780 reads, representing from 61.8% to 99.7% of the total reads. On the other hand, the polyclonal background was composed by high percentages of single reads, ranging from 0.01% to 2.04% of the total reads.

When all the sequences were aligned with IMGT tools for nucleotide analysis of immunoglobulin (IG), polymorphisms, and IG mutations,^21,24 clusters showing identical IGHV, IGHD, and IGHJ usage and CDR3 regions as the presumed monoclonal clusters were detected, ranging from 63 to 470 reads, with a range from 8.6% to 35.5% of the total reads.

Distinct clusters showed specific mutational profiles. They exhibited shared as well unique somatic mutations and sequences with different lengths for the presence of nucleotide deletion. This finding is indicative of intraclonal diversity and points to the emergence from a common ancestor. This phenomenon was particularly evident in all eBL cases, while only one case of sBL showed different subclusters without differences in length but only differences in point mutations Table 4 and Figure 2.

Figure 3 illustrates an example of the branching of the lymphoma clone (case eBL31) and shows that distinct subclones evolved along similar, although separate pathways, suggesting a selective pressure occurring in parallel in distinct subclones, thereby improving their BCR affinities.

To confirm that the presumed monoclonal clusters detected by NGS represented the major tumor cell population, we carried out Sanger sequencing of PCR products in 37 BL cases (26 eBL and 11 sBL samples) and analyzed the results. The VDJ identity of the presumed monoclonal cluster detected by NGS always corresponded to that inferred with Sanger direct sequencing.

In addition, we performed cloning of the PCR products and sequenced at least 20 clones per case. With this method, we were able to identify at least 69 different clones with the same VDJ and CDR3 region out of 150 detected by NGS (Table 4). This result can obviously be explained by the much larger amount of sequence data generated by NGS. In fact, cloning of the PCR products confirmed only the clones composed by the highest number of reads in NGS.

Targeting of Somatic Hypermutation

All eBL cases carried IGHV genes of highly mutated status with a germline identity ranging from 79.81% to 97.26% (median, 90.54%); in sBL cases, the mutational variation compared with the germline ranged from 90.88% to 97.42% (median, 94.6%). The eBL specimens demonstrated the presence of somatic mutations in the range of five to 42 mutations. In contrast, in sBL, the range of mutations was from three to 19 (Table 4). Evaluation of the range of somatic mutations disclosed a significant difference between the two subgroups (eBL and sBL) (P< .0001).

Total mutations in eBL cases were 583 (mean, 22 mutations per case), with a significant prevalence of replacement mutations (417/583), which were more than twice that of the silent mutations (166/583) (P < .0001). The total number of mutations in sBL was 100 (mean, nine mutations/case), with a slight prevalence of replacement mutations (82/100) over silent mutations (18/100) (Table 4).

In most of the eBL cases, the majority of mutations occurred in FRs rather than CDRs, leading to a complementary-determining region (CDR)/FR ratio greater than 0.3, which would indicate antigen selection. In addition, significant accumulation of replacement as opposed to silent mutations was found in the CDRs, which form the antigen binding sites.²⁷

Our results showed that replacement mutations most often determine substitutions of amino acids with dissimilar/very dissimilar physicochemical properties rather than similar/very similar properties, according to the standardized IMGT amino acid classes (Table 4).^20,25 Based on these results, it could be suggested that somatic hypermutation in BL induces conformational changes that allow the spatial structure of the antigen binding site to better contact antigens, thus improving affinity.²⁸

The pattern of amino acid substitutions was similar in both eBL and sBL, but in sBL, the low number of replacement mutations does not allow any significant conclusion.

Analysis of IghV Gene Usage and CDR3 Features

Clonal and in-frame VH gene sequences from BL cases derived from VH1, VH3, and VH4 families. In particular, in eBL, the comparison between the clonal VH sequences and the known germline revealed that 12 (46%) of 26 sequences were most homologous to the VH3 family (VH3-7, VH3-9, VH3-11, VH3-15, VH3-23, VH3-30, VH3-72, VH3-74) and, to a lesser extent, to the VH1 family (10/26, 38%) (VH1-2, VH1-3, VH1-8, VH1-69). The remaining four were most closely related to the VH4 family (VH4-34). In sBL, the usage of the three families was similar to that found in eBL (VH1 was found in one of 11 cases, VH3 in eight of 11 cases, and VH4 in two of 11 cases).

These results enabled us to test whether the cells from BL were derived randomly from the V gene repertoire expressed in normal peripheral blood lymphocytes or whether they exhibited a biased V gene usage. To this end, we compared the VH gene usage of this B-cell population with the rearrangements of normal peripheral B lymphocytes (PBLs).²⁹ Our results showed that in BL cases, the VH1,VH3, and VH4 family genes were significantly overrepresented compared with PBLs (P < .05). In addition, when we looked at the individual genes, we found a preferential expression of VH3-23 (9/37, 24%), the more represented region in memory B cells, and VH1-69 (6/37, 16%) and VH4-34 (6/37, 16%). This IGHV repertoire usage is remarkably biased in comparison to the IGHV repertoire of normal B cells, where VH3-23, VH1-69, and VH4-34 are used by 12.1%, 1.4%, and 3.9%, respectively (IGHV repertoire in eBL and sBL, respectively) Figure 4.²⁹

CDR3 from all cases was analyzed both at the nucleotide and amino acid levels. No identical or similar sequences were identified when comparing the sequences from one case with those of any other case.

The average VH CDR3 length of BL cells was similar to that of germinal center and marginal zone B cells (amino acid, 14.16 vs 14.75 and 15.01, respectively).

We determined the IGHD and IGHJ genes used by the IGHV-D-J sequences analyzed. The IGHD genes used were IGHD3 (48%) and IGHD2 (16%) and several others to a much lower extent, including IGHD1, IGHD4, IGHD5, and IGHD6. In particular, we demonstrated a bias toward IGHD3-10 (25%), which was considered statistically significant compared with the PBLs (P < .05).²⁹,³⁰ The IGHJ gene usage in all cases showed that IGHJ4 (56%) was used the most, followed by IGHJ3 (16%), IGHJ6 (16%), IGHJ5 (8%), and IGHJ1 and IGHJ2 (4%). In sBL, no significant difference in the distribution between the D and J segment usage has been found.