C57BL/6 (WT) and Rag1^–/–/VH81X mice were maintained in accordance with local and Home Office rules and ARRIVE guidelines under Project Licence 80/2529.

For each biological replicate, bone marrow from nine 6- to 8-week-old Rag1^−/−/VH81X mice (47, 48) or from fifteen 12-week-old male wild-type (WT) C57BL/6 mice was incubated with biotinylated antibodies against CD11B (MAC-1; ebioscience), Ly6G (Gr-1; ebioscience), Ly6C (Abd Serotec), TER119 (ebioscience), and CD3E (ebioscience) followed by incubation with streptavidin MACs beads (Miltenyi), to deplete macrophages, granulocytes, erythroid lineage, and T cells. Pre-B cells (surface IgM⁻CD25⁺B220⁺CD19⁺) were then flow sorted on a BD FACSAria in the Babraham Institute Flow Cytometry facility. Antibodies used were CD45R BV421 (B220, RA3-6B2, Biolegend), CD19 PerCP-Cy5.5 (1D3, BD Pharmingen), CD25 APC (PC61.5, eBioscience), and IgM PE (eB121-15F9, eBioscience). Sort purities were all greater than 92%.

BAM files of aligned V reads were loaded using default parameters into Seqmonk version 1.37.1 (Babraham Bioinformatics³), a tool for the visualisation and analysis of next generation sequencing data. Based on examination of the reads, two Vκ gene annotations were corrected (Supplementary Text S2). V reads on the same strand as the gene were then counted within windows extending 750 bp upstream (with respect to the gene orientation) from the 3′ end of each Vκ gene, for each Jκ gene separately. Read counts were normalised to the replicate with the median number of reads (replicate 1), either across all Js, or just for Jκ1-associated reads. For Jκ1-associated V sequences, the median of the normalised replicates for each Vκ gene was used as the recombination frequency in downstream analyses. The mappability (defined as the percentage of all possible 85 bp sequences over a given window that can be mapped uniquely) over 350 bp windows upstream of each Vκ gene 3′ end (within which the vast majority of VκJκ-seq reads are localised) was calculated, and where stated, genes with low mappability (below 70%) were excluded. Only 11 genes fell into this category, 10 of which actively recombine, and only 4 of these had mappability below 60%. Quantitated data and other information relating to each Vκ gene is provided in Table S1 in Supplementary Material.

To facilitate analysis by IMGT HighV-QUEST, a tool for the high-throughput analysis of V(D)J-recombined sequences (50), V and J sequences were first merged, and any gaps filled in, as described previously (14). Igκ rearrangements were then analysed by IMGT HighV-QUEST, with default parameters. Complementarity determining region 3 (CDR3) lengths and productive versus non-productive rearrangement data were obtained from the “Summary” file.

We defined active genes as those that were significantly enriched (padj < 0.01) for VκJκ-seq reads compared to the Vκ region as a whole, using a binomial test. Thus, the probability was defined as the length of the window in which reads were quantified, as a fraction of the total length of the Vκ region, n as the rounded median of normalised read counts for each gene, and N as the rounded median of normalised read counts across the entire Vκ region.

ChIP peaks (including DHS) were called using MACS peak calling algorithm (51). MACS 1.4, which performs better with broad peaks was used for H3K27me3 with pvalue cutoff = 1.00e−02. For all the remaining features, MACS2 was used with pvalue cutoff = 1.00e−05. The number of peaks over the Igκ locus for each dataset is shown in Table S3 in Supplementary Material.

For each Vκ gene, we calculated the distance from the centre of the gene to the closest upstream and closest downstream (defined with respect to the gene orientation) peak summits. If a peak summit was located upstream of the gene centre, and within 1 kb of the gene start site, it was labelled as promoter-associated; conversely any peaks with summits located downstream of the gene centre and within 1 kb of the gene 3′ end were labelled as RSS-associated.

Relative enrichment for ChIP-seq datasets was calculated as the number of reads over a given window, relative to the average number of reads within windows of identical size across the entire Igκ locus.

A phylogenetic tree of C57BL/6 mouse Vκ gene germline sequences, based on NCBI Reference Sequence: NG_005612.1⁴ but with the corrected annotations detailed in Supplementary Text S2, was constructed using the Phylogeny.fr tool, without including alignment curation (52, 53). Multiple sequence alignment was performed with MUSCLE (54), and the maximum likelihood used for tree construction (55, 56). The tree was visualised using the R package ggtree (57), and nodes comprising the majority of each Vκ gene family were collapsed.

