Characterizing adjuvants' effects at murine immunoglobulin repertoire level

Abstract

Generating large-scale, high-fidelity sequencing data is challenging and, furthermore, not much has been done to characterize adjuvants' effects at the repertoire level. Thus, we introduced an IgSeq pipeline that standardized library prep protocols and data analysis functions for accurate repertoire profiling. We then studied systemically effects of CpG and Alum on the Ig heavy chain repertoire using the ovalbumin (OVA) murine model. Ig repertoires of different tissues (spleen and bone marrow) and isotypes (IgG and IgM) were examined and compared in IGHV mutation, gene usage, CDR3 length, clonal diversity, and clonal selection. We found Ig repertoires of different compartments exhibited distinguishable profiles at the non-immunized steady state, and distinctions became more pronounced upon adjuvanted immunizations. Notably, Alum and CpG effects exhibited different tissue- and isotype-preferences. The former led to increased diversity of abundant clones in bone marrow, and the latter promoted the selection of IgG clones in both tissues.

Keywords: Genomic analysis; Immunology; Machine learning; Sequence analysis.

Graphical abstract

Figure 1

Sequencing library preparation protocol and data preprocessing pipeline We standardized the library preparation and data preprocessing pipeline by incorporating features and steps that had been proved to improve the quality and accuracy of RNA sequencing. (A) Each group of three mice received four doses of immunizations (PBS, OVA, OVA+CpG and OVA+Alum). One week after the fourth injection, spleen and bone marrow were harvested. Splenic B cells and bone marrow cells were prepared and total RNA was isolated. (B) The 5′RACE library prep protocol started with the TSOs containing 5′-biotin blockers and UMIs including spacers and RT-PCR cycles running at a higher temperature. The target Ig amplification PCR was achieved with a reverse isotype-specific primer containing short random nucleotide to increase the sequence diversity in the beginning of Read 2 sequencing. Both PCRs were monitored at real time to avoid over-amplification. Pair end sequencing was done on an Illumina HiSeq2500 system. (C) The data preprocessing pipeline does not require pair-read joining, and in all steps separated reads were used as input. We also developed and implemented new algorithm to de-duplicate UMI group reads and generate consensus sequences. Ig heavy chain VDJ annotation and clonal analysis were carried out using the Cloanalyst software. (D) The sequence logo plot for the TSO with UMIs. The UMI is different from a typical one by including spacers with reduced degeneracy. (E) Sequence yield of each sample at different steps between the pair-end joining protocol and non-joining pipeline in this work. Open bars are for the sequence yields of the typical pair-joining protocol and filled bars are for the non-joining separated reads pipeline. White bars are yields in number of sequences and orange yields in percentages. In the end, with the non-joining pipeline on average we obtained about 10-fold more sequences and 8-fold more clones for each sample. (Data are represented as mean ± SEM).

Figure 2

IGHV mean mutation frequency is affected by OVA immunizations and adjuvants in a tissue- and isotype-dependent fashion IGHV annotation was done by aligning to germline sequences with the software tool, Cloanalyst. Mutations were identified and expressed as per-nucleotide mutation frequency. The distributions of IGHV mutations and means and standard deviations were obtained and used for statistical tests. (A) IGHV mutation frequencies in different isotype compartments in non-immunized mice (PBS controls). (B) IGHV mutation frequencies in different tissue in non-immunized mice (PBS controls). (C) IGHV mutation frequencies were affected by OVA immunization and the effects were isotype- and tissue-dependent. (D) Mean mutation frequency difference between bone marrow and spleen. The difference was calculated as the BM mutation frequency subtracted by that of spleen for each OVA immunization group of different isotypes. A positive value means the BM compartment has a higher mutation frequency than the respective spleen compartment, a negative value means higher mutation frequency in the spleen, and zero means identical mutation frequency between two tissues. (Data are represented as mean ± SEM. Statistical p values were corrected by FDR and significant p values were indicated in figure).

Figure 3

