C57BL/6 (B6) and BALB/c (BC) female mice aged between 6-10 weeks (from the Australian Phenomics Facility, ANU) were used in the study. Animals were housed in a specific pathogen-free environment and used under strict adherence to protocols approved by the institutional Animal Experimentation Ethics Committee (AEEC), ANU, under protocol A2020/39. At experimental end points, animals were euthanized by cervical dislocation according to AEEC approved procedures.

The mammary carcinoma cell lines 4T1 (24) (ATCC), 4T1.2 (25) (kindly provided by Dr Robin Anderson, Olivia Newton-John Cancer Research Institute), 4T1Br4 (26) (kindly provided by Dr Normand Pouliot, Olivia Newton-John Cancer Research Institute), and AT-3-OVA (27) (kindly provided by Dr. Di Yu), the colorectal carcinoma cell lines CT26 (28) (ATCC) and MC38 (29) (kindly provided by Dr. Di Yu), and the melanoma cell line B16-F10 (30) (kindly provided by Dr Christopher Parish) were used in this study. Cell lines were confirmed clear of specific pathogens by Cerberus Sciences (ISO 9001 Licence No. AU843-QC). Cell lines were cultured and subcultured as described previously in supplemented (sRPMI) RPMI-1640 (11875093, ThermoFisher Scientific) (18).

Tumor cells (1 x 10⁵) were injected subcutaneously in the right hind flank (primary tumor) and then 3 days later in the left hind flank (secondary tumor) of syngeneic mice (cell lines 4T1, 4T1.2, 4T1Br4, and CT26 injected in BC mice, and cell lines AT-3-OVA, MC38 and B16-F10 injected in B6 mice) mixed across housing cages as described previously (18). Tumors were grown and monitored for 17-21 days. At study endpoint, mice were humanely sacrificed, and their tumor and spleens excised and weighed. The tumor burden was defined as the pooled tumor weights of each individual at end point.

Blood and spleens from mice were collected at experimental end point. Blood was collected and processed as described previously (18). Spleens were harvested and processed to single cells as described previously with the exception that the red blood cell lysis step was not performed (31).

Spleen cells from 3 (of the 3-5) replicate mice bearing the largest mass of 14-17-day-old tumor from either 4T1, 4T1.2, 4T1Br4, AT-3-OVA, CT26, MC38 or B16-F10 cell lines, or no-tumor controls (host B6 or BC mice) were pooled into 9 separate tubes in a total of 10 milliliters (mL) of phosphate buffered saline (PBS). The 9 spleen cell groups were adjusted to the equivalent of 2 spleen masses based on spleen weights (~equivalent to the mass of 2 normal spleens) by removing appropriate volumes of cell suspension from each tube. Cells volumes were then made to 10 mL with PBS and cells sedimented by centrifugation (300 x g for 10 min), supernatant aspirated and cells resuspended in 2.9 mL sRPMI. Each spleen cell group was then barcoded separately with a unique concentration of carboxyfluorescein diacetate succinimidyl ester (CFSE) and/or cell trace violet (CTV) ( Table S1 ) and all groups were then pooled into one sample as previously described (19). Cells were then suspended in 10 mL of sRPMI and passed through a 70 μm filter mesh and counted. A total of 400 x 10⁶ leukocytes was then suspended in 10 mL sRPMI, passed through a 70 μm filter mesh, sedimented by centrifugation (300 x g for 10 min) and supernatant aspirated, ready for immediate backbone antibody labelling.

Barcoded pooled cells were resuspended in 0.6 mL of Labelling Buffer (PBS containing 5 mM EDTA, 1% BSA [weight/volume]) with 5 μg/mL TruStain FcX™ (anti-mouse CD16/32) antibody (101320, Biolegend) for 15 min at 4°C. Samples were then incubated with a backbone panel of antibodies ( Table S2 ) by adding 0.6 mL of Labelling Buffer with 10% (v/v) Brilliant Stain Buffer Plus (BD) containing a 2x stock of each antibody ( Table S2 ), for 30 min at 4°C. The pooled barcoded and backbone antibody-labelled cells were then resuspended to 13.4 x 10⁶ cells per mL (i.e. 1x10⁶ cells/75 μL) in Labelling Buffer and passed through a 70 μm filter mesh ready for aliquoting into the wells of the LEGENDScreen plates.

