Accurate identification and quantification of commensal microbiota bound by host immunoglobulins

Identifying which taxa are targeted by immunoglobulins can uncover important host-microbe interactions. Immunoglobulin binding of commensal taxa can be assayed by sorting bound bacteria from samples and using amplicon sequencing to determine their taxonomy, a technique most widely applied to study Immunoglobulin A (IgA-Seq). Previous experiments have scored taxon binding in IgA-Seq datasets by comparing abundances in the IgA bound and unbound sorted fractions. However, as these are relative abundances, such scores are influenced by the levels of the other taxa present and represent an abstract combination of these effects. Diversity in the practical approaches of prior studies also warrants benchmarking of the individual stages involved. Here, we provide a detailed description of the design strategy for an optimised IgA-Seq protocol. Combined with a novel scoring method for IgA-Seq datasets that accounts for the aforementioned effects, this platform enables accurate identification and quantification of commensal gut microbiota targeted by host immunoglobulins. Using germ-free and Rag1−/− mice as negative controls, and a strain-specific IgA antibody as a positive control, we determine optimal reagents and fluorescence-activated cell sorting (FACS) parameters for IgA-Seq. Using simulated IgA-Seq data, we show that existing IgA-Seq scoring methods are influenced by pre-sort relative abundances. This has consequences for the interpretation of case-control studies where there are inherent differences in microbiota composition between groups. We show that these effects can be addressed using a novel scoring approach based on posterior probabilities. Finally, we demonstrate the utility of both the IgA-Seq protocol and probability-based scores by examining both novel and published data from in vivo disease models. We provide a detailed IgA-Seq protocol to accurately isolate IgA-bound taxa from intestinal samples. Using simulated and experimental data, we demonstrate novel probability-based scores that adjust for the compositional nature of relative abundance data to accurately quantify taxon-level IgA binding. All scoring approaches are made available in the IgAScores R package. These methods should improve the generation and interpretation of IgA-Seq datasets and could be applied to study other immunoglobulins and sample types. BHGKVvENxrBpamSANr5RyR Video abstract Video abstract


Background
The immunoglobulin isotype A (IgA) is secreted by plasma cells at mucosal barrier sites such as the gastrointestinal (GI) tract, which is also host to the densest community of commensal bacteria in the human body [1,2]. Intestinal homeostasis requires tolerance to the microbiota but a robust defence against pathogens, and growing evidence suggests that IgA binding is both taxon-dependent and context-specific. IgA in the mucus layer forms a protective barrier between the epithelium and the luminal microbiota [3][4][5]. Some commensal bacteria have adapted to this specialized niche to promote homeostasis. Bacteroides fragilis produces a surface polysaccharide that promotes IgA binding and facilitates mucus layer colonization [6] and the mucin-degrader Akkermansia muciniphila induces host IgA and IgG responses [7]. Other studies have suggested that lowaffinity IgA contributes to microbiota maintenance, whilst high-affinity IgA acts on pathogens to promote clearance and inhibit virulence [8][9][10][11][12][13][14][15]. Pathobionts are members of the microbiota that are typically harmless under homeostatic conditions but can drive disease given certain environmental stimuli [16]. IgA binding of taxa might also reflect such pathogenic potential and IgA bound gut bacteria have been implicated in inflammatory disease [17][18][19][20]. Accurate identification of commensal microbiota bound by IgA can therefore offer important insights into the mechanisms underlying host microbial mutualism and its dysregulation in disease.
Microbiota-wide profiling of IgA bound taxa can be achieved using IgA-Seq. Anti-IgA antibodies are used to stain a complex gut microbiota, then bound and unbound bacteria are separated by fluorescence-activated cell sorting (FACS) or magnetic-activated cell sorting (MACS) and identified by 16S rRNA gene sequencing [17,18]. We refer to the sorted IgA bound bacteria as the IgA + fraction, the sorted unbound bacteria as the IgA − fraction and the sample before sorting as the presort sample. IgA-Seq has been used to uncover antibody targeting of taxa by both T cell-dependent and T cellindependent pathways [7,10,12,21,22], as well as differential binding of taxa during disease. IgA-Seq was first described by Palm et al. to identify taxa differentially bound by IgA in inflammatory bowel disease (IBD) patients [17] and was adapted by Kau et al. to identify increased binding of Enterobacteriaceae by IgA in malnourished infants [18]. Variations of these protocols have since been applied to study IgA binding in other disease states, employing a range of IgA-Seq protocols utilising different reagents, sort modalities and configurations [19,20]. Some of these parameters have been validated individually in prior studies but there has yet to be a comprehensive description of how to design and benchmark an optimal IgA-Seq protocol.
More importantly, there has yet to be a formal consideration of the analytical approaches used to score taxon binding in IgA-Seq datasets. IgA-Seq does not quantify the affinity of IgA antibodies directly, but rather provides a measure of overall IgA binding. The likelihood of bacteria within a given taxon being bound by IgA will be influenced by a combination of the affinities of the IgA pool for the taxon, the spatial distribution of those bacteria and the expression of relevant epitopes by the taxon. Additionally, the IgA + abundance does not solely represent this likelihood of taxon binding. It is also a function of the number of bacteria present: a taxon highly targeted by host IgA could make up a low proportion of the IgA + fraction if only a few bacteria are present and vice versa. The principal aim of IgA-Seq is not to solely identify which taxa are most bound by IgA but rather to identify those taxa with a higher likelihood of binding relative to other bacteria. As a result, indices that account for the initial quantities of bacterial taxa must be applied to quantify IgA-Seq data.
In IgA-Seq studies to date, IgA binding has largely been scored using indices as described by either Palm et al. or Kau et al. (referred to here as the Palm and Kau indices for clarity) [17,18]. These scores compare the abundance of a taxon in the IgA + fraction to its abundance in the IgA − fraction. This accounts for the overall pre-sort quantities of taxa on the assumption that, as both fractions are drawn from the same starting pool of bacteria, the IgA + and IgA − abundances are representative of the proportion of a taxon's bacteria bound by IgA. However, these scores are influenced by an additional confounder-that these are relative abundances. As a result, an increase in one taxon will inherently reduce the percentage contribution of others. Thus, the IgA + and IgA − abundances of a given taxon are also influenced by the starting abundances and IgA binding of all other taxa. Scores like the Kau and Palm indices correlate with taxon-level IgA binding but the necessity to use relative abundances in IgA-Seq experiments abstracts their interpretation. We will show that, without careful consideration of these effects, it is possible to misinterpret results observed using these metrics.
Here, we aim to address these methodological and analytical limitations to improve the accessibility, quality and quantitative interpretation of IgA-Seq experiments. Using germ-free and immune-deficient mice as negative controls and a strain-specific IgA antibody for positive controls, we describe the methods that can be used to validate the quality of each stage of an IgA-Seq experiment and derive an optimised IgA-Seq protocol. In lieu of a gold standard for benchmarking IgA binding scores, we develop a platform for simulating an IgA-Seq experiment where the underlying relative IgA binding levels between taxa are known. Using this, we show that existing indices used to score taxon-level IgA binding are influenced by relative abundances, which needs to be taken into account when drawing conclusions. We present novel probability-based scores that circumvent this by providing direct quantification of the proportion of bacteria belonging to a given taxon that are bound by IgA. Finally, we demonstrate the utility of the protocol and novel probability scores by applying them to study IgA binding in mouse models of colitis and malnourishment. These methods facilitate generation of highquality IgA-Seq data and direct quantification of relative immunoglobulin binding between taxa.