Our computational approach comprises an unsupervised and a supervised step. In the unsupervised step, we set out to interrogate the chromatin landscape of the Igκ locus through integration of the histone modification, DNase hypersensitivity (DHS), germline transcription, and TF-binding profiles included in this study. The supervised step is constructed in two layers: first, we train a Random Forest Classifier (RF-C) C(X) to predict whether a given gene is active or not; second, we construct a Random Forest Regression (RF-R) model R(X) to predict the frequency of recombination of a given active gene. In what follows, we describe both steps in more details.

For the unsupervised step, we first trained a RF-C to predict whether a Vκ gene is “active” or not. We then constructed a RF-R model to predict the recombination level of an active gene. We chose RF since it is generally accepted to be superior in tackling high dimensionality (relatively high number of features with low number of samples for the training) and co-linearity between the features (61).

Read counts for DHS-seq, ChIP-seq, and RNA-seq data were generated using Seqmonk, within four distinct windows for each Vκ gene: “promoter,” extending from 500 bp upstream of the start of the gene to its centre; “RSS,” extending from the centre of the gene to 500 bp downstream of its 3′ end; and “upstream” and “downstream” windows, extending from 500 bp to 3 kb up- or downstream of the gene start or end, respectively. In addition to these four windows for each of the genome-wide datasets, giving a total of 76 chromatin features, we also included three genetic features. These were: the RSS RIC score, the orientation (or strand) of the Vκ gene, and the distance from the Vκ gene to Jκ1. All of these features, except for the gene orientation, were projected between 0 and 1. These 79 features were considered as the explanatory variable for both the RF-C and RF-R.

The response variable for RF-C was the binary recombination classes (active and inactive), which were defined as described above. For RF-R, the log2-transformed median of the normalised recombination frequencies of active genes was used as the response variable.

Both the RF-C and RF-R approaches were performed with 10-fold cross-validation: 10% of genes were assigned to the test set each time, with every gene included in a test set exactly once. The number of trees generated for each fold was 1,000. For the initial RF-C including all features, the number of features tried at each step was set to 20; for all other models default parameters were used. The average importance of each feature, and SE across the 10-folds, was recorded. For the classification model, the performance was assessed by calculating the percentage of correct predictions (accuracy) across all ten test sets: this was calculated overall, as well as for the active Vκ genes (giving a measure of the sensitivity with which we could identify an active gene) and inactive genes (which gives a measure of specificity). We also calculated the F1 score as a combined measure of sensitivity and specificity. To assess the performance of the regression model, we used the root mean squared error (RMSE) for the predicted recombination frequencies compared to the observed values across all ten test sets. The RMSE gives a measure of the SD in errors, thus 68% of our predictions are expected to have an error within this range. Since our recombination frequencies are log2-transformed, an RMSE of x corresponds to a 2^x-fold difference between the predicted recombination frequency and the observed recombination frequency.

For model selection, we focussed on the 16 most important features from the initial classification or regression model. We trained RF classification or regression models, with 10-fold cross-validation, for all possible combinations of the respective 16 features. These models were then compared using the performance metrics described above. Our analysis was performed using the R package randomForest (62).

Publicly available genome-wide datasets analysed during this study are available in the GEO repository; details including accession numbers are listed in Table 1. All were downloaded from GEO as raw short-read files (SRA) and realigned to NCBIM37/mm9 using Bowtie (49) or Bowtie 2 (63). The VκJκ-seq datasets generated in this study are available in the GEO repository with accession number GSE101606.⁵ Some of the quantitated data from this study is also provided in Table S1 in Supplementary Material.