A LEGENDScreen Mouse PE Kit (BioLegend) was used for spleen leukocyte screening for cancer-specific cell-surface markers. Plates from the kit were prepared according to the manufacturer’s instructions with lyophilized antibodies in each well of the assay plates being resuspended in 25 μL of deionized H₂O (dH₂O). The pooled barcoded and backbone antibody-labelled cells were added at 75 μL (i.e. 1x10⁶ cell) to each well containing the reconstituted antibodies and incubated in the dark for 30 min at 4^°C. Cells were then washed in Legend Screen Wash provided in the kit, pelleted and resuspended in 40 μL Labelling Buffer containing 1 μg/ml of the viability dye Hoechst 33285 (Invitrogen) and the equivalent of 500 Flow-Count Fluorospheres (7547053, Beckman Coulter) per 40 μL and stored at 4°C overnight before flow cytometry.

Blood samples (5 μL) were labeled with antibodies that included the backbone panel and the screen-identified antibodies ( Table S2 ) and prepared for flow cytometry analysis using methods described previously (18).

Flow cytometry was performed on a BD LSR-II X-20 (BD Bioscience) flow cytometer with FACSDiva software (version 8.0.1). Application Settings were applied to standardize fluorescence intensity readings between experiments, and fluorescence intensities monitored and calibrated using Sphero™ 8-peak Rainbow Beads (110620, BD Bioscience). Voltages were initially set up using unlabeled RBC-lysed blood leukocytes. BD CompBeads (552843, BD Bioscience) were used as compensation controls as previously described (31). Blood cell samples were acquired until a total of 2000 Flow-Count Fluorosphere beads were collected based on side scatter (log-scale) and forward scatter (linear-scale) plot gating. LEGENDScreen samples were acquired at 10,000 events/second using the sample fine adjust and on a low sample flow rate to collect a total of ~1-3 x 10⁵ live (Hoechst 33285⁻-gated) CD45⁺ cells. Every 36^th sample acquisition was followed by a 3 min run on a high sample flow rate with 10% sodium hypochlorite solution in dH2O, then a 2 min run on a high sample flow rate with dH₂O, and the stability of fluorescence signal of each channel assessed by acquiring 5000 Sphero™ 8-peak Rainbow Beads. Raw Flow Cytometry Standard (FCS) files (i.e., FCS3.1) of the data are available upon request at the ANU DATA COMMONS repository (https://dx.doi.org/10.25911/zrp3-nd51).

Flow cytometry analysis was performed using FlowJo v10 software (BD Bioscience), CytoExploreR version 2.0.0 (32), and cytoverse suite of R packages (33, 34). A combination of manual gating and unsupervised Pairwise Controlled Manifold Approximation Projection (PaCMAP) (35) analysis was used to delineate cell populations and assess for manual gate cell population segregation, and cell groups were then named based on marker expression represented by median fluorescent intensities (MedFI) of each marker plotted using heat map dot plots made using the tidyverse suite of R packages (36) (see Results section).

Supervised machine learning was performed using Orange 3 software. Random Forest (38) and CATboost (39) modelling used 100 trees for prediction and 500 trees for ranking the feature importance, with the maximum tree depth of 4 (for Random Forest) and 6 (for CATBoost). For Random Forest, the maximum number of features considered at each node was 5 and subsets smaller than 5 were not split. For CATBoost learning, the learning rate was 0.3, the regularization was lambda 3 and subsampling was 1. For classification of groups using monocytes, CATBoost was used and trained on 66% of randomly sampled data and tested on the remaining 34%. This was repeated 100 times and the results of predicted and actual classes were displayed as a confusion matrix. Feature ranking was done using both Random Forest and CATBoost (built into the models in Orange 3 software). For the learning curve as a function of decreased features (populations), CATBoost was used and trained on 66% of randomly sampled data and tested on the remaining data. This was repeated 100 times and the results were assessed using area under curve of the receiver operating characteristics (AUC; to assess separability of the classes), classification accuracy (CA; proportion of correct classification), precision (ratio of correct positive prediction to all positive prediction), recall (ratio of correct positive prediction to actual positive), and F1 score (weighted average of precision and recall). Final CATBoost prediction on the optimized feature subset was trained on 66% of randomly sampled data and performed on the remaining data.

To compare the means between non-tumor controls (Nil), CT26- and 4T1-burdened cohorts, data were transformed using the formula Y=Log(Y+1) to help normalize distributions and equalize variance, and then assessed by 2-way ANOVA using GraphPad Prism software. Analysis was corrected for multiple comparisons using the two-stage step-up method of Benjamini, Krieger and Yekuyieli (40). False discovery rate of 0.05 and p-values were reported to test the null hypothesis that the means are equal or distributions were from the same population. The overlap R package was used to calculate distribution overlap (41). PaCMAP used the pacmap python package through CytoExploreR. Multidimensional scaling was performed using the cmdscale function from the stats (v3.6.2) R package to visually inspect the improvement in group separability when running ML models using the highest ranked populations. Heatmap dot plots were generated through ggplot2, ComplexHeatmap and HeatmapR (42) R packages. Log ratio (M) log average (A) (MA) plots were constructed using the ggpubr, ggplot2 and ggrepel R packages. Pythagorean trees and confusion matrices were made in Orange 3 software. Circular bar plots were made using ggplot2 in R. The R version used was 3.6.2. Prism was also used for plotting data.

To identify cancer-specific leukocyte surface marker changes, a screening pipeline was established ( Figure 1 ). This involved: a), establishing tumors in mice; b), preparing spleen cells from tumor-bearing mice and their non-tumor-bearing counterparts: c), fluorescent barcoding of the spleen cells from each group with distinct concentrations of CellTrace Violet (CTV) and/or carboxyfluorescein succinimidyl ester (CFSE) vital dyes to give distinct fluorescent groups; d), pooling the barcoded cells and labeling them with a backbone of lineage marker-specific antibodies to identify the main leukocyte populations; e), aliquoting the pooled, barcoded and backboned samples equally into wells of LEGENDScreen assay plates to label them with the >250 PE-conjugated antibodies specific for cell-surface markers for analysis by flow cytometry; and f), assessing the samples for cancer-specific changes via an analysis pipeline.