Quality control standards for bacterial flow cytometry and IgA staining
We first determined the optimal gating strategy for isolating bacteria from faecal and colonic samples by FACS. The BD LSRFortessa X20 and FACSAria III flow cytometers were selected based on an identical configuration with 405, 488, 633 and an additional 561 nm laser, which eliminated PE-FITC spill over and increased PE sensitivity. Whereas eukaryotic cells have diameters ranging from 10 to 100 μm, bacteria generally range from 0.3-5 μm, though some can vary by two orders of magnitude and have variable width to height ratios. In stock configurations, both these cytometers use a sensitive photomultiplier tube (PMT) to detect side scatter (SSC), and therefore a minimum SSC threshold was used to fully capture a wide range of particle sizes.
Bacteria must then be distinguished from small debris particles in the buffers and the sample itself, particularly in faecal samples which contain substantial autofluorescent debris. For this, nucleic acid dyes have been recommended in flow cytometry applications for diverse hosts and environments [23][24][25][26][27][28], and help ensure collection of adequate DNA for IgA-Seq studies [12,[18][19][20][29][30][31]. It must be noted that nucleic acid staining may be variable amongst diverse bacteria due to differences in cell permeability and transport mechanisms, the size and number of copies of genomic DNA and/or asymmetric division. Studies comparing various DNA dyes have found that the cell-permeable SYBR Green I effectively stained phylogenetically distant bacterial strains, was the most sensitive and bound specifically to DNA; in contrast to DAPI which bound non-bacterial particles [25,26]. This reagent had also been used successfully in previous studies of immunoglobulin coating [30,31], and was thus selected for our IgA-Seq protocol. To accommodate bacteria that stained weakly with SYBR, we utilised a broad gating strategy to capture as much of this heterogeneity as possible (Fig. S1A). Our SYBR + gating strategy identified bacteria in samples from specific pathogen-free (SPF) mice but did not stain filter-sterilized FACS buffer or stool from germ-free (GF) mice (Fig. S1A). This approach defined SYBR + gating parameters that were specific for microbes whilst minimising contaminating debris (Fig. S1A).
Next, we sought to validate the specificity of the PEconjugated rat anti-mouse IgA (mA-6E1) antibody that would be used in downstream IgA-Seq experiments. We used IgA-deficient mice as negative controls and a bacterial strain-specific IgA antibody as a positive control. SPF Rag1 −/− mice have a complex microbiota but lack T cells and mature B cells responsible for the production of IgA. Essentially, no bacteria were stained by the PE anti-mouse IgA antibody in samples from Rag1 −/− mice compared to roughly 4% of the microbiota in wild-type SPF mice (Fig. S1B). To further confirm the specificity of the IgA antibody, we replicated the approach of Kau et al. [18], using an in vitro system with an artificial 2member bacterial community in which only one strain should be stained by PE anti-mouse IgA. Bacteroides thetaiotaomicron (B. theta) strain VPI-5482 was coated with mouse IgA antibody clone 225.4, which is specific to the capsular polysaccharide of strain VPI-5482 and no other B. theta isolates [32]. Escherichia coli was coated with a mouse IgG antibody as a control and mixed with the labelled B. theta prior to staining with the PE antimouse IgA antibody. Flow cytometry revealed that a proportion equivalent to the expected B. theta population was stained by the PE anti-mouse IgA antibody (Fig. S1C), confirming the specificity of the antibody to be used in downstream IgA-Seq experiments. The gating strategy used for the IgA-Seq protocol (Fig. 1a) was validated using these SYBR and IgA staining controls.

A benchmarked protocol for IgA-Seq
The goal of IgA-Seq is to isolate bacterial populations coated by endogenous host IgA and to identify them by 16S rRNA gene sequencing. To this end, it is critical to identify FACS parameters that maximise the integrity and purity of sorted samples. Using B. theta as a model commensal anaerobe, a culture-grown suspension of B. theta VPI-5482 was once again labelled with strainspecific IgA and PE anti-mouse IgA, and 100 cells were sorted as an array onto an agar plate to determine viable colony-forming units (CFU). The default cytometer settings for 100 and 85 μm nozzles were tested, in addition to a custom 30 psi configuration for the 85 μm. Intact bacterial recovery decreased with smaller nozzles sizes and higher pressures, and viability was enhanced by sorting at the lower 30 psi pressure on the 85 μm nozzle (Fig. 1b). Only 20% of bacteria could be recovered with the 70 μm nozzle (Fig. 1b), which was discontinued from further testing.
To test the sort fidelity of the fractions produced by each of the configurations, we sorted IgA + and IgA − bacterial populations from SPF mouse stool (Fig. 1a). These fractions were then analysed using the same parameters to determine the percentage of bacteria that were retained within defined sort gates. All three FACS configurations enabled > 90% sort purity for both IgA + and IgA − fractions (Fig. 1c). This is in contrast to magnetic-activated cell separation (MACS), an alternative sorting modality which performed poorly with purities of 37.1-42.9% for the IgA + fraction (Fig. 1d, e). The FACS configuration utilising the 100 μm nozzle produced the highest sort purities in both IgA + and IgA − fractions, Fig. 1 Benchmarking a protocol to isolate IgA + and IgA − faecal bacteria by fluorescence activated cell sorting. a IgA-Seq gating strategy. Faecal bacteria from SPF C57BL/6 mice stained with anti-mouse IgA (PE). Gating on a population based on FSC and SSC characteristics, followed by discrimination of SYBR + (FITC) bacteria from autofluorescent (PerCP Cy5.5) debris. IgA + and IgA − populations were sorted from SYBR + cells. b Viability of B. theta VPI-5482 (Bt) collected by various nozzle sizes (μm) and sheath pressures (psi) on the FACSAria III. Culture-grown Bt was stained with mouse anti-Bt IgA (clone 255.4) and rat anti-mouse IgA (PE, clone mA-6E1). One hundred cells were plated in single cell precision mode onto BHIS agar plates for enumeration of colony forming units (CFU). c IgA + and IgA − fractions were collected from SPF C57BL/6J mice (n = 3). Left: percentage of IgA + faecal bacteria detected in each collected fraction for a range of nozzle and pressure configurations. Right: sort purities for these configurations defined by the percentage of events that remained in the designated sort gates. Lines represent means with SEM. dA portion of the IgA-stained faecal suspensions from Fig. 2a were also separated into IgA + and IgA − fractions by magnetic activated cell separation (MACS) as a comparison. PE IgA-stained samples were incubated with anti-PE microbeads and fractions collected by LS column separation. Representative flow cytometry plots of the starting material in SPF stool and the IgA + and IgA − fraction isolated by FACS and MACS. e MACS purity of collected fractions as determined by flow cytometry, technical triplicates from the same stool suspension followed by the 85 μm nozzle at 30 psi (Fig. 1c). Practically, the low pressure required for the 100 μm nozzle reduces the flow rate and throughput by~70%. Given that a minimum of 400,000 cells are required per fraction for downstream sequencing applications, the 85 μm with 30 psi configuration is vastly more time effective whilst maintaining high viability and sort purities. We also found this configuration allowed for sensitive sort control that, if desired, allows for finer grain collection of IgA + subpopulations for instance IgA bright and IgA dim fractions (Fig. S2). This configuration was selected for downstream IgA-Seq.
As a final test of the precision of our sorting approach, we devised an experiment to specifically isolate an IgA + strain from a complex stool microbiota. We once again took advantage of our ability to coat B. theta VPI-5482 with mouse IgA in a strain-specific manner, which we stained separately before spiking into an unstained faecal suspension prepared from SPF wild-type C57BL/6 mice (Fig. 2a). Though some members of the faecal microbiota are coated by endogenous host IgA, because the faecal suspension was not stained, only the PE antimouse IgA stained B. theta VPI-5482 culture should be isolated with this approach. B. theta VPI-5482 was not found in the SPF stool without spike-in and comprised roughly 12.5% of the whole community in the spiked stool mixture, and high sort purities were achieved for both the IgA + (97.4-99.4%) and IgA − fractions (97.3-100%) across three replicate sorts (Fig. 2b). 16S rRNA gene sequencing was performed on the complex stool community prior to the addition of B. theta VPI-5482, and in addition to the IgA + fractions after three sorts from the same spiked sample. Genomic DNA from B. theta was included as a positive control. As a contamination check, sheath fluid and an extraction blank (no sample added) were also included as controls, but these contained undetectable quantities of DNA and did not produce a viable library. Whilst endogenous Bacteroides species are present in the faecal microbiota (Fig. 2c), amplicon sequence variant (ASV) level analysis confirmed that strain VPI-4582 is absent from C57BL/6 mice (Fig. 2d). The IgA + fractions were comprised almost entirely of B. theta VPI-5482 (ASV399, 98.4-99.6%; Fig. 2c, d) and lacked any additional unstained taxa observed in the complex community. This demonstrates that our optimised protocol can sort IgA bound bacteria from a complex gut microbiome with high accuracy.
Pre-sort differences in taxonomic composition and the relative nature of abundances influences existing IgA binding scores Once a sample has been sorted, amplicon-based sequencing of extracted DNA from the sort fractions is used to determine the relative abundances of the taxa within them. As discussed previously (see "Background"), the abundance of a taxon in the IgA + fraction is a function of both its likelihood of IgA binding and its overall abundance in the pre-sort community. As we are principally interested in the former, we must use scoring metrics to quantify the relative IgA binding of taxa whilst taking their pre-sort abundances into account. To date, this has principally been achieved using indices that compare taxon abundances in the IgA + and IgA − fractions and two such indices have been described-the Palm and Kau indices (see "Methods"). Using a toy example where two samples have identical IgA binding affinities for two species but differ in their starting abundances, we can see that these scores are still influenced by pre-sort abundances and produce different affinity scores for each sample (Fig. 3a). This is because these scores do not account for the fact that abundances are relative measures and thus the starting abundance and IgA binding of all other taxa will also influence the observed abundance of a given taxon in the sorted fractions.
If we instead consider the IgA + relative abundance as the probability of a taxon being within the IgA + fraction, it is possible to apply Bayes' theorem to derive a posterior probability of IgA binding for each taxon, i.e. the probability that a given bacterium will be bound to IgA given that it belongs to a given taxon (see "Methods" for a more detailed discussion). In our example, this posterior probability is not influenced by the effect of pre-sort taxonomic composition on relative abundances and provides a direct quantification of the likelihood of IgA binding for each taxon (Fig. 3a). It is also possible to apply the same approach to the IgA − fraction and estimate the probability a bacterium will be in the IgA − gate given that is belongs to a given taxon. Here, we describe two scoring methods based on these posterior probabilities: (1) the IgA + probability-direct use of the IgA + fraction probability; (2) the probability ratio-a composite score derived from the ratio of the IgA + fraction probability to the IgA − fraction probability (see "Methods").
Probability-based scores accurately recapitulate relative IgA binding in simulated data As the true level of IgA binding to each bacterial taxon cannot be measured directly in vivo, we created a platform to simulate an IgA-Seq experiment in silico to benchmark scoring methods. This allowed IgA-Seq data to be generated with pre-defined relative levels of IgA binding between species. In each simulation, we maintain the relative level of IgA binding for each species across all samples and only change the pre-sort abundances of species. This allows us to determine how robust scoring methods are to these effects. A perfect scoring method should produce identical IgA binding scores for each species across all samples, with a relative ranking of IgA binding between species that matches the pre-defined levels.
Ten species were simulated, each having a unique normal distribution of IgA binding. IgA binding was represented by an arbitrary value that just represents the relative levels of binding between taxa (Fig. 3b). Thirty samples were generated each containing 100,000 bacteria. The bacteria were assigned to each of the ten species randomly to produce different log distributed starting abundances for each sample (Fig. 3c). Each bacterium in each sample was then assigned an IgA binding value by sampling from the distribution of its species. Bacteria whose binding values were above four were considered IgA + and those below two in the IgA − fraction. These thresholds were chosen based on the overall distribution of IgA values to mimic the extremes selected in a typical IgA-Seq experiment (Fig. 3b). The number of bacteria in each sample's IgA + and IgA − fraction was used to determine the size of the fractions and the species counts were converted to taxonomic relative abundances (Fig. 3c). These values were in turn used to derive the various IgA binding scores.
Samples had different IgA + and IgA − fraction sizes, defined as the percentage of total bacteria within the fraction, due to their different initial taxonomic compositions and the differing IgA binding of the ten species (Fig. 3d). Given that each taxon's IgA binding was constant across samples, we would expect there to be no variation across samples if an IgA binding score properly adjusts for the differences in sample composition. Comparing the consistency of scoring methods Fig. 3 The influence of pre-sort abundances on existing IgA-Seq scores can be overcome using probability-based approaches. a A toy example where two samples have two species with the same IgA binding profiles but different starting abundances. b The distribution of arbitrary IgA binding values used for the ten species in the IgA-Seq simulations, dashed lines indicate the IgA + and IgA − fraction thresholds. Species coloured as in (c). c The species abundances in the pre-sorted, IgA + and IgA − samples for the 30 samples simulated. d The size of the IgA + and IgA − fractions in the simulated samples as a percentage of their total bacteria. e Comparison of variance in IgA binding estimates across samples using each scoring approach. Coefficient of variance was calculated for each species individually. **p < 0.01, ****p < 0.0001 in Mann-Whitney tests after FDR correction. Boxes represent the median and interquartile range (IQR) and the whiskers the largest and smallest values within one and half IQR of the upper and lower IQR limits. f Pearson and Spearman correlations between the true mean IgA binding scores used to simulate data and the scores estimated by each of the different indices across samples, the IgA + probability had a significantly lower coefficient of variation than the Palm index and the probability ratio had significantly lower variation than all other scores. The probability-based methods were more consistent than the Palm and Kau indices, and thus less influenced by the differences in pre-sort taxonomic composition between samples (Fig. 3e).
When considering the relative rankings of the species, all four scoring methods were significantly correlated with the true IgA binding affinities and ranked species in the correct order of IgA binding (Fig. 3f). However, the Palm and Kau indices had a lower correlation than either of the probability-based methods due to the larger variation between samples. Comparing the results of Pearson and Spearman correlations, it is apparent that the Kau and probability ratio scores produce a more linear distribution than the other two methods. This is due to their consideration of both the IgA + and IgA − fractions, which is not in the case in the IgA + probability, and taking the log to centre the score, which is not the case in the Palm index or IgA + probability. Overall in this initial simulation, probability-based methods proved robust to the influences of pre-sort taxon levels on postsort relative abundances, which introduced variability to the Palm and Kau indices. The probability ratio almost perfectly recapitulated the relative relationships of the true IgA binding values.

Probability-based methods overcome biases effecting the interpretation of between group comparisons
The Kau and Palm indices showed variation in scores due to pre-sort differences in taxonomic composition. This could be particularly problematic when comparing IgA binding between groups that have inherent and consistent differences in microbiome composition. This is likely, as a key use-case for IgA-Seq is to determine if previously established microbiome differences between groups, such as healthy controls and disease patients, are being driven by differential targeting of the IgA pool.
To determine the effects of pre-sort compositional differences on inter-group comparisons, we repeated the simulation but divided samples into a case and control group. The simulation was run with 30 samples per group each containing ten species with the same underlying IgA affinities as before (Fig. 3b). Random initial abundances were then generated as previously but additional Species 10 (with a moderate IgA affinity value) was introduced to the 30 case samples-simulating outgrowth of one specific taxon in a disease. In this simulation, as before, an optimal scoring approach will not detect any differences between the cases and controls whose species share the same IgA binding and only differ in their pre-sort composition.
As expected, Species 10 constituted a larger percentage of the relative abundances in both the IgA + and IgA − fractions of the cases compared to the controls (Fig. 4a). This led to a significantly higher IgA − fraction size in the control group when compared to cases but no difference in the IgA + fraction sizes (Fig. 4b). Comparing the overall distribution of scores, there was greater variability across samples in the Kau index than the probability ratio, and Species 10 IgA binding estimates were higher in cases using the Kau index (Fig. 4c). When testing differences between the simulated case and control groups, the Kau index inferred significant differences in the IgA binding of both Species 10 and Species 6 where there should be none (Fig. 4d). This is an artefact of measuring relative abundances. In the simulated disease cases, there was an increased abundance of Species 10 in the IgA + fraction. This reduces the relative contribution of Species 6 (the highest affinity taxa that would otherwise constitute the majority of the IgA + fraction) and alters the ratio when calculating the Kau score. No significant differences were observed using the Palm index but this only resolved the highly bound Species 6, which displayed increased variance in the controls (Fig. S3). The probability ratio resolved both high and low IgA binding of taxa and was more resistant to the influence of the pre-sort abundances, maintaining similar estimates between cases and controls (Fig. 4d).

IgA-Seq in a mouse model of colitis produces different outcomes across scoring methods
To demonstrate the utility of our benchmarked IgA-Seq protocol, we applied it to study commensal IgA binding in a mouse model of colitis. Helicobacter hepaticus (Hh) can inhabit the murine intestinal tract without harm but can emerge as a pathogen when host immune responses are dysregulated as a result of either genetic susceptibility or targeted immunomodulation [33]. Abrogation of anti-inflammatory IL-10 signaling by administration of an anti-IL10 receptor antibody (aIL10R) results in Hhdriven colitis [34]. We recapitulated this model in a gnotobiotic C57BL/6 mouse colony stably colonized with a 12-member microbiota (MM12) that recapitulates the taxonomic diversity of a more complex community [35]. This provides a defined minimal microbiota to specifically interrogate IgA binding.
Control MM12 mice treated solely with weekly aIL10R injections did not develop any pathology (Fig. 5a). In contrast, MM12 mice that received aIL10R in conjunction with H. hepaticus (Hh+aIL10R) developed marked inflammation and colitis after 14 days (Fig. 5a). Lipocalin-2, released by activated neutrophils and a biomarker for inflammation [36], was also elevated in the colon contents of Hh+aIL10Rmice but not control aIL10R only mice, whose lipocalin-2 levels were similar to steady-state mice lacking aIL10R injection (Fig. 5b). There was a significant increase in IgA-coating of colonic bacteria from colitic mice (Fig. 5c, d). Colonic contents from both control and colitis groups were subjected to IgA-Seq to identify taxa targeted during inflammation.
Taxa were profiled in the pre-sort samples and the IgA + and IgA − fractions using 16S rRNA gene Fig. 4 The probability ratio overcomes pre-sort abundance biases when carrying out between-group comparisons. a Simulations were carried out as in Fig. 3 with an additional 30 samples that were initiated with an exaggerated abundance of Species 10 (cases). This plot shows the species abundances in these case and control samples' pre-sort, IgA + and IgA − fractions using the thresholds as in Fig. 3b. b Differences in the IgA + and IgA − fraction sizes as a percentage of total bacteria between cases and controls. P value shown from Mann-Whitney tests. c Heatmaps showing the scores estimated by the Kau and probability ratios across all samples. d Comparison of IgA binding scores for the ten species in cases and controls when using either the Kau score or the probability ratio. Significant differences are highlighted with * (FDR adjusted Mann-Whitney p < 0.05) and species are ordered by true IgA binding value. In all boxplots the boxes represent the median and interquartile range and the whiskers the largest and smallest values within one and half IQR of the upper and lower IQR limits. In (d), points outside this range are shown Fig. 5 IgA-Seq in a mouse model of colitis with a defined 12-member microbiota. MM12 mice were administered weekly injections of anti-IL10R (day 0, 7). Controls received anti-1L10R alone (aIL10R, n = 4), and the colitis group were additionally infected with H. hepaticus at the start of the experiment (Hh+aIL10R, n = 4). a Experimental design and representative colon histology with haematoxylin and eosin staining. b Mouse lipocalin-2 ELISA as a measure of colonic inflammation. Lines at means and SEM. **p < 0.01 unpaired Mann-Whitney test, ns non-significant. SS represents steady state mice with no aIL10R exposure. c Representative flow cytometry plots of the percentage of colonic MM12 bacteria coated by IgA. d Percentage of bacteria coated by host IgA and subjected to IgA-Seq. P values represent Mann-Whitney tests. e Relative abundance plots showing the MM12 species in each sample. Each row represents an individual animal. f Non-metric multi-dimensional scaling of pre-sort, IgA + and IgA − fractions from Bray-Curtis distances. g Left: means and 95% confidence intervals of IgA binding scores significantly different between groups. Significance determined by complete permutation of group labels where, exact permuted p < 0.1). Right: Mean difference of scores between groups (colitic -control). Mean is normalised by its standard deviation for comparability between scoring methods (using the strictly standardised mean difference SSMD). h Left: boxplots of A. muciniphila relative abundance. Centre: boxplots of IgA + fraction size. Right: boxplots of adjustment for the difference in fraction size (abundance in the IgA + fraction size multiplied by A. muciniphila abundance). In the boxplots in (d) and (h), boxes represent the median and IQR and the whiskers the largest and smallest values within one and half IQR of the upper and lower IQR limits sequencing. Analyses were carried out at the level of species summarised ASVs. Using the V4 16S rRNA gene region, we were able to distinguish almost all members of the MM12 and Hh either by species level assignment or sole membership of a genus (Table S1). The exception was two closely related Clostridium species that were both assigned to Clostridium XlVa. A blank sequencing control was used to identify potential contaminants introduced at the DNA extraction and library preparation stage. Given that we used mice with a defined microbiota, we were also able to determine the presence of contaminants introduced elsewhere. After removing taxa observed in the blank control, several unexpected (non-MM12 or Hh) species remained (Fig.  S4A). The majority of these were at a lower abundance than the expected species across all samples. These were removed from later analyses using an abundance threshold chosen to remove rare taxa whilst retaining the majority of total observations (Fig. S4B). It is of note that contaminant taxa were far more predominant in the sorted IgA + and IgA − fractions. This may be due to the lower biomass in these samples but could also be due to species being introduced during sorting or library preparation. These were accounted for by additional screening to only retain species observed in the pre-sort sample for each mouse, which led to no unexpected taxa being observed in any sample (Fig. S4C). Thus, we recommend using abundance thresholding followed by validation against pre-sort samples to filter taxa prior to scoring in IgA-Seq analyses.
After quality control, Hh was detected in the whole and IgA + fractions of all four Hh+aIL10R mice and samples clustered clearly by both experimental condition and sort fraction (Fig. 5e, f). Hh was also detected in the pre-sort sample in one of the aIL10R only mice, but as this was not detected in either the IgA + or IgAfractions, it ultimately did not influence the IgA binding scores for this mouse. We compared the different IgA binding scores across all mice, focusing on the methods that included both the IgA − and IgA + fractions. All methods broadly ranked species IgA binding in the same order and there was a significant positive correlation between all scores (Fig. S5A). The Kau index and probability ratio shared the strongest correlation, with the main difference between the two being a more even distribution of values within the probability ratio than in the Kau index, which contained isolated peaks of scores where taxa were predominantly observed in the IgA + or IgA − fraction (Fig. S5A).
When comparing species IgA binding scores between groups, the method used had important and distinct effects on the interpretation of the results. All methods identified significant differences in IgA binding of Clostridium between the Hh+aIL10R and aIL10R only treatments (Fig. 5g, Fig. S5B). This was the largest effect observed and was clearly detected by all three methods in the same direction with an increased IgA binding of Clostridium in colitic mice. Some taxa, such as Enterococcus, were only observed in the IgA + and/or IgA − fractions in the colitic mice. It was therefore not possible to compare their binding relative to the control mice. This highlights the need to additionally observe presence/absence patterns of taxa in the fractions alongside using quantitative scoring methods.
The Kau index and probability ratio additionally detected a significant difference in IgA binding to A. muciniphila between the groups. However, the Kau and Palm scores were lower for A. muciniphila in colitic mice, whereas the probability ratio was higher (Fig. 5g). This discrepancy results from the large contribution of Hh to the relative abundances in the IgA + fraction of the colitis group. This inherently reduced the observed relative abundance of A. muciniphila in this fraction even in the absence of a change in the absolute levels bound to IgA. This is less pronounced in the IgA − fraction (where less Hh is observed) and thus the overall IgA + /IgA − ratio of A. muciniphila was reduced (Fig. 5g). However, even though A. muciniphila accounts for a reduced percentage of the IgA + fraction, the total IgA + fraction is much larger in the colitic mice (Fig. 5h). A smaller percentage of this much larger total IgA + fraction actually represents a larger total amount of A. muciniphila bound to IgA in the colitic group (Fig. 5h). Thus the probability ratio, which accounts for the fraction sizes, infers a positive rather than negative shift in IgA binding of A. muciniphila. This highlights the benefit of using probabilitybased metrics that provide a direct estimate of the likelihood of IgA binding for each taxon in comparison to solely abundance-based indices, which might lead to a conflicting interpretation without careful consideration of what the scores represent.
The probability ratio increases power to detect differences in IgA binding between conditions As an additional test of the probability ratio, we applied it to a published IgA-Seq dataset with extensive experimental evidence to support differences in IgA binding between conditions. Huus et al. used FACS-based IgA-Seq to study bacterial binding of GI IgA in a mouse model of nutrient restriction [20]. Mice on a conventional (CON) or moderately malnourished diet (MAL) were tracked from weaning to adulthood. Applying the Kau index, they found that mice on a conventional diet develop increased binding of Lactobacillus as they age from three to seven weeks old, whereas malnourished mice do not. They observed this effect to be consistent across different sites in the GI tract and confirmed it was driven by differential binding of IgA to Lactobacillus using an enzyme-linked immunosorbent assay (ELISA). This provided an ideal dataset to test the robustness of the probability ratio.
Huus et al. provided processed 16S rRNA gene sequencing data generated from faecal samples used in their study and the associated IgA + and IgA − fraction sizes, which had been recorded for a subset of the original data [20]. This enabled the probability ratio to be calculated for samples from four CON and MAL mice at each of the three ages. The Kau score (used in the original study) was also calculated as a comparator.
Comparing between diet and ages using the Kau index, we identified five taxa with differential IgA binding (Fig. 6a).
Lactobacillus was the only taxon significantly different by the interaction of age and diet and reflected the patterns observed previously. Three taxa, Lactococcus and unclassified Lachnospiraceae and Erysipelotrichaceae, displayed increased IgA binding in MAL mice. All three of these taxa were observed with the same association in the original analysis of the dataset. We additionally observed a reduced IgA binding of Roseburia in 3week-old MAL mice. This was not reported in the Huus et al. analysis using the Kau index. However, due to the smaller sample size available we explicitly targeted our testing to between group effects and the prior analysis did not.
Carrying out similar modelling using the probability ratio to score IgA binding, we reproduced all of the associations identified using the Kau index with the exception of the unclassified Erysipelotrichaceae taxon, which was no longer significantly different between groups  (Fig. 6b). Across all of the other associations detected as significant by either the Kau Index or probability ratio, the probability ratio consistently produced lower p values (Fig. 6c, left). This reflects an increased magnitude in the differences between the mean scores of the MAL and CON groups in most comparisons when using the probability ratio compared to the Kau index (Fig. 6c,  right). This can in part be attributed to a reduced variation in taxon scores within groups when using the probability ratio (Fig. 6d), as was also observed in the simulated data (Fig. 3e). This increased power is an additional benefit of accounting for the size of the IgA + and IgA − fractions when scoring IgA-Seq data.

Discussion
The development of IgA-Seq has enabled powerful, unbiased investigation of immune interactions between host and commensal microbiota. As this technique gains wider adoption, it is critical to benchmark experimental and analytical approaches to ensure accurate quantification and interpretation of IgA targeting. To this end, we provide an optimised IgA-Seq protocol and novel scoring approaches that facilitate accurate identification of IgA-coated bacteria. We developed a simulation platform that enabled the first formal testing of IgA-Seq scores. This revealed that existing scoring approaches are influenced by relative abundances and, without careful consideration, can lead to the inference of alternate IgA binding patterns solely as a result of pre-sort abundance shifts. Instead, we have developed new probability-based scoring methods that overcome these limitations by directly quantifying the likelihood of a taxon being bound by IgA. Importantly, as we have showed in both simulated and in vivo datasets, these probability-based scores provide increased power to detect differences in IgA binding between groups and promote accurate interpretation of IgA targeting of taxa in health and disease states.
Analysis of in silico and in vivo datasets revealed that it is critical to properly account for the pre-sort taxonomic composition of a sample given its effect on the relative nature of post-sort fraction abundances when comparing IgA targeting between case and control groups. Scoring indices that neglect to do so can infer artificial and/or conflicting shifts in IgA binding estimates between experimental groups. Simulating consistent group-wise shifts in pre-sort abundances resulted in significant differences in taxon level binding estimates when using the Kau index. It is not that the Kau scores are incorrect, but rather that they represent an abstract combination of both the IgA binding and abundance of the taxon scored, and the IgA binding and abundance of all other taxa in the sample. This inherently complicates their interpretation. In contrast, the probability-based scores adjust the taxonomic abundances by the number of bacteria in IgA + and IgA − fractions. This overcomes the relative nature of quantification, which removes the influence of other taxa, and allows for proper adjustment for pre-sort composition. As a result, probability-based scores provide a direct measure of the likelihood that a bacterium is bound to IgA given that it belongs to a specific taxon. As shown with the increased binding of A. muciniphila in the colitis model, this simplifies the interpretation of between group changes in IgA binding when using the probability ratio.
Aside from artefacts that arise in between-group comparisons, all scoring indices were correlated when considering global patterns across samples. This suggests that taxa ranked as generally highly bound in previous studies are unlikely to change if quantified using the novel probability-based scores. However, we found that the influence of pre-sort taxonomic composition meant there was higher variability between samples when using the Kau and Palm scores compared to the probabilitybased methods. This increased variance reduces the power to detect effects in experimental data and resulted in a moderate increase in power when applying the probability ratio to the Huus et al. data. Maximising power by using probability-based metrics is particularly desirable given the high investment of experimental time required to generate IgA-Seq data, which can introduce a practical limit on sample sizes.
Our work highlights the importance of selecting scores based on the underlying biological question.
If the primary interest is to simply identify which taxa have the most bacteria bound by IgA, then the IgA + abundance of a taxon can be used directly. If the purpose is to identify which taxa are being more/less preferentially bound over other taxa, which is arguably more relevant when studying IgA targeting, then probabilitybased scores that account for the relative nature of taxonomic abundances should be used (Fig. S6). To select between using the IgA + probability or probability ratio: if the sole interest is to identify taxa most likely to bind IgA then the IgA + probability is sufficient and benefits from its direct interpretability as a probability of IgA binding (ranging from 0 to 1); otherwise, the probability ratio provides a more generally applicable score that also provides better resolution of relative binding between taxa with lower IgA binding. For example, detecting the loss of Lactobacillus binding in undernourished mice detected in the Huus et al.'s study [20]. The probability ratio additionally benefits from being centred on zero, so a positive or negative score indicates if the majority of a taxon's bacteria fall into the IgA + or IgA − fraction respectively. Score choice might also be based on available sequencing data. The IgA + probability requires sequencing of the IgA + fraction and the pre-sort sample and the probability ratio requires sequencing of the IgA + and IgA − fraction. However, sequencing of both the IgA + and IgA − fractions and the pre-sort sample has often been applied in previous studies [17,18]. As we have shown here, this provides additional benefits in the detection of contaminants introduced by sorting and overall compositional differences between experimental groups.
We have also presented an optimised IgA-Seq protocol that enables precise separation of IgA + and IgA − gut microbiota, and their cultivation to support subsequent interrogation in gnotobiotic disease models and/or strain phenotyping. Our FACS-based protocol outperformed MACS in sort purity and benefits from innate quantification of the IgA + and IgA − fraction sizes required to calculate the probability ratio. We employed a protocol design strategy that combined individual approaches applied in previous studies including the controls used for validation of the staining and gating strategy [17,18]. Whilst access to specific flow cytometers may vary, the concepts in experimental design presented herein can be applied to validate IgA-Seq protocols on different sort modalities. To minimize technical variation between experiments, regular instrument calibration alongside staining controls and gating templates should be employed.
This framework can additionally be used to evaluate reagent selection, a process essential to obtaining robust IgA-Seq results. To differentiate bacteria from intestinal debris, we utilised a sensitive SYBR I dye that efficiently stains diverse bacteria and has been used in other IgA-Seq studies [30,31], whilst others have used SYTO BC dyes [12,18,20]. More comprehensive studies are needed to elucidate the potential for taxon-specific false positives and negatives that occur with the usage of different DNA dyes in different environmental contexts. We selected anti-mouse IgA antibody clone mA-6E1 for our studies, which has been used by Palm et. al., Huus et. al. and a large number of other IgA-Seq studies [6,29,[37][38][39][40]. We utilised the same antibody lot for all experiments within this study, but staining controls and gating strategies need to be verified in instances of lot-to-lot variation. Other studies have also used polyclonal antibodies and biotin-streptavidin staining approaches [12,18]; future work applying similar benchmarks will help shed light on how these may differ in specificity and/or affinity and the subsequent impact on IgA-Seq results. The use of germ-free and Rag1 −/− mice as negative controls and our novel positive control approach to protocol evaluation, isolating pure IgA-bound B. theta from a complex microbiota, provide benchmarking standards to validate these various reagents in the context of specific experimental setups.
16S rRNA gene sequencing has enhanced our ability to identify microbes in a variety of environments, including low biomass samples for IgA-Seq. A major challenge to accurately characterize low biomass communities is contamination from external sources, as we observed in the enrichment of low abundance contaminant reads in our sorted samples. Contamination could potentially occur during sorting, DNA extraction and library preparation [41]. Whilst our negative FACS sheath controls were blank, environmental contamination could still occur whilst sorting IgA + and IgA − fractions as sample and collection tubes remain open for the duration of the sort. In our in vivo model, using a library prepared from a blank kit, we were able to trace some but not all contaminants to the DNA extraction step, a well-established source of contaminants [41]. In addition, index hopping and stochastic sequencing effects have been shown to contribute to noise in low biomass samples on the Illumina MiSeq platform [42,43]. In our study, even utilising a defined MM12 community along with the inclusion of sheath and extraction controls, there were inevitably low abundance contaminants that needed screening from the sorted fractions. This highlights the need to carryout rigorous controls for the low biomass sequencing runs. In particular, we found it invaluable having the pre-sort sample sequencing available as a reference to remove taxa not observed in this higher biomass source.
Similar approaches to IgA-Seq have also been used to study other isotypes, such as IgM and IgG [22,[44][45][46], and immunoglobulins isolated from other sites, for example quantifying GI bacteria bound by serum IgA [29]. Expanding beyond IgA targeting of the gut microbiota, in future work, our protocol design strategy and scoring approaches could be applied to investigate other sample types and immunoglobulins. Further work could also be done to explore other sequencing variables influencing IgA binding scores. For example, amplicon-based taxonomic abundances are inherently influenced by PCR biases, differences in 16S rRNA gene copy number and the taxonomic resolution of the variable region amplified [47][48][49]. However, these are more general limitations of 16S rRNA gene sequencing and their impact should be minimised assuming they occur uniformly in the IgA + and IgA − fraction. Mock sorts have previously shown that the sorting process itself has little influence on community structure [17].
Technical variation in the isolation of IgA + bacteria will also influence IgA binding estimates. We have not explored the influence of sort purity on binding estimates. We assume that using a less pure sorting approach such as MACS would reduce the accuracy of IgA binding indices across the various scoring approaches. This is not likely to influence study outcomes if sort purity is uniform across an experiment but it may be valuable to ensure that this is the case; in particular to check that sort purity is not correlated with experimental conditions. Similarly, using a more sensitive anti-IgA antibody or setting wider IgA + gating in FACS will increase the size of the IgA + fraction and inherently increase the likelihood a taxon will be IgA + . This will not influence relative scoring between taxa within a sample but it is important to maintain these variables across samples in an experiment and/or consider randomisation of samples across technical batches to avoid confounding of experimental groups. For similar reasons, absolute score values from any method are unlikely to be analogous between experiments; however, the relative rankings of taxa should remain comparable. IgA-Seq is additionally limited to scoring relative, overall, IgA binding between taxa. Understanding what drives the observed differences, such as alterations in the spatial distribution of bacteria or changes in the affinities of host IgA, requires additional targeted experiments, such as fluorescence in situ hybridisation and enzyme-linked immunosorbent assays [10,20].
We have developed an optimised experimental protocol and novel scoring methodologies for IgA-Seq. These enable robust and accurate identification of commensal bacteria targeted by the host immune response. The IgAScores R package facilitates the application of these scores and we have presented a guiding framework for selecting appropriate metrics. Furthermore, its simulation platform will aid the future development of IgA-Seq analytics. This work enhances the accessibility and quality of IgA-Seq studies, which should consequently uncover novel mechanistic insights into host-microbiota interactions and the microbial drivers of health and disease.

Conclusion
Existing methods used to score taxon-level IgA binding in IgA-Seq experiments are influenced by the initial taxonomic composition of samples and the relative nature of abundance quantification. This inflates the variability in these scores and can have important consequences in between-group comparisons, where scores might infer false or contrasting associations as a result of systematic differences in microbiome composition between groups. Probability-based IgA scores can be used that overcome this by providing a direct estimate of the likelihood a bacterium will be bound by IgA given its taxonomy. These scoring methods are available in the IgAScores R package. Combined with the practical approaches described here, these analytical developments enhance our ability to study taxon-level immunoglobin binding of commensal microbiota; by increasing the power to detect associations through reduced variability across samples and by providing an accurate, direct, quantification of taxon binding from IgA-Seq data.

Methods
Indices for scoring relative IgA binding of taxa in IgA-Seq experiments Palm index In Palm et al.'s initial description of IgA-Seq, they simply used the ratio of a taxon's IgA + to IgA − abundance in a score termed the IgA coating index [17]. Here, for clarity, we refer to this as the Palm index (Palm). If IgA + and IgA − represent taxon abundances in the respective fractions, for taxon i in sample j the Palm index is defined as in Eq. 1, where if IgA − ij is equal to zero a pseudo count (c) is added to negate division by zero. Throughout, we use a pseudo count of a similar magnitude but below the smallest abundance observed in any sample.

Kau index
Kau et al. described the IgA index, which is also solely based on IgA + and IgA − relative abundances [18]. It uses log transformation to make scores more comparable between taxa and considers the difference of the IgA + and IgA − abundances relative to their total. This has the benefit of centering the score such that positive and negative values reflect overall higher abundances in the IgA + or IgA − fractions respectively. Using notation as previously, the Kau index (Kau) is defined as in Eq. 2.

IgA+ probability
The relative abundances of taxa in a fraction or pre-sort sample sum to one. These can be considered as probabilities. For example, taxa abundances in the pre-sort sample represent the probability that any given bacterium in the sample belongs to a given taxon. The abundances in the IgA + and IgA − fractions can be considered as conditional probabilities (i.e. the probability of event A given that B has occurred). For instance, IgA + abundances represent the probability a bacterium belongs to a taxon given that it is sufficiently bound to IgA to be in the IgA + gate. In the case of IgA-Seq, we are interested in the inverse: the probability that a bacterium will be IgA + if it belongs to a given taxon. This is a posterior probability and can be calculated following Bayes' theorem (Eq. 3).
Where P(A| B) is the posterior probability that a bacteria is IgA + given that it is from a chosen taxon. P(A) is the probability that any bacterium is bound to IgA sufficiently to be in the IgA + fraction-measured as the percentage of bacteria that fall within IgA + FACs gate. P(B) is the probability that a bacterium will belong to a given taxon-measured as the taxon's relative abundance in the pre-sorted sample. P(B| A) is the probability a bacteria belongs to the taxon given that it is in the IgA + fraction-measured by the taxon's IgA + abundance. In which case, the posterior probability (Prob IgA þ ) for taxon i in sample j is defined as in Eq. 4.1. Where FracSize equals the fraction of sorted bacteria in the IgA + gate and PreSort equals taxon abundances in the pre-sorted sample.
This measure will be referred to as the IgA + probability. The numerator provides a direct estimate of the proportion of initial bacteria bound to IgA avoiding the relative influence of other taxa in the IgA + fraction abundance. The numerator must be normalised to the total amount of that taxon in the original sample prior to sorting (PreSort) to make it comparable across taxa. As this is the pre-sort level of the taxon, the denominator should theoretically always be larger or equal in value. However, due to amplification and other technical biases, this may not practically be the case in IgA-Seq data, thus Eq. 4.2 can be used to account for this and ensure probability estimates fall below 1.
This will assign a probability of 1 when the estimate of IgA + bacteria in a taxon is larger than the observed total abundance in the pre-sort sample, a reasonable assumption being that very high levels of IgA binding would be necessary to produce these results.
This score only requires two additional pieces of data than the Kau and Palm indices: the taxonomic abundances in the pre-sort fraction and the percentage of sorted bacteria that are in the IgA + gate in the FACS run. The latter value is already generally recorded during sorting and most IgA-Seq studies to date have profiled the microbial communities from pre-sorted sample s [17,18]. Furthermore, if the only interest is to identify highly bound taxa, this method can negate the need to sequence the IgA − fraction.

Probability ratio
The IgA + probability only provides a measure of the likelihood a taxon is sufficiently bound to IgA to be in the IgA + fraction. The converse, the IgA − probability, can also be calculated by simply substituting the IgA − abundances and fraction sizes into the IgA + probability equation. The IgA − probability thus represents the probability a bacterium is in the IgA − fraction given it belongs to a given taxon. If we calculate both the IgA + and IgA − probabilities from a sample, we can be combine the estimates by comparing the ratio of the IgA + probability to IgA − probability, a score we term the IgA probability ratio (Eq. 5). This provides finer granularity of estimates of the relative binding between taxa with low levels of IgA binding than using the IgA + probability alone. Additionally, as both the IgA + and IgA − probabilities are normalised to the taxon's abundance in the pre-sort fraction, this can be removed from the equation (Eq. 5) negating the need to sequence the pre-sort sample whilst still accounting for pre-sort compositional differences.
Pseudo counts (c) must be included in the probability ratio to negate division by zero where taxa are not observed in the IgA − fraction. Using the probability ratio has the benefit of capturing both extremes of IgA binding within a single score (whereas the IgA + or IgA − probabilities are limited to either direction). Additionally, we take the log of the ratio, which, similar to the Kau index, centres values on zero. This directionality provides a clear indication of whether bacteria in a taxon have a higher probability of being found in the IgA + or IgA − fractions.
One limit of using the ratio as in Eq. 5 is that the range of possible values returned is determined by the pseudo count used. The magnitude of these effects can be reduced by dividing ProbRatio by a scaler based on the pseudo count as in Eq. 6. This will bring all values within the range − 1 to 1, with these values representing the extreme case where all bacteria are in a single fraction and belong to a single taxon. Scaled probability ratios are used throughout.
The Palm, Kau, IgA + probability and probability ratio scores were all calculated as defined above using the igascores function that we have created as part of an R package IgAScores, which is available at https://github. com/microbialman/IgAScores.

Simulating an IgA-Seq experiment
A simulation was used to generate IgA-Seq data where the underlying relative taxonomic binding of IgA binding was known. Ten species were assigned a mean IgA binding value by random sampling from an exponential distribution (raising two to the power of the generated values to increase their spread). These values are arbitrary numbers that simply provide a relative IgA binding measure between the species and are used as the mean to generate a normal distribution of IgA binding for each species (with a standard deviation of one). For initial simulations, 30 samples were then generated with random pre-sort abundances for the ten species sampled from a log normal distribution.
To simulate an IgA-Seq experiment, each sample was assumed to have profiled 100,000 bacteria (a reasonable magnitude for 16S rRNA gene sequencing). For each sample, the 100,000 counts were assigned to one of the ten species based on the sample's species abundances. Then, for each of the 100,000 bacteria in each sample, a value of IgA binding was generated at random from the normal distribution defined by its species mean. The distribution of IgA scores across all species was used to identify values that cut off the upper and lower tails of IgA binding and these were used to define the bacteria within the IgA + and IgA − fractions respectively. The bacteria within the IgA + and IgA − fractions for each sample were then used to generate the IgA + and IgA − relative abundances and all other measures required for generating the various IgA-Seq scores.
A secondary simulation was performed to model the case where two groups are compared with a consistent difference in starting abundances between them. For this, the simulation was run using the same ten species as before but now considering 60 samples. Additionally, prior to generating bacterial IgA binding values, 0.5 was added to the abundance of species ten in half of the samples before renormalizing the abundances to sum to one. Samples with additional species ten were considered cases and the unaltered samples considered controls.
All stages of the IgA-Seq simulation were carried out using the function simulateigaseq, which is also part of the IgAScores package. All IgA-Seq indices were calculated as defined previously using the igascores function. For simulated data experiments, a value of 10 −6 was used for pseudo counts when required as this was the same magnitude as the minimum non-zero abundance value observed in the simulated data. Correlations were assessed using the cor function in R [50]. Between group comparisons were performed using the wilcox.test function, adjusting p values using p.adjust. Code for all simulation analyses can be found at https://github.com/ microbialman/IgAScoresAnalyses.

Mouse strains
All experiments were conducted in accordance with the UK Scientific Procedures Act (1986) under a Project License (PPL) authorized by the UK Home Office. Animals were housed in accredited animal facilities at the University of Oxford and provided with sterile water, food ad libitum and environmental enrichment. Both males and females were used in experiments after 8 weeks of age. Specific pathogen-free (SPF) C57BL/6J and C57BL/6J Rag1 −/− mice were routinely tested and negative for Helicobacter spp. and other known intestinal pathogens. Germ-free (GF) C57BL/6J mice were maintained in sterile flexible film isolators and routinely screened for bacterial contamination by culturing, Gram stain and 16S rRNA gene PCR. GF C57BL/6J mice were stably colonized with a defined consortium of 12 bacterial members, see below, of murine intestinal microbiota to generate a sDMDMm2 [35,51] (MM12) colony that was maintained in sterile flexible film isolators. All GF and MM12 mice were maintained on sterile food and water.

Helicobacter hepaticus-induced colitis
To induce colitis using Helicobacter hepaticus (Hh), NCI-Frederick isolate 1A strain 51449 was grown as previously described [52]. Littermate control mice with MM12 microbiota were randomly assigned to groups. Mice were fed 1.0 × 10 8 CFU by oral gavage on day 0 and 1. On day 0 and 7, 1 mg IL-10 receptor blocking antibody (clone 1B1.2) was delivered by intraperitoneal injection. Control mice were housed in a separate flexible film sterile isolator and were given IL-10 receptor blocking antibody in the absence of Hh administration.
Animals were routinely monitored for weight loss and onset of diarrhoea. Mice were sacrificed at day 14 and processed aseptically in a class II microbiological safety cabinet. Colon contents were collected for IgA-Seq and lipocalin-2 measurements (EMLCN2, Life Technologies) and stored at − 80°C. Formalin-fixed paraffin-embedded cross-sections of proximal, middle and distal colon were stained with haematoxylin and eosin to assess histopathology. Sections were examined for hallmarks of inflammation: epithelial hyperplasia, goblet cell depletion, leukocyte infiltration, crypt abscess formation, submucosal leukocyte infiltration and interstitial edema.

Lipocalin-2 enzyme-linked immunosorbent assay
To prepare intestinal supernatants, 50-100 mg of colon contents were added to screw-capped O-ring tubes filled with 1 ml of cOmplete ULTRA EDTA-free protease inhibitor (Roche). Samples were homogenized on the Tis-sueLyser II (QIAGEN) for 10 min at 30 Hz and debris separated by centrifugation (10,000×g for 5 min at 4°C). Supernatant was carefully collected and either used immediately or stored at − 20°C. Lipocalin-2 was quantified using a mouse LCN2 enzyme-linked immunosorbent assay (ELISA) kit (ThermoScientific) according to manufacturer's instructions. Each sample was run in duplicate and measurements reported as ng of lipocalin per gram of sample.

Preparation of samples for bacterial flow cytometry
Only sterile reagents were used to process samples. Frozen stool and intestinal contents (10-50 mg of sample) were rehydrated in 1 ml of PBS and suspended by vortexing. Large debris were separated by centrifugation (50×g, 15 min, 4°C) and the clarified supernatant was passed through a 70 μm filter. A 500 μl aliquot of the suspension was transferred into a new 1.5 ml microcentrifuge tube and an additional 1 ml of PBS added. Bacteria were pelleted (8000×g, 5 min, 4°C) and blocked with 100 μl of 10% rat serum (10 min, ice). An additional 1 ml of PBS was added to the blocking suspension and centrifuged, and pelleted bacteria were washed with 1 ml of FACS buffer (PBS + 1% bovine serum albumin). Following centrifugation, the resulting pellet was stained with PE-conjugated rat anti-mouse IgA (1:100, eBioscience clone mA-6E1) in a 100 μl volume (30 min, covered on ice). Clone mA-6E1 was titrated testing concentrations of 1:10, 1:20, 1:50, 1:100, 1:200, 1:500 and 1: 1000. Samples were washed twice with FACS buffer and resuspended in 1 ml buffer containing SYBR Green I nucleic acid gel stain (1:400,000, Life Technologies) in the dark at room temperature. The optimal concentration of SYBR was determined by titration testing a range from 1:100,000 to 1:500,000. Samples were then placed in a covered ice box and immediately acquired on the flow cytometer (LSRFortessa X20 and FACSAria III, BD Biosciences).
To validate the IgA-Seq protocol, 1 ml of an overnight B. theta culture was first stained with mouse IgA clone 255.4 (25 min, ice), washed with 1 ml PBS, followed by staining with PE-conjugated rat anti-mouse IgA (1:100, eBioscience clone mA-6E1). Stained bacteria were washed three times in FACS buffer before resuspension in 1 ml FACS buffer. A 500 μl aliquot of IgA-coated B. theta was spiked into a faecal suspension prepared from specific pathogen-free (SPF) C57BL/6J mice. SYBR was added to the final mixture.

Fluorescence-activated cell sorting of bacterial populations
Samples were sorted using a FACSAria III (BD Biosciences) instrument with external aerosol management. Both sample and collection tubes were maintained at 4°C . Sterile, preservative-free PBS without calcium or magnesium (Invitrogen) was used for sheath fluid and was autoclaved prior to sorting. The flow cytometer was sterilized according to the manufacturer's Aseptic Sort protocol to decontaminate the entire sheath and sample paths. In addition, both before and after sorting, freshly made 10% bleach was run for 20 min on max flow rate followed by sterile PBS for 10 min. Sheath fluid from the droplet stream was collected prior to sorting as a control. Sample Line Backflush was run between samples to prevent carryover.
For finer resolution of small particles including bacteria, the neutral density filter was removed and threshold set to 200 for side scatter (SSC). Drop delay values were determined with the Accudrop feature and were updated with each new sort setup. To sort and culture live B. theta, various nozzle sizes (70, 85, 100 μm) and sheath pressures (20,30,45, 70 psi) were tested. A 10 × 10 grid was programmed into the automated cell deposition unit (ACDU) stage, and B. theta suspended in PBS was sorted in Single Cell precision mode onto BHIS agar plates for enumeration of colony-forming units (CFU). Nozzle sizes and sheath pressures were also tested for sort purity of IgA + and IgA − populations of faecal bacteria.
For stool and intestinal samples, the SYBR nucleic acid stain (FITC channel) helped to differentiate bacteria from autofluorescent debris (BV421 or PerCP Cy5.5 channels). Samples were diluted in fluorescenceactivated cell sorting (FACS) buffer containing SYBR in order to achieve an event rate of 5000-6000 events per second. To begin, 100,000 events were recorded from the initial sample (pre-sort) in order to calculate the percentage of bacteria in the IgA + and IgA − gates used in the probability-based scores. Gates were then drawn to delineate bacterial populations that were bound with IgA (IgA + ) or unbound (IgA − ). The 4-way purity precision mode was used to capture 400,000 events for each IgA + or IgA − fraction. Fractions were collected into autoclaved 2 ml screw cap tubes (ThermoFisher Scientific) pre-coated with 200 μl FACS buffer. Sort purity of each fraction was recorded from 100 events, and DNA was immediately extracted from the remainder. BD FACSDiva and Flowjo software (TreeStar) were used for instrument setup and data analysis.

Magnetic-activated cell separation
As a comparison with FACS, IgA-stained samples were incubated for 30 min at 4°C with anti-PE microbeads (1:50, Miltenyi) in magnetic-activated cell separation (MACS) buffer (0.5% bovine serum albumin and 2 mM EDTA in PBS). Microbeads were titrated by testing concentrations of 1:5, 1:20, 1:50, 1:100, 1:200, 1:500. Cells were pelleted by centrifugation (10,000×g, 5 min, 4°C) and washed twice in MACS buffer before resuspension in 1 ml. LS columns (Miltenyi) were placed on a magnetic stand (Miltenyi) and pre-washed with 10 ml MACS buffer before the sample was loaded. Unlabelled "MACS neg" cells were collected in 3 × 3 ml washes. The LS column was then removed from the magnet and the labelled "MACS pos" cells were collected in 6 ml MACS buffer using the LS column plunger. MACS fractions were then analysed for purity by flow cytometry.

DNA extraction
Only sterile molecular-grade reagents and filtersterilized buffers were used to extract DNA. Sorted fractions were centrifuged (16,000 × g, 20 min, 4°C) and the supernatant removed such that only 300 μl remained in the tube. For unsorted (pre-sort) samples, 100 μl was transferred into an autoclaved 2 ml screw cap tube containing 200 μl FACS buffer. The following was then added to each tube: 250 μl autoclaved 0.1 mm zirconia/ silica beads (Biospec), 300 μl Lysis buffer (200 mM NaCl; 200 mM Tris; 20 mM EDTA; pH 8), 200 μl 20% SDS and 500 μl of phenol:chloroform:isoamylalcohol (PCI 25: 24:1, pH 7.9, Sigma). Samples were bead beaten (6500 rpm for 4 min, 3 min on ice, repeat 2×, Precellys tissue homogeniser), then centrifuged (6000×g, 10 min, 4°C) to separate phases. The aqueous phase was transferred into 5PRIME Phase Lock Gel Light tubes (Quantabio) and an equal volume of PCI added. After mixing by inversion, samples were centrifuged (16,000×g, 10 min) and the aqueous phase transferred to a new Eppendorf tube. DNA was precipitated with 1/10 volume 3 M NaOAc (pH 5.5) and 1 volume of isopropanol overnight at − 20°C. After precipitation, DNA was pelleted (16, 000×g, 30 min, 4°C) and the supernatant carefully removed. The pellet was washed in 500 μl 70% ethanol, centrifuged (16,000×g, 10 min, 4°C) and the wash repeated again. Residual ethanol was removed from pelleted DNA and the tube placed on a 37°C heat block until the pellet was completely dry. Cleaned DNA pellets were resuspended in 30 μl TE buffer at 50°C for 30 min, then stored at − 20°C until 16S rRNA gene library preparation.
All 16S rRNA gene sequencing data were quality checked using FastQC and processed using DADA2 within the NGSKit pipeline (https://github.com/ nickilott/NGSKit) to generate an amplicon sequence variant (ASV) table [55]. In brief, DADA2 was run on each sample individually. The filterAndTrim function was first used on the paired-end files with the following parameters: trunclen 240,140; maxN 0; maxEE 2,2; truncQ 2. These filtered forward and reverse read files were then used for generation of the error models using the learnErrors function over 10,000,000 bases; for dereplication using the derepFastq function; and for sample inference using the dada function. The reads were then merged using mergePairs and chimeric reads screened using removeBimeraDenovo. The assignTaxonomy and addSpecies functions were used to assigned taxonomy to ASVs using a database built from a combination of both the RefSeq and RDP databases [56]. Unless otherwise specified, the default parameters were used for all functions. ASV counts were then summarised at the species level and all analyses carried out on these relative abundances (summing to one within a sample).

IgA-Seq in a mouse model of colitis
ASV counts were generated from 16S rRNA gene sequences as described above, producing 54,091 ± 7970 (mean ± SD) ASV counts per sample (excluding the blank control which had 5119 counts). These were summarised at the species level. Given the samples contained a known MM12 microbiota, the presence of unexpected species in samples and species observed in a sequencing control blank were used to find an abundance threshold that removed the majority of contaminants. Non-metric multi-dimensional scaling was used to visualise the relationship between samples from each mouse in each fraction using the metaMDS function from the vegan R package [57]. Comparisons of IgA-Seq data were performed between the four anti-IL10R only controls and the four Hh-treated mice for the Palm, Kau and the probability ratio scores. These were generated using the igascores function as for the simulated data. A pseudo count of 10 −3 was added where necessary as this was the magnitude of the lowest non-zero abundance value observed in the data.
Scores were compared between the colitic and control groups for taxa where there were at least three samples scored in both groups (taxa were scored NA if they had zero abundance in both the positive and negative fractions). Given the small sample size of each comparison, significant differences were determined by carrying out all possible permutations of group labels and observing the difference in group means for each taxon for each score. The true difference in means was compared to the differences across permutations to derive an exact p-value, considering significance where p < 0.1. This more permissive threshold was used given the small sample size (min three per group). To generate comparable effect sizes between scores, which are on different scales, the difference in means between groups was normalised to the standard deviations within the groups using the strictly standardised mean difference, calculated using the ssmd function from the phenoDist R package [58]. The species level data for these analyses are made available within the IgAScores package and all code for the a n a l y s e s a r e a v a i l a b l e a t h t t p s : / / g i t h u b . c o m / microbialman/IgAScoresAnalyses.

Analysis of Huus et al.'s data
Sample metadata, IgA + and IgA − fraction sizes and the pre-processed taxon count table used in the original study were provided by the authors. The counts table was filtered to remove low abundance taxa using code published with the original study and then subset to just the samples for which the fraction sizes were available. The Kau Index and probability ratio scores were then generated using the igascores function as previously. For each taxon, a two-way analysis of variance (ANOVA) was performed, using the aov function in R, to determine if there was a difference in IgA binding score between diet treatments or the interaction of diet and mouse age. This was carried out for both the Kau Index and probability ratio, considering a taxon significant if either diet or the diet-age interaction had a nominal p < 0.05. Taxa were only tested that had scores (i.e. were not absent in both the IgA + and IgA − fractions) in at least six samples spread across all age points and both conditions. The magnitude of effects was estimated as the difference in mean scores between MAL and CON for each taxon and was calculated at each of the three ages individually. Mean differences were normalised by the standard deviation within groups, using the strictly standardised mean difference as for the colitis analysis (to make the values comparable between the Kau index and probability ratio). Coefficient of variation was calculated within each taxon within each experimental group (age and diet combined) for each scoring method and compared using a Mann-Whitney test. All code for the replication analysis can be found at https://github.com/ microbialman/IgAScoresAnalyses.
Additional file 1: Figure S1. Identification of IgA-coated bacterial populations by flow cytometry. Figure S2. Sorting configuration enables sensitive sort control. Figure S3. Between group comparisons of the Palm index in an IgA-Seq simulation. Figure S4. Quality control to filter spurious taxa prior to analyses. Figure S5. Comparison of IgA-Seq binding scores between groups in a mouse model of colitis. Figure S6. Summary of the representation of the different scores and what questions they address. Table S1. Mapping of bacterial strain names to their mappable identifiers in the 16S rRNA database used.