To quantify the usage of Vκ genes in an unbiased way from DNA, we adapted our previously reported mouse Igh VDJ-seq assay (14) for the mouse Igκ locus (Figure S1 and Supplementary Text S1 in Supplementary Material). We generated VκJκ-seq data for three biological replicates in wild-type (WT) bone marrow pre-B cells (B220⁺/CD19⁺/CD25⁺/IgM⁻), and one replicate in pre-B cells from a Rag1^−/− mouse with a rearranged Igh transgene (VH81X) (47). These Rag1^−/−/VH81X cells lack the RAG1 recombinase, precluding V(D)J recombination, but progression to the pre-B cell stage is permitted through expression of the VH81X transgene. Thus, they serve as a negative control in which reads mapping to Vκ genes, indicating a VκJκ recombination event, should not be detected, giving a measure of the spurious incorporation of these reads into our libraries. Indeed, while 91.8% of unique, Jκ bait-associated reads for this library mapped upstream of unrecombined Jκ genes, only 9 reads (0.0002%) mapped to Vκ genes. Conversely, in WT pre-B cells over 30% of reads mapped to Vκ genes for all replicates. This equated to a total of 400–530,000 unique VκJκ recombined fragments for each replicate (Figure S2A in Supplementary Material), an order of magnitude greater than previous high-throughput assays of the Vκ gene repertoire (37, 38). In the Rag1^−/−/VH81X library, we noted a slight bias towards Jκ2 (Figure S2B in Supplementary Material), suggesting minor preferential priming of the Jκ2 gene. However, within the VκJκ recombined fragments, over 35% of reads were associated with Jκ1, which usually recombines first (11), indicating that we are capturing a large proportion of the primary Vκ gene repertoire. Analysis of the sequences using IMGT/HighV-QUEST revealed that the ratio of productive:non-productive rearrangements was approximately 37:63 (Figure S2C in Supplementary Material), close to the expected two-thirds of non-productive rearrangements. Consistent with previous reports (8, 9, 38), the vast majority of functional rearrangements had a CDR3 of 9 amino acids (Figure S2D in Supplementary Material).

Using normalised frequencies for each Vκ-Jκ gene combination, we clustered the dataset based on both Vκ and Jκ genes, taking each replicate separately. While Vκ genes with poor RIC scores clearly recombine infrequently, we did not observe any relationship between the repertoire of Vκ genes and their distance from, or orientation with respect to, the Jκ genes. For each Jκ gene, the three replicates were highly correlated (Pearson correlation coefficients >0.992) and clustered more closely to each other than to the repertoires of other Jκ genes, indicating that they are associated with distinct Vκ gene profiles (Figure 1). Notably, the pattern of recombination to Jκ1 segregated from the other Jκ genes. This is consistent with the preferential usage of Jκ1 in the generation of the primary repertoire, while the other Jκ genes are subsequently used for receptor editing (11). Since we aimed to assess the features driving recombination of the germline Igκ locus in the formation of the primary repertoire, we chose to focus on the Jκ1 repertoire for further analyses (Table S1 in Supplementary Material).

The Vκ-Jκ1 repertoire varied widely across the locus, with no clear geographical pattern (Figure 2A). The RIC scores of Vκ genes from the same family were quite homogeneous, while their recombination frequencies could vary by more than 10-fold (Figures 2B,C). Comparing the families, we noted some patterns: for example, Vκ1 genes recombine quite frequently compared to several other families, even when their median RIC score was similar (e.g., Vκ2) or higher (e.g., Vκ4). For all genes with a RIC score > −38.81, which are theoretically considered capable of recombination (12), usage of Vκ genes on the forward and reverse strands was not significantly different (Figure 2D), despite the significantly lower RIC scores of forward compared to reverse strand genes (Figure 2E). This contrasts with observations from an RNA-based assay (37) in which recombination to Jκ1 was biased towards inversional rearrangements. Moreover, while 6 out of the 10 Vκ genes that were most frequently represented in their expressed Jκ1 repertoire (37) are included in the top 20 of our Vκ-Jκ1 repertoire, 4 are not, and the DNA repertoire is not dominated by a small number of genes. Our assay also reveals numerous Vκ genes that are more highly represented at the DNA level, including 14 pseudogenes that were not detected in the RNA repertoire. Conversely, all Vκ genes present in the expressed repertoire were detected in our assay, albeit in some cases with very low read counts. This highlights the significant contribution of transcription and posttranscriptional processes to the expressed repertoire, which would confound the aim of this study to interrogate the pre-recombination chromatin state.