To assess how the leukocyte subpopulation profiles might change in the presence of cancer, seven murine tumor cell lines representing breast (including, 4T1 (24), 4T1.2 (25), 4T1Br4 (26), AT-3-OVA (27)), colorectal (including CT26 (28) and MC38 (29)) and skin (including B16-F10 (30)) cancers were used to establish tumors in syngeneic strains of mice (B6 and BC), yielding a total of nine groups (seven cancer and two no-cancer control samples). A primary and secondary tumor was established on each hind flank (left and right; to mimic disseminated disease) of each animal for each cell line and grown for 17 days. The endpoint masses showed variable degree of growth, with the tumors on the BC background typically weighing more than those on B6 background ( Figure 2A ). Although most animals had viable tumors, MC38 tumors had negligible mass at the end point ( Figure 2A ). The weight of the spleens extracted from these mice for the LEGENDScreen pipeline also varied considerably between groups, lending support to the theory that certain tumor types may promote spleen enlargement (38) ( Figure 2B ).

After processing the leukocytes through the analysis pipeline, we investigated how the standard lineage backbone antibody labelling characterized the immune cells in the pooled spleens and blood of the animals. The backbone panel comprised primarily antibodies to delineate key leukocyte subsets. Manual gates were created to identify ~24 main cell subsets ( Figure S1 ). The general population segregation was confirmed by overlaying the manually gated populations on graphs of unsupervised Pairwise Controlled Manifold Approximation Projection (PaCMAP) dimensional reduction of all the backbone markers in 2-dimensional space ( Figures S2A, B ). The marker expression of the segregated populations was summarized by dot plot heat maps annotated with pooled cell proportions ( Figure S2C ) and used to confirm population name designations ( Table S3 ). Manual gates were also used to delineate the barcoded groups ( Figure S1 ), and the proportion of CD45⁺ live leukocytes comprising each of the identified leukocyte subpopulations was calculated across the groups ( Figure 3 ). Major changes were seen in the leukocyte composition of all animals with tumors on BC backgrounds in both spleen and blood, whereas more a subtle change in the leukocyte composition was observed in animals with tumors on a B6 background. The major tumor-specific changes appeared to be the proportional increases in myeloid cells, most notably monocytes and neutrophils, particularly in BC mice burdened with 4T1 tumor variants.