OVA immunization affects IGHV gene usage and adjuvants modify the tissue distribution of the effects IGHV genes were annotated and gene usage was expressed as the percentage of each IGHV gene segment in each sample. The data were transformed into compositional data format. PCA was applied to reduce data dimensions and group correlating gene usages. Statistical analyses were done on top PCs with large contributions. (A) The PCA scree plot showed the individual and accumulative PC contributions. The top six PCs with largest contributions (each explaining >5% total variation) have an accumulative contribution of about 58%. We also plotted the Cronbach Alpha value of each PC (gray dashed line), which is a measure of internal PC consistency. (B) The PC loading plot. It showed correlations of IGHV gene usages projected on the top two PC plane. Genes are highly correlated if they point to the same directly, and contribute more to a PC component if they have larger loadings. (C) Sample distribution on the first two PC space. The samples were labeled by their isotype (point shape) and immunization group (color). (D) The PC1 distributions by different isotype, tissue and immunization status. PC1 explained about 15% of total variations and showed a significant immunization effect in a tissue- and isotype-dependent way. (E) The raw IGHV gene usage patterns to verify the PC1 distribution. We plotted the top 10 raw gene usages (centered by the mean of each group) with largest contributions to PC1. Both positively (top) or negatively (bottom) contributing gene usages were shown. (F) The PC3 distributions by different isotype, tissue and immunization status. PC3 explained about 9.5% of total variations and showed a significantly different pattern induced by the OVA+CpG immunization in the spleen IgG compartment. (G) The raw IGHV gene usage patterns to verify the PC3 distributions. We plotted the top 10 raw gene usages (centered by the mean of each group) with largest contributions to PC3. Both positively (top) or negatively (bottom) contributing gene usages were shown. (Data are represented as mean ± SEM. Statistical p values were corrected by FDR and significant p values were indicated in figure).

Figure 4

CDR3 length is affected by OVA adjuvant immunizations in a tissue- and isotype-dependent fashion IGHV genes were annotated and CDR3 regions was identified. The CDR3 length distribution of each sample was determined, and the mean CDR3 lengths were compared and plotted by different compartments and immunization status. (Data are represented as mean ± SEM. Statistical p values were corrected by FDR and significant p values were indicated).

Figure 5

Clonal diversity as Hill numbers Ig sequence data were annotated and partitioned into clones. Diversity profiles as Hill numbers were estimated based on clone abundance data with orders between 0 and 4. They were also scaled to control the effect of clone richness and focus on clone evenness (relative clone abundance/size). (A) Diversity profiles of different isotype compartments under the OVA+Alum immunization. Each line is for one isotype compartment in different tissue. Standard errors are plotted for diversities at the order of 0, 1, 2, 3 and 4. (B) Example clone size distribution pie charts. Here, it shows the data for one mouse immunized with OVA+Alum, and results for four compartments of different isotype and tissue combinations were shown. (C) Scaled diversity profiles of different isotype compartment at the steady state. Each line is for one isotype compartment in different tissues. Standard errors are plotted for diversities at the order of 0, 1, 2, 3 and 4. (D) Scaled diversity profiles of different tissue samples at the steady state. Each line is for one tissue sample in different isotype compartments. Standard errors are plotted for diversities at the order of 0, 1, 2, 3 and 4. (E) Scaled diversity profiles of samples under different immunization status. Each line is for one samples under different immunization in different tissue-isotype compartments. Standard errors are plotted for diversities at the order of 0, 1, 2, 3 and 4. (F) Scaled diversity distribution and comparison of diversity values at the order q = 4 for samples under different immunization conditions in different tissue-isotype compartments. (Data are represented as mean ± SEM. Statistical p values were corrected by FDR and significant p values were indicated. ∗, p < 0.05; ∗∗, p < 0.01, ∗∗∗, p < 0.001, ∗∗∗∗, p < 0.0001).

Figure 6

Evaluation of clonal selections and their changes by OVA immunization with adjuvants We developed a new metric to quantify the changes of clone phylogenetic trees to reflect the effect of selection. It involves two measurements, mean clonal mutation frequency and pairwise intra-clonal dissimilarity. Under no selection, a clone evolves to increase these two parameters in proportion, and two parameters show as a straight line on an x-y 2D plot. Otherwise, some parts of a clone tree are selected, and the linear relationship between the two parameter is distorted. To analyze and compare the selection strength metric, we plotted the two parameters of the largest 40 clones of each sample (A) and defined a threshold (the green dashed line), above which are clones showing some level of selection. For each sample, on the basis of such threshold we can obtained the percentage of selected clones among the largest clones. The percentages were compared statistically and results were shown in (B). (Statistical p values were corrected by FDR. Significant p values are shown for comparing OVA immunization effects).

Figure 7

Unsupervised machine learning to characterize Ig repertoire changes Multiple factor analysis (MFA) was carried out on the groups of variables including IGHV gene usage, mutation frequency, CDR3 length, clone diversity (expansion) and clone selection level. Sample data were then projected on the first two dimensions of MFA arranged by their isotype (A), tissue location (B), and immunization conditions (C). MFA can also summarize and project supplementary category variables (isotype, tissue and immunization) into the new PC dimensions (D) to show their correlations. Furthermore, the correlations between the variable groups and the first PC dimensions were shown in (E). Finally, the partial point graphs of supplementary category variables in MFA were shown in (F), in which supplementary categorical variables (such as isotype, tissue and immunization) were plotted as points with different shapes at the barycenter and variable groups were connected by lines with each center. The figure was arranged into sub-plots by isotype, tissue (F) and immunization (Figure S26) for easy visualization.