To facilitate further investigation of the Vκ-Jκ1 repertoire, we performed a binomial test to distinguish Vκ genes that are significantly recombining (padj < 0.01, Figure 3A; Table S1 in Supplementary Material). Out of 162 genes, 105 (64.8%) passed the binomial test and were labelled “active” to denote “actively recombining”; these genes were detected with a minimum of 59 reads, and included 15 pseudogenes, and are hereafter referred to as active genes. The remaining 57 genes had insufficient evidence of activity, with the median read count for each below 39, and were labelled “inactive.” These inactive genes included eight Vκ genes that are considered to be functional, suggesting that they contribute little to the primary repertoire. The usage of active genes was weakly correlated with the RIC score (R = 0.42; Figure 3B); however, genes with a similar RIC score could recombine at markedly different frequencies. Moreover, some inactive genes have RIC scores that are comparable to those of active genes (Figure 3C). Importantly, a linear regression model revealed that only 17.7% of the variability in Vκ gene usage could be explained by RIC score alone (Figure 3B), highlighting the need to explore whether other mechanisms, such as chromatin features, contribute to shaping the repertoire.

We next sought to understand how genetic and chromatin features regulate Vκ gene usage. Thus, we trained a RF-C to assess the power of each feature to correctly predict whether a gene is active or inactive; we have previously used this approach for analysis of Igh recombination (14). The RF-C takes features relating to each sample (here, each Vκ gene), and generates a large number of decision trees that vote on the response (here, Vκ gene recombination activity) for a training set of samples. At each step, one feature out of a random subset can be chosen, such that each tree will be unique. Feature importance is gauged by comparing the prediction accuracy for trees that include or exclude a given feature. The overall accuracy is then assessed by predicting the activity of an independent test set of genes. We chose the RF approach because it performs well both with a large number of features relative to the number of samples, and with highly correlated features (61).

We considered the signal intensity over the promoter separately from the signal intensity over the RSS, taking windows from the centre of each gene extending 500 bp upstream (promoter) or downstream (RSS) of the gene. We also calculated the signal intensities in windows extending a further 2.5 kb upstream of the promoter and downstream of the RSS (Figure 7A). In addition, we systematically gauged the contribution of the orientation of each Vκ gene and their genomic distance from Jκ1 to facilitating recombination. Thus, the input for each gene included three genetic features: the RIC score, strand, and distance from Jκ1; and four separate features for each chromatin dataset (hereafter referred to as promoter, RSS, upstream and downstream).

A 10-fold cross-validation approach including all features revealed a mean prediction accuracy of 93.3% (SD 5.3%), with an F1 score of 0.949 (SD 0.0405). The prediction accuracy for active genes (97.2%, SD 6.2%) was better than that for inactive genes (85.7%, SD 14.5%), indicating a high sensitivity but slightly lower specificity in detecting active genes. The RIC score was by far the most important feature in distinguishing active from inactive genes (Figure 7B). The second most important feature was PU.1, which contributes to activation and recombination of the Igκ locus (21). Our model specifically suggests that while PU.1 binding at the promoter is of no consequence, its binding at the RSS is an important driver of recombination.

H3K4me2 enrichment within both the promoter and the RSS windows was also important for prediction accuracy, in addition to MED1 binding at the promoter and a number of E state-associated chromatin features at the RSS (Figure 7B). Notably, PAX5, CTCF, and RAD21, which are key drivers of Igh recombination (14), were not identified as important in promoting Vκ gene activity. Several up- or downstream features, such as DHS upstream and IRF4 binding downstream, were, however, ranked quite highly (Figure 7B), suggesting that in addition to the chromatin state of the Vκ genes themselves, the surrounding chromatin also influences the capacity of each gene to recombine.

Next, we performed a model selection analysis, considering all possible combinations of the 16 most important features. RF-C models that included the RIC score were all more accurate than those that did not (Figure 7C). This clear distinction was driven by the much lower prediction accuracy for inactive genes in RF-C models in which the RIC score was excluded (Figure 7D). There was also a striking contribution of PU.1 binding at the RSS to prediction accuracy, particularly in models that excluded the RIC score (Figure 7E); indeed, the bimodal distribution in prediction accuracy observed for models excluding the RIC score appears to be primarily dependent on this feature. This was apparent for both active and inactive genes (Figures S7A,B in Supplementary Material). We also observed a slight shift in prediction accuracy for models that included several other important features, such as IKAROS enrichment at the RSS (data not shown).