To resolve the leukocyte subpopulation differences at higher resolution, we assessed each of the backbone antibody panel-identified splenic leukocyte populations across all tumor groups for changes in the median fluorescent intensities (MedFI) of the >250 PE-conjugated LEGENDScreen antibodies (herein, termed MedFI-PE) relative to background (no tumor) cell levels ( Figure S3 ). This was summarized in a dot plot heat map where an increase in the dot size reflected an increase in the absolute MedFI-PE change in tumor-bearing animals relative to the background control (either increase or decrease). In addition, the plots showed increases and decreases in the relative marker MedFI-PE from the background by the dot color tending toward darker red or darker blue, respectively ( Figure S3 ). This global analysis demonstrated that while most of the cell surface markers showed little change in the presence of cancer, certain leukocyte population markers underwent clear changes across several tumor types ( Figure S3 ).

To narrow the scope of the LEGENDScreen data and inform the next step in the experimental validation, our analysis focused on two cancer types with consistent growth and diverse immune perturbation. These were the CT26 colorectal cancer model, which generated relatively mild immune perturbation, and the 4T1 breast cancer model, which generated marked immune perturbation (18). In addition, markers associated with metastatic potential were assessed using the 4T1 metastatic variant, the 4T1.2 line (25) ( Figure 4 ). To focus on the strongest cancer-associated marker changes, a threshold was set to display markers that were within the top 30% of the increased or decreased MedFI-PE change relative to background no-tumor controls ( Figure 4 ). This was an arbitrary range set to filter for the most extreme changes in surface marker labeling. For quality control, the expression of each marker after subtracting the MedFI-PE of isotype control antibodies ( Figure 4 , left bar plot) was also determined, which revealed that all markers had a positive signal above matching isotope control binding except CD357 and CD28 markers, which had negative fluorescence relative to their isotype antibody controls, suggesting a technical anomaly. CD357 and CD28 markers were thus excluded from further analysis. Among the remaining surface markers, CT26 tumors clearly induced fewer marker changes than 4T1 tumors ( Figure 4 ). While the tumors from 4T1 and their metastatic variant 4T1.2 resulted in a similar number of marker changes, some of these were clearly variant-specific, potentially relating to the metastatic tendency. Since only up to 6 additional antibodies could be added to the backbone antibody panel for further studies given the flow cytometer instrumentation restrictions, the sum of all absolute MedFI-PE changes was calculated across the groups for each marker to identify those with the greatest cumulative changes ( Figure 4 , right bar blot). To ascertain the cell population numbers associated with marker changes, enumeration of each cell population based on % of total CD45⁺ leukocytes and total leukocyte population per spleen ( Figure 4 top annotations) and cell leukocyte population amounts per milligram of spleen ( Figure 4 bottom annotations) were also calculated. These demonstrated that CD24, CD45RB, and CD44 had relatively strong changes across all tumor types, rendering them the preferred candidates for panel inclusion. Also notable was CD62L downregulation on CD8⁺ T cells in 4T1 and 4T1.2, but not CT26, models; CD66a upregulation on myeloid subpopulations in 4T1 and CT26, but not 4T1.2, models; and surface IgD downregulation on B cells in the metastatic 4T1.2 model but not its parental line, 4T1. This combination of markers therefore might reflect the changes that are broadly associated with cancer presence, specific to cancer type and its metastatic potential. While there were several other marker changes in the models, many of these were already in the backbone panel, and others, while potentially of value, were not included in the expanded labelling panel due to restrictions on the total number of markers based on flow cytometry infrastructure. The expression of these markers on the leukocyte populations in which they underwent the most apparent change was also displayed in raw histogram format ( Figure S4 ). This demonstrated that the degree of the changes ranged from subtle to obvious between the CT26, 4T1 and no-tumor control groups.

