For targeted surveillance of rare pathogens, screenings must be performed on a large 2 number of individuals from the host population to obtain a representative sample. For 3 pathogens present at low carriage rates of 1% or less, a typical detection scenario 4 involves testing hundreds to thousands of samples before a single positive is identified. 5 Although advances in molecular biology and genomic testing techniques have greatly 6 lowered the cost of testing, the large number of negative results still renders any of specimens in order to detect just a few thousand cases. The large number of negative 16 tests struck Dorfman as being extremely wasteful and expensive and he proposed that 17 more information could be gained per test if many samples were pooled together and 18 tested as a group [5] . If the test performed on the pooled samples was negative (which 19 was very likely), then all individuals in the group could be cleared using a single test. If 20 the pooled sample was positive, it would mean that at least one individual in the 21 sample was positive and further testing could be performed to isolate the positive 22 samples. This procedure had the potential to dramatically reduce the number of tests 23 required to accurately screen a large population and it sparked an entirely new field of 24 applied mathematics called group testing. 25 Due to practical concerns, Dorfman's group testing approach was never applied to 26 syphilis screening because the large number of negative samples had a tendency to 27 dilute the antigen in positive samples below the level of detection [6] . Despite this, 28 sample pooling has proven to be highly effective when using a sufficiently sensitive, 29 often PCR-based, diagnostic assay. In fact, ad hoc pooling strategies have long been 30 used to mitigate the costs of pathogen detection in disease surveillance programs. For 31 example, surveillance of mosquito vector populations in the U.S. involves combining 32 multiple mosquitoes of the same species (typically 1 -50) into a single pool, prior to 33 testing for the presence of viral pathogens [7] [8] [9] [10] . Elsewhere, such pooling techniques 34 have been successful in reducing the total number of tests in systems ranging from 35 birds [11] , to cows [12] , to humans [13] [14] [15] . In many wildlife/livestock surveillance 36 programs, sample pooling is used to simply determine a collective positive or negative 37 status of a population (e.g. a herd or flock) without identifying individual positive 38 samples. While this is often appropriate and sufficient for small-to-medium scale 39 research experiments or surveillance programs, a well designed pooling scheme can 40 easily provide this valuable information with little additional cost. For the purposes of 41 this paper, we will focus on pooling methods that provide accurate classification of each 42 sample so that infected individuals can be identified. 43 Group testing theory primarily focuses on minimizing the number of tests required 44 to identify positive samples and many nearly-optimal strategies for sample pooling have 45 been described. From a combinatorial perspective, a testing scheme begins by 46 examining a sample space which includes all possible arrangements of exactly k positive 47 samples in N total samples. Because the positive samples are indistinguishable from 48 negative samples, a test must be performed on a sample or a group of samples in order 49 to determine their status. The test is typically assumed to always be accurate, even 50 when many samples are tested together (in practice, this is often not the case and 51 approaches that consider test error and constraints on the number of samples per pool 52 have been examined [16, 17] ). In the worst case, all of the samples would need to be 53 tested individually requiring N tests. The goal of group testing is to devise a strategy 54 which tests groups of samples together in order identify the positive samples in fewer 55 than N tests. Group testing methods are generally more efficient when positive samples 56 are sparse. As the number of positive samples increases, the number of tests will 57 eventually exceed individual testing for all of the methods. This point has been 58 previously estimated to be roughly when the number of positives is greater than N 3 for 59 sufficiently large N [18, 19] . In order to establish the most optimal testing procedure, non-adaptive and adaptive pooling approaches and, in each case, we assume that the 106 test applied to the pools is noiseless (the test will always be positive if a positive sample 107 is present in the pool and negative otherwise) and it produces only a binary or two-state 108 outcome (e.g. positive/negative or biallelic SNP typing).
number of times any two samples are included in the same pool [20] . This is achieved by 117 staggering the samples that are added to each pool in different sized windows or 118 intervals ( Fig 1) ; importantly, the size of the windows must be greater than √ N and 119 co-prime to minimize the intersections between samples. The number of different 120 pooling windows (the weight) should be one greater than the expected number (upper 121 bound) of positive samples, w =k + 1, to ensure accurate results.
122 Fig 1. DNA Sudoku pooling example. In this example, there are a total of N = 96 samples. The 96-well plates show which samples are combined into each pool for the two different window sizes (W 1 = 10 and W 2 = 11 which are greater than √ N and co-prime). By using two different window sizes, the weight of this pooling design is w = 2 meaning that k = w − 1 = 1 positive sample can be unambiguously identified in a single step using T = W 1 + W 2 = 21 tests. The positive samples are decoded by finding the samples that appear most often in the positive pools. For example, if G10 is the only positive sample, we can detect this from the pooling results by noticing that G10 was added to both of the positive (red) pools while other samples in those pools were added to only one or the other. Alternatively, if both G10 and D4 are positive, four samples occur with equal frequency (D4, G10, E12, and F2) in the positive pools (red and purple) and it is impossible to determine which are the true positive samples. This ambiguity is introduced because the test was designed to handle only one positive sample. Multidimensional pooling is another non-adaptive approach that is generally easier to 142 perform than DNA Sudoku but can be more prone to producing ambiguous results. As 143 the name implies, this procedure can be extended to many dimensions [21, 22] , however 144 it becomes more difficult to perform without robotics when more than two dimensions 145 are used. In the two dimensional (2D) case, N samples are arranged in a perfectly 146 square 2D grid or in several smaller but still square sub-grids [23] . For example, when 147 testing 96 samples (as in Fig 2) , this could be achieved through a single 10x10 grid or Dorfman's original pooling design for syphilis screening was an adaptive two-stage test. 167 Following this method, samples are partitioned and tested in g groups of size n. All of 168 the samples in groups with negative results are considered to be negative and all of the 169 samples in groups with positive results are tested individually. Ignoring the constraints 170 of the actual assay, the optimal group size that minimizes the number of tests depends 171 on the number of positive samples, k. Specifically, there should be roughly √ N k groups 172 of size N k [5, 24] . Dorfman's two-stage approach was later generalized to any number 173 of stages using Li's S-Stage algorithm [24] , which can reduce the number of tests 
Sobel and Groll [25, 26] , introduced several adaptive group testing algorithms based on 189 recursively splitting samples into groups and maximizing the information from each test 190 result. They demonstrated that this class of algorithm is robust to inaccurate estimates 191 of k, particularly in the case of the Binary Splitting by Halving algorithm which can be 192 performed without any knowledge of the number of positive samples. Binary Splitting 193 by Halving (Fig 4) begins by testing all of the samples in a single pool. If the test is 194 negative, all of the samples are negative and testing is complete, if the test is positive, 195 the samples are split into two roughly equal groups and only one of the groups is tested 196 in each step. If the tested half is negative, we know that all of the samples in the tested 197 group are negative and testing is now complete for those samples. We also know that testing [25, 26] . This is the only approach discussed here that does not rely on an Step 2). If the tested half is negative, then all of the samples in the tested half are considered to be negative and at least one negative sample is known to be present in the other non-tested half of the samples. If the tested half is positive, then it contains at least one positive sample and no information is gained about the other untested half. In either case, the method continues by halving and testing whichever group is known to contain a positive sample until a single positive sample is identified (either by individual testing, as seen in Step 7, or by elimination, as seen in Step 16) . Once a single positive sample is identified, the remaining unresolved samples (non-grey wells) are pooled and tested to determine if any positive samples remain and the process continues until all positive samples are identified. Only one test is required per round, and in this example, it takes 17 sequential rounds to recover both positive samples. fewer tests [27] . As the ratio of samples to positive samples ( N k ) increases, the number 228 of tests required to identify k positive samples approaches k log 2 n k which is nearly 229 optimal; however, like Binary Splitting by Halving, the Generalized Binary Splitting 230 approach requires many sequential steps to complete testing. Here we are introducing a new approach that we developed with the goal of finding a 233 good balance between the number of tests, the number of steps, simplicity, and 234 robustness. We found that many of the methods described previously focus on 235 optimizing only one of these features usually to the detriment of the others. Instead of 236 attempting to perform the best in a single area, we wanted to take a more balanced 237 approach and find tradeoffs that allow good performance across each of these areas. Our 238 Modified 3-Stage approach (Fig 5) is based on the S-Stage approach but it is modified 239 so that the number of steps is constrained to a maximum of three. At three steps, this 240 approach requires only one additional step than ambiguous non-adaptive approaches 241 that require two steps for complete validation. Because the S-Stage algorithm is already 242 fairly robust, constraining the number of steps does not have a large impact on the 243 number of tests required. We also modified our method to be simpler and easier to 244 perform by borrowing the recursive subdividing used in the binary splitting approaches. 245 In the S-Stage approach, the remaining samples in each step are arbitrarily redivided 246 into pools. Not only does this make it difficult to keep track of the remaining samples 247 spread across the plate, it can also make it more difficult to collect the samples for a 248 pool using a multichannel pipette (e.g. Step 2 in Fig 3) . Instead, we opted to recursively 249 subdivide the samples from positive pools. This makes it easier to keep track of the 250 samples that should be pooled at each stage and, because the samples are always in close 251 proximity, they are easier to collect using a multichannel pipette (compare 3 and 5). 
is the number of groups tested at each step. The number of pipettings 300 for a single channel pipette was equal to the number of samples in each of the pools that 301 were tested. For multichannel pipettes, the number of samples in each pool was divided 302 by the number of channels and rounded up. In cases where the samples in the pool were 303 not in adjacent wells, additional pipettings were required. Experimental validation of modified 3-stage approach 336 We set up rare pathogen detection experiments in complex microbiome backgrounds to 337 test our Modified 3-Stage approach. We used a total of 768 samples (eight 96-well 338 plates) that contained a background of 2 µL of DNA extraction from cow's milk and 8 339 µL of molecular grade water. These samples originated from 24 distinct cow milk 340 samples and were replicated (32 replicates each) to fill eight 96-well plates -a total of 341 24 unique microbiome backgrounds. C. burnetti DNA (1 µL) was added to 10 randomly 342 chosen background samples (∼ 1.3% carriage rate) as we verified that the spike-in was 343 successful using a highly sensitive Taqman assay designed to target the IS1111 344 repetitive element in Coxiella burnetti [28] . Using the same Taqman assay, we also 345 verified that the target pathogen was not present in any of the 24 unique microbiome 346 backgrounds prior to the spike-in. To ensure a consistent amount of background DNA, 347 the milk extractions were tested to determine the amount of bacteria with a real-time 348 PCR assay that detects the 16S gene and compares it to a known standard [29] .
The pooling procedure was carried out by a typical researcher looking to identify The number of sequential steps is one of the major factors that differentiates pooling 397 methods. The major benefit of non-adaptive pooling methods is that, in some cases, all 398 of the tests can be run at the same time which means that testing can be completed 399 faster. Clearly, the nonadaptive tests required the fewest number of steps even when the 400 results were ambiguous, necessitating a second round of validation 7. For 96 samples, 401 the highest weight that we tested was 6 which meant that any simulation with 5 or The number of samples that are combined in a single pool is a very important practical 443 concern because it can determine whether the assay can produce accurate results. In most methods, using a multichannel pipette reduced the number of pipettings by 477 an order of magnitude in some cases. Compared to the 8-channel pipette, the likely because the method requires many samples to be pooled at each step for many 489 steps. The performance is slightly improved when multichannel pipettes are used but it 490 is still the least efficient in many cases. Using a single pipette, DNA Sudoku was not 491 the most inefficient compared to the other methods. However, because the samples that 492 are combined in each pool are spaced out in different intervals instead of in consecutive 493 groups, the number of pipettings did not improve by using multichannel pipettes. This 494 means that, in the best case, a laboratory technician would need to correctly pipette Table 1 which shows that method is more sensitive 513 to overestimation than to underestimation. Of the methods that depended on an 514 estimate of the number of positive samples, the S-Stage (Fig 10, left) and our Modified 515 3-Stage approach (Fig 10, second from left) were the most robust to misestimations of k. 516 The number of steps was more robust in the Modified 3-Stage approach than the 517 S-Stage due to the 3-Step constraint; however, the Modified 3-Stage was more sensitive 518 in the number of tests in some cases (Table 1) . DNA Sudoku (Fig 11, left) was the most sensitive method overall. Overestimating of 520 the number of positive samples caused the weight of the pooling design (w =k + 1) to 521 be set higher than it needed to be. When this happened, all of the positive samples were still unambiguously identified but each unnecessary increase in the weight required 523 more than √ N additional tests. When the number of positive samples was 524 underestimated, fewer tests were performed but the pooling scheme was no longer able 525 to unambiguously identify the positive samples in a single step and a second round of 526 verification was required. A similar pattern occurred in the 2D Pooling simulations (Fig 527  11, right) . While the grid dimensions did not directly depend on k, generally larger 528 grids were more efficient when the number of positive samples was low and smaller grids 529 reduced ambiguous results when the number of positive samples was high but at the 530 cost of many more tests. However, because 2D Pooling was constrained to two 531 dimensions, the number of tests did not vary as drastically as DNA Sudoku. Table 2 . Although the expected number of positive 541 samples per plate was ∼ 2 given the 1.3% carriage rate, the actual number of positives 542 ranged from 0 to 3 and none of the plates had exactly 2 positive samples ( Table 2 ). The 543 TaqMan assay was able to accurately identify the positive pools without any false 544 positives or false negatives even during the first step when the number of samples per 545 pool was the largest at 16. Using an 8-channel pipette where appropriate, a total of 180 546 pipettings was required to pool the samples. A total of 120 TaqMan assays were 547 performed which is ∼ 84% fewer than would be required to individually test 768 548 samples.
Picking the right pooling approach for a given pathogen surveillance campaign can be a 551 complicated decision, which is often driven by a set of conflicting constraints and complexity. DNA Sudoku, however, is far from optimal for monitoring rapidly changing 566 pandemics due to its extreme sensitivity to misestimation of the carriage rate of the 567 pathogen in population.
A good middle ground between the adaptive and non-adaptive pooling approaches is 569 the Modified 3-Stage approach -our preference in our own surveillance applications.
While it is never the absolute best in any one category, it is always nearly optimal in 571 terms of number of serial steps (2nd best), complexity (2nd best), number of tests (4th 572 best), and extremely resilient to misestimation of the carriage rate (2nd best). The 573 latter is particularly important, as it allows this approach to be useful for surveillance in 574 situations with rapidly changing pathogen carriage rates (e.g. in pandemic or seasonal 575 outbreaks), while keeping number of serial steps as low as possible for an adaptive 576 method.