While the RF-C approach identifies the features that are most highly predictive of active Vκ gene recombination, and thus likely facilitate recombination, it does not directly show whether their levels contribute to the frequency with which an active gene is used. To address this, we trained a RF-R model for active Vκ genes, using the same set of features used for the classification model. This extends beyond the RF-C approach, giving a numerical prediction of the recombination frequency, with the log2-transformed recombination frequencies of active genes as the response variable. We used the RMSE to measure model performance. An initial RF-R using all features, with 10-fold cross-validation, achieved an RMSE of 1.57, indicating that 68% of the models’ predictions fall within 2^1.57-fold (since recombination values are log2-transformed) of the observed recombination. This revealed the RIC score to be the most important feature for predicting the frequency of active Vκ gene recombination, as in the RF-C model (Figure 8A). However, there was far less separation between the importance of the RIC score and the other features. PU.1 binding at the RSS was much lower in importance compared to the RF-C model, while several features that were unimportant for classification contribute much more to the prediction of recombination frequency: these include PAX5 and CTCF binding upstream, and H3K4me3 and E2A binding at the RSS.

We next performed a model selection analysis, considering the top 16 features, to identify a minimum subset that best predict Vκ gene usage. While several combinations of 6–9 features have similarly low RMSE (Figure 8B), the minimum RMSE of 1.36 (Figure 8C) is achieved with a combination comprising: RIC score, H3K4me2, and H3K4me3 within the RSS window, IKAROS binding at both the promoter and the RSS, and PAX5 and PU.1 binding upstream. These seven features, in addition to H3K4me2 within the promoter window, were also the most highly represented in all models that had a low RMSE (<1.39), with the RIC score and H3K4me3 at the RSS being present in all models (Figure 8D). While this does not necessarily mean that these features are the most important individually (compare to Figure 8A), it suggests that together they are able to explain the largest proportion of the variability in our data. In models that excluded the RIC score, the best combination had an RMSE of 1.44 (Figure 8E), and comprised H3K4me3 at the RSS, IKAROS binding at both the promoter and the RSS, and PAX5 and PU.1 binding upstream, in addition to MED1 binding at the promoter. These six features, in addition to CTCF binding upstream, were the most frequently represented in models with low RMSE (<1.48) that excluded the RIC score, with MED1 binding to the promoter in all models (Figure 8F); this is particularly noteworthy since this feature was rarely present in the best combinations that included the RIC score. The inclusion of H3K4me3 at the RSS in a model was accompanied by a shift towards lower RMSE; this shift was more pronounced for the left-hand peak, corresponding to models that also include the RIC score (Figure 8G). Conversely, a shift towards lower RMSE for models that included MED1 binding at the promoter was only evident for the right-hand peak, representing models excluding the RIC score (Figure 8H). We also noted a shift towards a lower RMSE for models that included several other important features, including IKAROS binding at the promoter or RSS, or PAX5 or PU.1 binding upstream (Figure S7C in Supplementary Material).

These analyses revealed that we are able to predict the recombination frequency of Vκ genes from a combination of 6 or 7 features with a mean error rate of less than threefold, even when the RIC score is excluded (2.71-fold compared to 2.57-fold when including the RIC score). This model performance is highly significant when noting the variability of greater than 150-fold in recombination frequency across all active Vκ genes.

In order to further dissect the influence of the chromatin features shown to be important in our RF models, we split all active genes with a good RIC score (>−14) into three groups based on their recombination frequency, and examined the enrichment of features of interest at the location(s) in which they were important (Figure 9). Genes that recombine more frequently tended to display higher levels of active histone modifications and TF binding at both the RSS (Figure 9A) and the promoter (Figure 9B). These trends were particularly striking for H3K4me2 (RSS and promoter windows) and IKAROS (RSS), in addition to IRF4 and E2A binding at the RSS, all of which displayed a significant positive relationship between recombination and enrichment. For PU.1 binding at the RSS, which was highly important in distinguishing active from inactive genes, the difference between genes that recombine at low and high frequency was more subtle (p > 0.05), consistent with the lower importance of this feature for predicting the frequency of recombination. We also noted a subtle, but non-significant, trend towards higher enrichment of MED1 and IKAROS at the promoter, and H3K4me3 at the RSS, of more highly recombining genes: the importance of these features in the RF-R thus suggests that a more complex relationship exists between these features and recombination frequency. Conversely, we noted a slight, negative association between the recombination frequency and the binding of some TFs upstream of the gene, including PAX5, PU.1, and CTCF, with a significantly greater enrichment of PAX5 upstream of genes that recombine at a low level compared to those that recombine at a medium level (Figure 9C).