With the inclusion of CD24-, CD45RB-, CD44-, CD66a-, CD45RB-, and IgD-specific antibodies to the existing leukocyte-profiling backbone panel, we analyzed the blood leukocytes from animals bearing CT26, 4T1 or no tumors across three independent experiments. While the tumor growth was variable, all mice in the CT26 and 4T1 groups had palpable tumors by the endpoint ( Figure S5A ). Manual gating was used to identify leukocyte populations across all experimental groups ( Figures S5B, C ); gated populations overlaid onto a 2-dimensional PaCMAP embedding ( Figures S6A, B ); and phenotype assessed on a dot plot heat map displaying MedFI of marker expression and population proportions ( Figure S6C ). Up to 39 populations were identified with the inclusion of the additional antibodies ( Figures S6B, C, Table S4 ). Separating the animal groups into no-tumor control (Nil), CT26 and 4T1 models ( Figures 5A, B ) showed obvious proportional ( Figure 5B ) and total count ( Figure 5C ) changes in the delineated blood leukocyte populations across the groups. There were obvious fold changes in the counts of numerous leukocyte subsets in both tumor types relative to normal levels, which was more pronounced in animals with 4T1 tumors ( Figure 5C ). The counts of nine populations of myeloid cells were significantly increased by >2-fold above control in CT26 samples ( Figure 5D ); the increases were more dramatic in 4T1-bearing mice with most myeloid cells showing a significant >2-fold increase above no-tumor control and CT26-bearing animals ( Figures 5E, F ). In addition, the 4T1-bearing mice had significant >2-fold increases in most counts of lymphoid populations compared to the other groups ( Figures 5C, E, F ). In both tumor models, NK cell subpopulation counts were significantly decreased at >2-fold levels compared to the no-tumor control animals, and 4T1 also showed a decrease in CD62L⁺ CD45RB⁺ CD8⁺ T cells compared to the other groups ( Figures 5D–F ). This illustrates the wide scope of the blood leukocyte count changes and their complexity in relation to the cancer growth, which also highlights the utility of adding markers to identify the leukocyte population changes based on a screening approach.

While fold changes and p values can help test the group similarities in leukocyte counts, we sought to generate a statistic that reflects the populations that are most capable of distinguishing between the tumor-bearing and healthy groups. This would be particularly useful in identifying the blood cell biomarkers that correlate with the presence and type of cancer in a clinical setting. For this, we assessed pairwise the measures of differences in the blood leukocyte population between CT26- and 4T1-bearing and no-tumor groups, including the metrics based on standard fold changes ( Figure 6A , |LgFC| column) and p-values ( Figure 6A , -Lg (p) column), but also an overlap statistic (41) to measure the overlap of the blood count distributions of a given leukocyte population between groups ( Figure 6A , OV & 1-OV columns). The rationale here is that the distribution of a particular leukocyte count that showed the least overlap between groups would be better at distinguishing the groups. Leukocyte counts with the largest fold change and smallest overlap between groups would be an ideal biomarker for distinguishing cancers and could be presented as a single statistic, with the ratio of fold change divided by overlap ( Figure 6A , |LgFC|/OV column). From this analysis, the top distinguishing population counts appeared to be from the monocyte subsets that had large |LgFC|/OV values ( Figure 6A ). Among these monocyte populations, the counts of monocytes that expressed higher levels of CD62L (Mo F+L+; monocytes expressing SiglecF and CD62L) had large fold changes and small p-values (larger -log2(p)) between the groups, and notably had minimal count overlap between the tumor-bearing groups and no-tumor controls, thus representing a suitable marker to distinguish cancer presence. In addition, monocytes that downregulate CD62L (Mo F+L- monocytes expressing SiglecF and lower levels of CD62L) also displayed similar statistical trends, but these were for distinguishing 4T1-bearing mice from both CT26-bearing and no-tumor control mice ( Figure 6A ). These two monocyte populations, therefore, harbor the potential to distinguish all groups. Indeed, plotting the counts of these populations ( Figures 6B, C ) showed a clear pathway to distinguish no-tumor animals from cancer-bearing animals and also differentiate between the two cancer types based on a decision tree analysis ( Figure 6D ). To highlight this further, these two monocyte populations were used alone or together to train and test a decision tree ML model to assess the capacity to predict the presence and type of cancer from blood leukocyte profiling ( Figures 6E–G ). The counts of monocytes expressing higher levels of CD62L (Mo F+L+) could classify the presence or absence of tumors but could not distinguish the tumor type ( Figure 6E ). The counts of monocytes expressing lower levels of CD62L (Mo F+L-) could distinguish 4T1 tumor from the other groups ( Figure 6F ). Using both monocyte populations, all groups could be distinguished ( Figure 6G ), highlighting the utility of overlap statistics in rapid identification of the group-delineating populations.

Overlap statistics measure differences in distributions pairwise, and therefore would be best suited to identify key distinguishing features between a few groups. In cases where larger number of groups are being defined or they have larger distribution overlap, more complex models might be required. To establish a generic approach for detecting leukocyte population counts that are most useful for group classification, we used a decision tree-based feature ranking built in the Random Forest and CATBoost tree-based ML models ( Figure 6A ML rank column). While the monocyte population remained one of the highest ranked features according to the model ranking, several other leukocyte populations of significance were also identified. A cumulative ranking profile based on Random Forest rank, CATBoost rank, p values and the fold-change over overlap ratio was generated to illustrate the overall feature importance across the measures ( Figure 6A right bar plot annotation). To narrow down the most important features, we created a learning curve to score the classification performance as a function of diminishing leukocyte groups added to the model, and removed the lowest CATBoost ranked populations first ( Figure 7A ). As can be seen, the model performed above 95% of all performance measures with 3 or more top-ranked features. However, there appeared to be variability in performance depending on the number of features, wherein the performance dropped with a large number and enhanced from 8 to 6 top-ranked populations, suggesting that the top 6 features ( Figure 7B ) might be the most informative. Of note, these 6 populations did not include monocyte (Mo F+L-) or neutrophil (N B-6a+; neutrophils expressing low levels CD45RB and higher levels of CD66a) populations that scored the highest on our cumulative ranking profile ( Figure 6A right bar plot annotation). Given that these two populations are potentially key leukocyte biomarkers for group segregation, we added them to the top 6 CATBoost ranked populations, making a total of 8 populations consisting of various monocyte, neutrophil, CD8⁺ T cell and B cell subsets ( Figure 7B ). To better understand how these 8 populations may perform in cancer classification in one another’s presence, they were compared against all populations for their segregating capacity by multidimensional scaling (MDS) analysis aimed at summarizing the metric distances across all of the 8 leukocyte subsets in a 2-dimensional space ( Figure 7C ). This demonstrated that the selected 8 populations improved group separation relative to all populations, particularly of the CT26 from the no-tumor control (Nil) groups. Indeed, the final CATBoost model trained with the blood counts of these 8 leukocyte populations yielded definitive prediction of all 3 groups ( Figure 7D ). A summary of how the counts of each population change across the groups and act as a collective biomarker profile is displayed as their log2 fold change relative to one another ( Figure 7E ). This clearly illustrates the differential changes in the myeloid (annotated Mye) and lymphoid (annotated Lym) populations in different types of cancer. Thus, our screening pipeline approach has allowed for identification and selection of cancer-specific leukocyte biomarker signatures that are predictive of the cancer presence and its type from a single blood sample.