Rafael A Lopes
email: rafael.lopes@gmail.com

Samuel Thibault
email: samuel.thibault@labri.fr

Alba C M A Melo

MASA-StarPU: Parallel Sequence Comparison with Multiple Scheduling Policies and Pruning

Keywords: Parallel sequence comparison, parallel programming environment, dynamic programming

Sequence comparison tools based on the Smith-Waterman (SW) algorithm provide the optimal result but have high execution times when the sequences compared are long, since a huge dynamic programming (DP) matrix is computed. Block pruning is an optimization that does not compute some parts of the DP matrix and can reduce considerably the execution time when the sequences compared are similar. However, block pruning's resulting task graph is dynamic and irregular. Since different pruning scenarios lead to different pruning shapes, we advocate that no single scheduling policy will behave the best for all scenarios. This paper proposes MASA-StarPU, a sequence aligner that integrates the domain specific framework MASA to the generic programming environment StarPU, creating a tool which has the benefits of StarPU (i.e., multiple task scheduling policies) and MASA (i.e., fast sequence alignment). MASA-StarPU was executed in two different multicore platforms and the results show that a bad choice of the scheduling policy may have a great impact on the performance. For instance, using 24 cores, the 5M x 5M comparison took 1484s with the dmdas policy whereas the same comparison took 3601s with lws. We also show that no scheduling policy behaves the best for all scenarios.

I. INTRODUCTION

In the last decades, we have observed an astonishing evolution of the sequencing methods, which allowed the rapid assembly of genetic sequences by a huge number of laboratories around the world. Even though this is clearly a great benefit, it created the so-called data deluge and, thus, genetic sequences are being produced in a rate that is much higher than the rate of their analysis [START_REF] Goodwin | Coming of age: ten years of next generation sequencing technologies[END_REF].

One of the first steps in biological sequence analysis is pairwise sequence comparison, where a newly obtained sequence is compared to sequences which have been catalogued, in search of similarities. Smith-Waterman (SW) [START_REF] Smith | Identification of common molecular subsequences[END_REF] is a wellknown algorithm that provides the optimal result. It uses dynamic programming and has quadratic time and space complexities. This leads to high execution times when the sequences compared are long (Megabase comparison).

In order to accelerate sequence comparison algorithms, parallel platforms composed of multicores or multicores and accelerators such as GPUs (Graphics Processing Units), Intel Xeon Phis and FPGAs (Field Programmable Gate Arrays) have been used. CUDAlign 4.0 [START_REF] Sandes | CUDAlign 4.0: Incremental speculative traceback for exact chromosome-wide alignment in GPU clusters[END_REF] is a tool that uses a variant of the SW algorithm to compare huge DNA sequences in GPUs. Its code has been re-structured into the MASA [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF] architecture and now it runs in multicores, GPUs and Intel Xeon Phis with various programming environments (CUDA, OpenCL, OmpSs and OpenMP). In addition, MASA incorporates the block pruning capability, which accelerates considerably the execution when the sequences compared have a high degree of similarity. One drawback of block pruning is that the pruning behaviour is determined during execution and, for this reason, scheduling issues may occur [START_REF] Figueiredo | Parallel comparison of huge dna sequences in multiple gpus with block pruning[END_REF].

StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] is a general-purpose task-based parallel programming environment which provides multiple task scheduling policies and runs in several parallel platforms. In StarPU, the programmer composes a graph of task dependencies and StarPU keeps track of the tasks that become ready. Then, the ready tasks are executed according to the selected task scheduling policy. Currently, StarPU offers more than 10 different scheduling policies [START_REF]StarPU handbook for StarPU 1.3.3[END_REF].

Related work in the area of Megabase DNA sequence comparison show that impressive performance can be attained with accelerators [START_REF] Sandes | CUDAlign 4.0: Incremental speculative traceback for exact chromosome-wide alignment in GPU clusters[END_REF] [8] [START_REF] Liu | SWAPHI-LS: Smith-Waterman algorithm on Xeon Phi coprocessors for long DNA sequences[END_REF]. The performance results for multicores (CPUs) are less impressive but this is the platform used at most Bioinformatics laboratories in developing countries and, for this reason, this is the target platform used in this paper. In the literature, very few CPU tools based on Smith-Waterman are able to align Megabase sequences. The popular Water tool (www.ebi.ac.uk/Tools/psa/emboss water) restricts the sizes of the sequences to less than 1M. MASA-OpenMP and MASA-OmpSs are part of the MASA framework [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF] and are able to align Megabase sequences, achieving the best performance of about 4.80 and 5.88 GCUPS (Billions of Cells Updated per Second), respectively, in a multicore platform with 12 cores when comparing two Megabase DNA sequences. ASW [START_REF] Tang | ASW: Accelerating Smith-Waterman algorithm on coupled CPU-GPU architecture[END_REF] is another tool able to align Megabase sequences, which obtained 7.2 GCUPS in a platform containing an APU (Accelerator Processing Unit) composed of 512 GPU cores and 4 CPU cores. In the CPU-only execution, ASW attained 0.46 GCUPS [START_REF] Tang | ASW: Accelerating Smith-Waterman algorithm on coupled CPU-GPU architecture[END_REF]. As far as we know, there is no Megabase DNA sequence comparison tool in the literature that supports multiple allocation policies.

This paper proposes and evaluates MASA-StarPU, a bio-logical sequence aligner that runs on top of StarPU and takes advantage of its multiple scheduling policies. The contributions of the paper are twofold. First, we integrate a generic parallel programming environment (StarPU) to a domainspecific programming environment (MASA), creating a tool which has the benefits of StarPU (i.e., multiple task scheduling policies) and MASA (i.e., fast sequence alignment with the SW algorithm). Second, we assess the behaviour of multiple task scheduling policies for sequence comparisons under block pruning and draw some directions about how to choose the most appropriate policy for a particular alignment scenario. MASA-StarPU was evaluated in two different multicore platforms: (a) a notebook Acer with 4 cores; and (b) a node (2 Intel Xeon E5-2680 -24 cores) in the Miriel machine from the PlaFRIM platform (www.plafrim.fr/en/home). Real DNA sequences whose lengths ranged from 10 KBP (Thousands of Base Pairs) to 5 MBP (Millions of Base Pairs) and 7 scheduling policies from StarPU were used in the tests.

MASA-StarPU was able to attain 18.41 GCUPS (24 cores) when comparing 5 MBP x 5 MBP sequences with the dmdas (dequeue modeling data aware sorted) task allocation policy [START_REF]StarPU handbook for StarPU 1.3.3[END_REF], with a total execution time of 24 minutes. If we compare the same sequences but change the scheduling policy to lws (local work stealing), the GCUPS drops to 7.68, with an execution time of 59 minutes. We also show that there is no scheduling policy which provides the best execution time for all pruning scenarios and sequence lengths in both platforms.

The remainder of this paper is organized as follows. In section II, we describe DNA sequence comparison algorithms, their data dependencies and the block pruning strategy. MASA and StarPU are explained in Section III. The design of MASA-StarPU is presented in Section IV. Section V presents and discusses experimental results. Finally, Section VI concludes the paper.

II. DNA SEQUENCE COMPARISON

A DNA sequence is viewed as an ordered set of characters that belong to the alphabet Σ = {A, T, G, C}. Pairwise sequence comparison is a core operation in Bioinformatics that determines the similarity between two sequences by producing an alignment, which highlights regions of coincident and distinct characters. Since the goal is to maximize the similarity score, positive values are assigned to coincident characters (match) and negative values are assigned to different characters (mismatch) or to gaps included in one of the sequences.

An alignment may be (a) global, if it contains all characters from the sequences; or (b) local, if it only contains a subset of the characters from the sequences. The focus of this paper is on local alignments and Figure 1 illustrates a local alignment between two DNA sequences. The values assigned for matches, mismatches and gaps are +1, -1 and -2.

The Smith-Waterman (SW) algorithm [START_REF] Smith | Identification of common molecular subsequences[END_REF] computes optimal local alignments with Dynamic Programming (DP) in quadratic time and space (O(mn)), where m and n are the lengths of the sequences. It executes in two phases: (1) obtain the optimal score and (2) traceback. In Phase 1, a DP matrix H is calculated using the recurrence relation in Equation 1, where ma, mi and g are the values assigned for matches, mismatches and gaps, respectively, and S[i] is the i-th character of sequence S. When the whole matrix is computed, this phase ends by returning the optimal score (H[i, j] -similarity score) and its (i, j) position.

S 0 A C C - T G C C S 1 A C C T T G C C +1 +1 +1 -2 +1 +1 +1 +1 score = 5
H i,j = max          H i-1,j-1 + (if S 0 [i] = S 1 [j] then ma else mi) H i,j-1 +g H i-1,j +g 0 (1)
A C C T T G C C A T 0 0 0 0 0 0 0 0 0 0 0 A 0 1 0 0 0 0 0 0 0 1 0 C 0 0 2 1 0 0 0 1 1 0 0 C 0 0 1 3 1 0 0 1 2 0 0 T 0 0 0 1 4 2 0 0 0 1 1 G 0 0 0 0 2 3 3 1 0 0 0 C 0 0 1 1 0 1 2 4 2 0 0 C 0 0 1 2 0 0 0 3 5 3 1 G 0 0 0 0 1 0 1 1 3 4 2 A 0 1 0 0 0 0 0 0 1 4 3 G 0 0 0 0 0 0 1 0 0 2 3 Fig. 2. SW DP matrix. The aligment path is shown in bold.

Phase 2 (traceback) starts at position i, j and follows the traceback path, retrieving the alignment from the end to the beginning, until a zero-valued DP cell is attained. Figure 2 presents the DP matrix for the sequences shown in Figure 1.

SW computes alignments using the linear gap model (i.e. each gap has the same value). However, the affine gap model obtains more biologically relevant alignments and, thus, Gotoh [START_REF] Gotoh | An improved algorithm for matching biological sequences[END_REF] proposed an algorithm that favours sub-sequences of gaps to the detriment of isolated ones. In this algorithm, three DP matrices are computed, instead of one.

The quadratic space complexity imposes a considerable limitation on the lengths of the sequences. To overcome this, a divide and conquer algorithm is proposed [START_REF] Hirschberg | A linear space algorithm for computing maximal common subsequences[END_REF] that recursively computes DP cells which belong to the optimal alignment in linear space O(min{m, n}). Then, Myers-Miller (MM) [START_REF] Myers | Optimal alignments in linear space[END_REF] combined [START_REF] Hirschberg | A linear space algorithm for computing maximal common subsequences[END_REF] and [START_REF] Gotoh | An improved algorithm for matching biological sequences[END_REF] into an algorithm that retrieves alignments according to the affine gap model in linear space.

It must be noted that, in SW and its variants, each cell (i, j) depends on three previously calculated cells: (i -1, j -1), (i, j -1) and (i -1, j) (Equation 1). So, each antidiagonal of the DP matrix can be computed in parallel, using the wavefront method (Figure 3(a)), which requires a synchronization barrier at the end of the antidiagonal computation. This requirement may be relaxed if we use the task graph built from the data dependencies, leading to a generic processing (Figure 3(b)).

When comparing long sequences, it is usual to have several areas in the DP matrix with very low scores. Since the goal is to maximize the score, there are some areas which will never lead to a score higher than the maximum score calculated up to the moment. Therefore, the calculation of these areas can be skipped, accelerating the computation. This is the idea of the block pruning optimization [START_REF] Sandes | Formalization of block pruning: Reducing the number of cells computed in exact biological sequence comparison algorithms[END_REF]. If the sequences are similar, the current maximum score grows quickly, leading to a big pruning area. On the other hand, dissimilar sequences have small scores and, consequently, small pruning areas.

Block pruning was implemented MASA [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF] for two types of parallelization (Figure 3). It was observed that the generic way (Figure 3 (b)) provides a better pruning area, since the shape of the task graph may be adapted during execution to a square shape, that leads to better pruning results. Also, the values assigned to matches, mismatches and gaps have a great impact on pruning. More details about block pruning analysis can be found in [START_REF] Sandes | Formalization of block pruning: Reducing the number of cells computed in exact biological sequence comparison algorithms[END_REF].

Block pruning transforms a diamond-shaped task graph into an irregular graph and this is challenging for the task scheduling policies. Since the graph's shape depends on the pruning area, which is built at execution time, we claim that no scheduling policy behaves better in several pruning scenarios. As far as we know, there is no proposal in the literature that provides multiple allocation policies for SW executions with pruning and this is the challenge addressed in this paper.

III. MASA AND STARPU PROGRAMMING ENVIRONMENTS

One of MASA's great virtues is the block pruning strategy (Section II) which has a great impact on reducing the execution time. In its most recent version [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF], MASA does not use block pruning in runs with more than one device because this optimization prunes the DP matrix at runtime, introducing a scheduling issue.

MASA [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF] is freely available and it executes comparisons in 5 stages, where stage 1 implements Phase 1 of the Gotoh variant of the SW algorithm (Section II) and stages 2 to 5 implement Phase 2 (traceback) with a modified version of the Myers-Miller algorithm (Section II).

The MASA architecture [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF] has 6 modules (Figure 4). The data management module stores data necessary to the execution, such as the input sequences and user parameters. The statistics module collects information about the execution, including execution time, percentage of blocks pruned, etc.

The coordination of the execution is done by the stage management module. Stages 1 to 3 execute in the device (multicore, GPU or IntelPhi) and the stage management module divides the DP matrix into partitions, which are subdivided into blocks, and dispatches the computation of each partition to the module that contains the platform-specific code. MASA offers 2 types of block pruning in the Block Pruning module and the programmer can choose one of them. The parallelization strategy is also user-selectable. Finally, the DP matrix computation is part of the platform-specific code, which is entirely provided by the programmer. Currently, there are several MASA versions such as MASA-CUDAlign (GPU), MASA-OpenMP (multicores and Intel Phi) and MASA-OmpSs (multicores). Among the MASA versions, the one which uses the OmpSs [START_REF] Duran | Ompss: a proposal for programming heterogeneous multi-core architectures[END_REF] parallel programming environment implements the dataflow parallelization and has achieved the best performance for CPUs [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF].

In order to generate a MASA version, the programmer uses the MASA-API. In this API, the IAligner class is the interface between the platform-independent part of MASA and the MASA versions (platform dependent). A MASA version has an entry point, provided in the MASA::EntryPoint call, which is inherited from the platform-independent part of the MASA library.

Each MASA version follows a class hierarchy. The class Aligner implements the IAligner methods to do initializations, align a partition and finalize the computation. The class AbstractAligner is a child class of IAligner that initializes block pruning, receives the first row and column of the partition and outputs the last row and column computed, as well as the highest score. There are two classes for block pruning: GenericBP and DiagonalBP, which will be invoked depending on the programmer's choice. The class that actually computes the cells of the DP matrix is AbstractBlockAligner, a child class of AbstractAligner which will schedule the blocks for execution. This parent class has two child classes: AbstractDiagonalAligner, which will compute blocks by diagonal, and AbstractBlockAligner, which will compute blocks using the dataflow model. The AbstractBlockAligner has methods to ScheduleBlocks and AlignBlock. The whole class diagram of MASA can be found in [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF].

StarPU [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF] is a freely available parallel programming environment. A StarPU application is a set of tasks which dependencies are expressed in a task graph. StarPU schedules tasks to workers when they become ready, using one of its several scheduling policies [START_REF]StarPU handbook for StarPU 1.3.3[END_REF]. StarPU has two classes of scheduling policies: regular and performance model based ones. The performance model based polices use codelets representing platform-specific code to enhance scheduling decisions. A codelet is a piece of code that has a performance model associated to it and is used to estimate the task's duration with regressions generated with execution histories [START_REF]StarPU handbook for StarPU 1.3.3[END_REF]. In the following paragraphs, we will describe shortly the StarPU policies used in MASA-StarPU.

Five regular StarPU policies were integrated to MASA-StarPU. The eager policy uses a central queue and workers retrieve tasks from the unique queue in a self scheduling [START_REF] Yew | Processor self-scheduling for multiple nested parallel loops[END_REF] basis. The prio policy also uses a central queue but schedules tasks using priorities provided by the programmer. In StarPU 1.2.9 used for this article, the priority range is [-5, 5] and 5 is the highest priority. In StarPU 1.3 and beyond, the range is arbitrary. Random uses a local queue for each worker and distributes tasks randomly. There are 2 work stealing [START_REF] Blumofe | Scheduling multithreaded computations by work stealing[END_REF] based policies: (a) ws uses a local queue per worker and assigns tasks to the worker which releases the task -if a worker is idle, it steals tasks from the most loaded worker; and (b) lws behaves as ws but steals tasks from neighbour workers first, and then from workers running farther in the topology.

Two performance based policies were integrated to MASA-StarPU. Dmda (dequeue modeling data aware) is a HEFT (Heterogeneous Earliest Finish Time) [START_REF] Topcuoglu | Performance-effective and low-complexity task scheduling for heterogeneous computing[END_REF] based policy that aims to minimize the tasks termination time and schedules tasks when they become available and taking into account data transfer time. The policy dmdas (dequeue modeling data aware sorted) sorts tasks by priority and uses dmda to schedule tasks of the same priority. Unlike prio, there is no pre-defined range of priority values.

IV. DESIGN OF MASA-STARPU

The main goal of MASA-StarPU is to provide multiple task scheduling strategies which are appropriate to block pruning executions. In other words, we propose the integration of the frameworks MASA and StarPU, taking advantage of multiple scheduling policies in block pruning executions.

The scope defined for MASA-StarPU is the execution of MASA's stage 1. This stage has the greatest performance impact on sequence comparisons [START_REF] Liu | SWAPHI-LS: Smith-Waterman algorithm on Xeon Phi coprocessors for long DNA sequences[END_REF] [8] [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF], also allowing greater parallelization.

StarPU defines the application as a task graph and this model is adequate for the dataflow parallelization shown in Figure 4. In MASA-StarPU, the partition (DP matrix) is divided into blocks and, due to data dependencies (Equation 1), it is not possible to calculate a block without previously calculating its adjacent ones (up, left, diagonal). In order to express the dependencies, we used the system of TAGs provided by StarPU where each task has a TAG id and a pointer to TAGs from other tasks indicating dependencies.

Respecting the data dependencies, there may be a great number of blocks that are ready at a given time. Since some StarPU task allocation policies (e.g., prio and dmdas) use priorities, we set priorities to favour the square processing, since this way leads to a better pruning area [START_REF] Sandes | Formalization of block pruning: Reducing the number of cells computed in exact biological sequence comparison algorithms[END_REF]. Our version of StarPU allowed a small range of priorities for prio, so we adjusted the priorities to the [-5, 5] interval. For dmdas, we set the priorities for block (bx, by) as min(-bx, -by), using the same scheme of MASA-OmpSs [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF].

The general structure of MASA-StarPU is shown in Algorithm 1. Each MASA implementation has an entry point. In line 24, it receives args, which are the paths to the files that contain the input sequences, and the procedure which will execute the aligner (StarPUAligner). Procedure StarPUAligner (Algorithm 1 line 1) first initializes its structures and then does the computation. InitializeStructures() (line 2) is a MASA procedure which was modified to initialize StarPU structures and start its execution (Algorithm 2). Then the procedure scheduleBlocks is called (Algorithm 1 line 3) and the MASA-StarPU execution terminates in line 4.

In Algorithm 1 lines 5 to 22, the procedure scheduleBlocks is presented. The loops in lines 6 to 8 implement the square shaped execution [START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF]. Parameters width and height represent, respectively, the width and height of the DP matrix in number of cells, and maxT hreads represents the number of threads used in the execution.

Procedure alignBlock (lines 12, 15) is responsible for creating a task to align a block and uses i and j, which are the position of the block in the DP matrix. The details of alignBlock(i,j) are given in Algorithm 3. The call starpu_task_wait_for_all() (Algorithm 1 line 20) indicates that StarPU needs to wait for all tasks to finish and, finally, the call starpu_shutdown() finalizes StarPU and updates StarPU parameters, such as the performance model.

Algorithm 2 shows the pseudocode of procedure initializeStructures. Line 3 indicates which callback function StarPU should perform if a task is submitted to the codelet initializes in line 2. As the focus of this work is multicore execution, only the CPU function has been defined. StarPU uses the C language, so the parameter cpu_func from the codelet expects a C function. However, the framework MASA was written in C++, where a static method is considered a C function. Therefore, parameter cpu_func is defined as a callBack method, whose function is to call the MASA method processBlock (bx, by) which in turn is responsible for actually computing the DP matrix in a block. Algorithm 3 presents the procedure alignBlock, which is responsible for creating the parallel tasks. Variables i0, i1, j0 and j1 define the upper left (i0, j0) and the bottom right (i1, j1) corners of the block and are defined by StarPUAlignerParameters [START_REF] Augonnet | StarPU: a unified platform for task scheduling on heterogeneous multicore architectures[END_REF]. Parameters bx and by are the coordinates of the block in the matrix.

Line 2 defines the task's priority. If the block is not pruned, a task is created to perform its computation (lines 3 to 13). Initially, the codelet parameters are initialized (line 4). Then, the task is created and its parameters are updated (lines 5 to 10). In line 9, function TAG returns a hexadecimal string which represents the tag id in StarPU. Then the express_deps(bx, by) (line 11) function is called to define the dependencies. For a generic block, the dependencies are expressed as starpu tag declare deps(T AG(i, j), 2, T AG(i-1, j), T AG(i, j -1)). Note that the dependency at the diagonal (i-1, j -1) is covered by dependencies (i-1, j) and (i, j -1). Finally, the task is submitted for execution (line 12). Figure 6 presents a flowchart that summarizes Algorithms 1, 2 and 3, executed by MASA and StarPU in MASA-StarPU.

Figure 5 presents the placement of MASA-StarPU among the other MASA versions. MASA-StarPU was developed for CPUs, uses the dataflow parallelization and generic block pruning.

Algorithm 3 MASA-StarPU -alignBlock 1: procedure STARPUALIGNER::ALIGNBLOCK(bx,by) 2: priority = min(-i0, -j0); 3:

if not isBlockPruned(bx, by) then 4:

initializeParameters(bx,by,i0,j0,i1,j1,this) 5: task = starpu task create() 6:

task.cl = cl; 7:

task.cl arg = params 8:

task.use tag = 1 9:

task.tag id = TAG(bx, by) 10:

task.priority = priority 11:

express deps(bx,by) 12:

starpu task submit(task) 13:

end if 14: end procedure 5M. The optimal local scores obtained with MASA-StarPU are shown in Table I since they indicate the similarity between the sequences and the potential for block pruning. For instance, sequences 50K x 50K have low similarity whereas sequences 5M x 5M have high similarity.

The SW parameters used were: match: +1; mismatch -3; first gap: -5; extension gap: -2 and the size of the block was empirically set to 2048 2048. The seven policies explained in Section III were used in the tests.

First, we evaluate the impact of the scheduling policies, execution platform and block pruning rate in the execution time/GCUPS for the seven comparisons shown in Table I. GCUPS is computed by dividing the size of the DP matrix (m * n) by the execution time in seconds * 10 9 . Since the environments were dedicated to our experiments, the standard deviation in the measurements was very low.

Figure 7 shows the GCUPS for MASA-StarPU in the notebook (a) and PlaFRIM (b) platforms. Since MASA-StarPU has quadratic time complexity, we would expect that the GCUPS grow if the lengths of the sequences augment, until the maximum parallelism provided by the device is attained [START_REF] Sandes | Cudalign 3.0: Parallel biological sequence comparison in large gpu clusters[END_REF]. However, block pruning does not compute some blocks in the DP matrix and, when the sequences are very similar, more than 50% of the matrix may be pruned [START_REF] Sandes | Formalization of block pruning: Reducing the number of cells computed in exact biological sequence comparison algorithms[END_REF]. So, in Figure 7 (a) and (b) we can notice a small dip in the 3M region in the graph. This happens because the 3M sequences have very low similarity (Table I) whereas 1M sequences have medium similarity and 5M sequences have high similarity. So, the dips in both graphs show the block pruning effect.

Even though the shapes of the graphs in both platforms are similar, the behaviour of the policies is different. In Figure 7(a), we observe that, with the exception of the 5M comparison, the policies have a close GCUPS. This happens because the number of cores ( 4) is small. Even so, there is a significant difference among the scheduling policies in the 5M comparison. In this comparison, two factors were combined: a big DP matrix (5M x 5M) and a big pruning rate. This led to considerable differences: dmdas and prio have the highest GCUPS whereas ws, lws and random have the lowest ones.

In Figure 7(b), the scheduling policy has a bigger impact in the performance. For all comparisons, the difference in the GCUPS is higher than 55%, if we consider the lower and higher value for each execution. The 5M comparison presented the highest GCPUS variation (139.7%) and the best GCUPS (18.41) was obtained with the dmdas policy. Even though prio had a very good performance in the notebook (Figure 7 Table II shows the GCUPS and pruning ranges (lowesthighest) as well as the policy which obtained the best GCUPS for platforms Notebook and PlaFRIM. The gray scale in the rows indicate the pruning rate. Pruning rates lower than 10% are white and the higher the pruning rate, the darker the color. We can see that, with the exception of the 10K comparison, there is no great variation on the pruning rate when we change the computing platform (Notebook or PlaFRIM). A big variation occurred in the 10K comparison because the length of the sequences compared is close to the size of the block (2K). So, a slight variation on the number of blocks pruned led to a big percentage variation. In the comparison of longer sequences, this big variation in the pruning range between platforms did not occur.

Concerning the best policy, we can see that eager and random did not achieve the best result for any comparison in any platform. This was expected since these policies are very simple and the tasks which compose our application produce a complex task graph. When considering each platform alone (Notebook or PlaFRIM), we can see that there is no policy which provides the best result for all comparisons. For sequences smaller than 1M and with low pruning rate (50K, 150K and 500K), the policy ws had high GCUPS in the notebook. However, ws did not have good results in PlaFRIM -e.g, for the 150K comparison, the ws GCUPS was 3.86 and the best GCUPS was 4.93 (Table II).

In comparisons of sequences longer or equal to 1M with pruning rate higher than 10% (1M and 5M), the dmdas policy has a very good performance in both platforms. For long sequences with low pruning rate (3M), the best policies for smaller sequences (lws and ws, for the PlaFRIM and Notebook platforms, respectively) still have very good performance. When comparing 3M sequences in the notebook, the policy dmdas also provide very good results.

Figure 8 presents the pruning rates and GCUPS in two platforms for comparisons which have pruning rates higher than 1%. It can be seen that the scheduling policy has significant impact on the pruning rate. With the exception of Figures 8(a) and 8(d), the pruning results in both platforms are comparable. In the 10K comparison using the notebook (Figure 8(d)), the reduced number of cores had a negative impact on the pruning rate, which dropped to half the rate obtained in the PlaFRIM platform (Figure 8(a)). In both platforms, the best pruning rate for the 10K comparison was achieved by prio.

Although there is a variation on the pruning rate for the 1M comparison in both platforms (Figures 8(b) and 8(e)), this variation is small (less than 2%). For the 5M comparison, the best pruning rates were achieved by dmdas and prio. Not surprisingly, these are the policies which use priorities. Work stealing policies (ws and lws) did not provide good pruning results in either platform.

Looking at the results presented at this section, we can make some general observations. When the pruning rate is expected to be low, simple scheduling heuristics which do not take priorities into account will be fine. The overhead of taking priorities into account is not worth the little pruning optimization that could be obtained. Nevertheless, when the pruning rate is expected to be high, one really needs to use a scheduling strategy which takes priorities into account. Additionally, on platforms with a large number of cores such as PlaFRIM, the prio and eager policies, which use a single task queue, suffer from contention, so the dmdas policy should be used instead since it employs one task queue per core and ). The execution times measured here are the wallclock times. So, they include the initialization time (MASA, StarPU and OmpSs), as well as reading the sequences from files and writing the output.

MASA-StarPU attains the smallest execution time for all the comparisons (Table III). The speedup over the second best tool ranged from 1.4x to 1.8x. MASA-StarPU was able to achieve 18.41 GCUPS and, in the same comparison, MASA-OmpSs, which also follows the square shape, attains 12.72 GCUPS.

VI. CONCLUSION AND FUTURE WORK

In this work, we proposed and evaluated MASA-StarPU, a solution that uses parallel computing to compare long DNA sequences with support for multiple task scheduling strategies. StarPU was integrated to MASA, with the objective of efficiently computing the optimal score and its coordinates in the SW matrix. The main features of MASA, such as dataflow processing and block pruning (BP) have been incorporated into MASA-StarPU, contributing to very good performance.

MASA-StarPU was executed in two platforms, obtaining, in general, varied performances according to the task scheduling policy, the lengths of the sequences, the pruning rate and the number of cores. The impact of the these factors in the execution time/GCUPS for the 7 comparisons was evaluated. In the notebook, with the exception of the 5M comparison, the scheduling policies obtained a close GCUPS and dmdas and prio had the highest GCUPS whereas ws, lws and random had the lowest ones. In PlaFRIM, we observed that the scheduling policy has a bigger impact on performance, and the central queue design of prio and eager affected their performance compared to the distributed dmdas and dmda policies. For all comparisons, the difference in GCUPS was higher than 55%, if we consider the lower and higher value. The 5M comparison presented the highest GCPUS variation (139.7%) and the best GCUPS (18.41) was obtained with the dmdas policy.

As future work, we intend to retrieve also the optimal alignment between the two sequences. To achieve this objective, the other MASA stages should be incorporated to MASA-StarPU. We also intend to analyze the effect of the block size on the behavior of the scheduling policies.

In this paper, we observed that it was not possible to determine the best scheduling policy for all cases. However, we think that it is possible to evaluate specific cases and then determine the best policy for such cases. To achieve this, we plan to create an intelligent mechanism that tries to learn the correlation among platform, scheduling policy and pruning rate and acts accordingly.

Fig. 1 .

 1 Fig. 1. Example of local alignment between two DNA sequences.

Fig. 3 .

 3 Fig.3. Antidiagonal and generic way to calculate the DP matrix[START_REF] Sandes | MASA: A multiplatform architecture for sequence aligners with block pruning[END_REF].

Fig. 4 .

 4 Fig. 4. The MASA Architecture [4].

Algorithm 2

 2 MASA-StarPU -initializeStructures 1: procedure INITIALIZESTRUCTURES 2: starpu codelet init(cl) 3: cl.cpu func := callBack 4: cl.model := masa perf model 5: masa perf model.type := STARPU HISTORY BASED 6: masa perf model.symbol := "masa perf model" 7: starpu init(NULL); 8: end procedure Line 4 initializes the parameter model and line 5 indicates that the model should be updated using the execution history. If this file does not exist, it is created after the first call to starpu_shutdown(), with the information from the last execution.

Fig. 5 .Fig. 6 .

 56 Fig. 5. MASA-StarPU extension in the MASA architecture

Fig. 7 .

 7 Fig. 7. GCUPS for several sequence comparisons. Higher is better. it had an average performance in PlaFRIM (Figure 7(b)). As expected, random had bad results in all comparisons.TableIIshows the GCUPS and pruning ranges (lowesthighest) as well as the policy which obtained the best GCUPS for platforms Notebook and PlaFRIM. The gray scale in the rows indicate the pruning rate. Pruning rates lower than 10% are white and the higher the pruning rate, the darker the color. We can see that, with the exception of the 10K comparison, there is no great variation on the pruning rate when we change the computing platform (Notebook or PlaFRIM). A big variation occurred in the 10K comparison because the length of the sequences compared is close to the size of the block (2K). So, a slight variation on the number of blocks pruned led to a big percentage variation. In the comparison of longer sequences, this big variation in the pruning range between platforms did not occur.Concerning the best policy, we can see that eager and random did not achieve the best result for any comparison in any platform. This was expected since these policies are very simple and the tasks which compose our application produce a complex task graph. When considering each platform alone (Notebook or PlaFRIM), we can see that there is no policy which provides the best result for all comparisons.

Fig. 8 .

 8 Fig. 8. Pruning rate for the 10Kx10K, 1Mx1M and 5Mx5M comparisons. GCUPS are shown in the bottom between brackets.

  , with Linux Ubuntu 16.04. The machine in (b) was chosen because many Bioinformatics laboratories in developing countries use notebooks and we want to evaluate MASA-StarPU in this adversary scenario.We compared real DNA sequences (TableI) retrieved from the National Center for Biotechnology Information (NCBI) at www.ncbi.nlm.nih.gov. The sequence sizes vary from 10K to

	V. EXPERIMENTAL RESULTS
	MASA-StarPU	was	implemented	in	C/C++
	using MASA-core release 1.3.9.1024 (available at
	github.com/edanssandes/MASA-Core) and StarPU release
	1.2.9 (available at gforge.inria.fr/projects/starpu/ ).
	Two execution platforms were used: (a) PLaFRIM (Feder-
	ative Platform for Research in Computer Science and Math-
	ematics), a platform hosted at INRIA/France, composed of
	several parallel computing environments -in this paper, we
	used one multicore (24 cores) at the Miriel machine (2 x
	12 core Haswell Intel¨Xeon¨E5-2680v3, 2.5 GHz, 64 GB
	RAM, with Linux CentOS Release 7.1.1503); and (b) one
	notebook Acer (4 cores) with Intel i7-7700HQ, 2.8 GHz,
	16 GB RAM

TABLE I SEQUENCES

 I USED IN THE TESTS.

	Comparison	Sequence 1 Accession Size	Sequence 2 Accession Size	Similarity Score
	10K	AF133821.1	10K AY352275.1	10K	5091
	50K	NC 001715.1 57K AF494279.1	57K	52
	150K	NC 000898.1 162K NC 007605.1 172K	18
	500K	NC 003064.2 543K NC 000914.1 536K	48
	1M	CP000051.1	1M AE002160.2	1M	88353
	3M	BA000035.2	3M BX927147.1	3M	4226
	5M	AE016879.1	5M AE017225.1	5M	5220960

TABLE II GCUPS

 II , PRUNING RANGES AND BEST POLICY FOR 7 COMPARISONS IN THE NOTEBOOK AND PLAFRIM PLATFORMS

	Comp.		Notebook			PlaFRIM
			(4cores)			(24cores)
		Pruning	GCUPS	Best	Pruning	GCUPS	Best
		range	range	policy	range	range	policy
		(%)
	10K	13.8-20.8	0.7-1.1	prio	33.4-43.7	0.3-1.8	dmda
	50K	0-0	0.9-1.3	ws	0-0	1.8-4.5	dmda
	150K	0-0	1.2-1.4	ws	0-0	2.4-4.9	dmda
	500K	0-0	1.3-1.4	ws	0-0	3.3-6.2	dmda
				lws
	1M	10.6-12.2	1.5-1.6	dmdas	10.9-12.1	4.0-7.1	dmda
				prio
	3M	0.1-0.1	1.3-1.4	dmda	0.1-0.1	3.9-6.3	dmda
				ws
	5M	50.1-66.4	2.8-4.1	dmdas	43.7-66.0	7.6-18.4	dmdas

  Concerning the work-stealing-based lws and ws policies, they seem to fail at distributing over many cores the very irregular pattern of the task graph.Finally, we compared the performance of MASA-StarPU with the state-of-the-art (MASA-OpenMP and MASA-OmpSs) in PlaFRIM. TableIIIshows the execution times of the best MASA-StarPU policy and the execution times of these two tools, as well as the speedup obtained by MASA-StarPU. The speedup was calculated by dividing MASA-StarPU execution time by the second best execution time (bold in the figure

			TABLE III
	COMPARISON BETWEEN MASA-STARPU, MASA-OPENMP AND
			MASA-OMPSS (24 CORES)
	Comp.	MASA-	MASA-	MASA-	Speedup	GCUPS
		StarPU	OpenMP	OmpSs		MASA-
		(s)	(s)	(s)		StarPU
	10K	0.05	0.07	0.08	1.40x	1.86
	50K	0.71	1.37	1.28	1.80x	4.54
	150K	5.65	9.87	9.20	1.62x	4.93
	500K	46.80	76.73	77.17	1.63x	6.21
	1M	156.91	247.16	258.23	1.57x	7.14
	3M	1621.42	2472.80	2637.42	1.52x	6.37
	5M	1484.08	2824.90	2147.63	1.44x	18.41
	distributes work evenly thanks to the performance models.

ACKNOWLEDGMENT

The authors would like to thank the PLaFRIM administrators for providing access and support to use their platform.

This work is partially supported by project Capes/Procad 183794.


V Sanchis
email: vsanchis@pasteur.fr

Occurrence and Linkage Between Secreted Insecticidal Toxins in Natural Isolates of Bacillus thuringiensis Sylvain Espinasse, 1 Josette Chaufaux, 1 Christophe Buisson, 1 Ste ´phane Perchat, 1 Michel Gohar, 1,2 Denis Bourguet, 1 Vincent Sanchis 1,3 1 Unite ´de Recherches de Lutte Biologique, INRA La Minie `re, 78285 Guyancourt Cedex, France 2 Aventis Crop Science, Jozef Plateaustraat 22-B 9000 Ghent, Belgium 3 Unite ´de Biochimie Microbienne, Institut Pasteur, 25 rue du Dr. Roux, 75724 Paris Cedex 15, France Abstract. Little is known about the occurrence and linkage between secreted insecticidal virulence factors in natural populations of Bacillus thuringiensis (Bt). We carried out a survey of 392 Bt strains isolated from various samples originating from 31 countries. The toxicity profile of the culture supernatants of these strains was determined individually against Anthonomus grandis (Coleoptera) and Spodoptera littoralis (Lepidoptera). We analyzed ␤-exotoxin I production and searched for the genes encoding Vip1-2, Vip3, and Cry1I toxins in 125 of these strains. Our results showed that these insecticidal toxins were widespread in Bt but that their distribution was nonrandom, with significant linkage observed between vip3 and cry1I and between vip1-2 and ␤-exotoxin I. Strains producing significant amounts of ␤-exotoxin I were more frequently isolated from invertebrate samples than from dust, water, soil, or plant samples.

Although microbial populations display apparent overall genotypic and phenotypic homogeneity, many microbial species exhibit considerable diversity, which has an impact on the biological properties of individual strains [START_REF] Lan | Intraspecies variation in bacterial genomes: The need for a species genome concept[END_REF]. This is especially important for microorganisms with pathogenic properties, such as Bacillus thuringiensis (Bt). Bt, which has a worldwide distribution [START_REF] Martin | Worldwide abundance and distribution of Bacillus thuringiensis isolates[END_REF], is a sporogenic bacterium used in biological control that synthesizes specific insecticidal toxins [START_REF] Schnepf | Bacillus thuringiensis and its pesticidal crystal proteins[END_REF]. To date, most studies on Bt insecticidal toxins have focused on ␦-endotoxins, which are produced as parasporal inclusions, whereas little is known about the pathogenic potential of the compounds secreted into the culture medium [START_REF] Barjac | Isolation of an identifiable nucleotide from the thermostable toxin of Bacillus thuringiensis var Berliner[END_REF][START_REF] Kostichka | Cloning of a cryV-type insecticidal protein gene from Bacillus thuringiensis: The cryV-encoded protein is expressedearly in stationary phase[END_REF]. Many Bt strains secrete vegetative insecticidal proteins (Vips) during vegetative growth [START_REF] Warren | Vegetative insecticidal proteins: Novel pro-teins for control of corn pests[END_REF]. Two classes of Vip toxins have been described. The first consists of a binary system composed of two proteins, Vip1 and Vip2, which are 100 kDa and 52 kDa in size, respectively. These proteins are highly toxic to certain coleopteran species [START_REF] Warren | Vegetative insecticidal proteins: Novel pro-teins for control of corn pests[END_REF][START_REF] Warren | Auxiliary proteins for enhancing the insecticidal activity of pesticidal proteins[END_REF]. The second class consists of an 82.5-kDa pro-tein, Vip3, which is active against a wide spectrum of lepidopteran insects [START_REF] Estruch | Vip3A, a novel Bacillus thuringiensis vegetative insecticidal protein with a wide spectrum of activities against lepidopteran insects[END_REF]. These two classes of protein do not display sequence homology. Cry1I (also referred as to CryV in the literature) is another insecticidal toxin secreted during early stationary phase. It is active against certain lepidopteran and coleopteran insect larvae [START_REF] Kostichka | Cloning of a cryV-type insecticidal protein gene from Bacillus thuringiensis: The cryV-encoded protein is expressedearly in stationary phase[END_REF]. Finally, several Bt strains also produce a thermostable exotoxin, known as ␤-exotoxin I or thuringiensin [START_REF] Barjac | Isolation of an identifiable nucleotide from the thermostable toxin of Bacillus thuringiensis var Berliner[END_REF]. This compound is an adenine nucleotide analogue, highly toxic to a wide range of insect species [START_REF] Farkas | The structure of exotoxin of Bacillus thuringiensis var. gelechiae[END_REF], that is thought to inhibit RNA polymerase [START_REF] Sebesta | Mechanism of inhibition of DNA-dependent RNA polymerase by exotoxin of Bacillus thuringiensis[END_REF], thereby affecting insect molting and pupation, in some cases having teratogenic effects [START_REF] Burgerjon | Teratology of the Colorado potato beetle, Leptinotarsa decemlineata, as provoked by larval administration of the thermostable toxin of Bacillus thuringiensis[END_REF]. ␤-exotoxin I is also toxic to mammalian cells [START_REF] Beebee | Differential inhibition of mammalian ribonucleic acid polymerases by an exotoxin from Bacillus thuringiensis. The direct observation of nucleoplasmic ribonucleic acid polymerase activity in intact nuclei[END_REF][START_REF] Mackedonski | Inhibition of ribonucleic acid biosynthesis in mice liver by the exotoxin of Bacillus thuringiensis[END_REF] and is very persistent in the environment [START_REF] Benz | On the chemical nature of the heat stable exotoxin of Bacillus thuringiensis[END_REF]. It has therefore been banned for public use in accordance with WHO recommendations [START_REF] Who | Guidelines specification for bacterial larvicides for public health use[END_REF]. If we are to improve the efficacy and ensure safety of Bt products, we need to increase our understanding of the genetic variability of natural Bt populations with respect to the expression of these secreted virulence factors. In this study, we analyzed the toxicity of culture supernatants of natural Bt isolates to larvae of Spodoptera litto-ralis (Lepidoptera: Noctuidae) and Anthonomus grandis (Coleoptera: Curculionidae). We then explored the genomes of 125 natural Bt isolates for the genes encoding the four known insecticidal toxins secreted by Bt. The prevalence of these toxins in this species was determined and their distribution in relation to sample origin analyzed. We found that the secreted insecticidal toxins were not randomly distributed among Bt populations.

Materials and Methods

Bacterial cultures and supernatant preparation for toxicity assays.

The INRA collection contains 1260 natural Bt strains, representing isolates from 101 countries over the world, isolated from various sources (soil, insects, dust, plants, animal waste, etc.). These strains were primarily characterized as B. thuringiensis on the basis of their ability to produce crystals during sporogenesis [START_REF] Chaufaux | Research on natural strains of Bacillus thuringiensis in different biotopes throughout the world[END_REF]; they are conserved as spores in 12% vol/vol glycerol suspension at Ϫ20°C. 392 Bt strains, representing isolates from France and from 31 other countries across the five continents, were randomly taken from our collection. The strains were grown on LB nutrient agar plates at 30°C. Liquid cultures, in 100-mL LB, were inoculated from a single colony and shaken for 24 h at 30°C. The cultures were harvested at mid-sporulation (before cell lysis) by centrifugation at 14,000 g for 10 min at 4°C. The supernatants were filtered twice through 0.2-m Nalgene filter units (Nalgene). The resulting extracts, to be used for toxicity assays and ␤-exotoxin I determination, were stored at Ϫ20°C until bioassay.

Insect bioassays. The toxicity of the Bt supernatant preparations was determined using neonate larvae of S. littoralis and eggs of A. grandis, using a free ingestion technique. S. littoralis larvae were fed with an artificial diet [START_REF] Poitout | Elevage de plusieurs espe `ces de Le ´pidopte `res Noctuidae sur milieu artificiel riche et sur milieu artificiel simplifie ´[END_REF], dispensed into 50-well plastic plates (1.65-cm 2 surface each). Supernatant (25 L) was applied uniformly over the food surface and allowed to dry. One neonate larvae was placed on each of 35 wells and incubated for 10 days at 25°C, with a photoperiod of 16:8 (L:D) hours and 70% relative humidity (RH). Mortality was recorded at 3, 6, and 10 days. For A. grandis, the supernatant was incorporated into an artificial diet to a final concentration of 8% and dispensed into sterile 24-well plates. The diet consisted of a mixture of soya flour, cotton flour, wheat germ oil, cholesterol, sugar, vitamins, and agar. One egg of A. grandis was placed on the food in each well. Control larvae developed on untreated food for three weeks after hatching at 25°C. Percentages of mortality were calculated at three weeks after Abbott correction [START_REF] Abbott | A method for computing the effectiveness of an insecticide[END_REF].

PCR amplifications for cry1I, vip3, and vip1-2 detection. The oligonucleotides used in this study are listed in Table 1. The specific primers used to detect the cry1I, vip2, and vip3 genes were designed from the cry1I, vip2, and vip3 gene sequences present in the databases using Primer Select from the DNAStar software (DNAStar Inc). DNA was obtained from cells grown overnight on LB-agar plates and resuspended in 100 L of distilled water. Cell membranes were disrupted by freezing at Ϫ70°C and immersion in a 98°C water bath (heat shock). The cells were centrifuged and the DNA suspension (2 L) was mixed with 0.5 L of each primer (0.1 M), 0.5 U Taq polymerase, water (50 L), MgCl 2 , and buffer, according to the manufacturer's instructions (Gibco-BRL). PCR conditions were as follows: denaturation for 5 min at 94°C, followed by 35 cycles of 1 min at 94°C, 1 min at 50°C, 1 min and 30 s at 72°C, and a final elongation for 10 min at 72°C. Amplifications were recorded as positive, only when specific fragments of the expected size were found after migration on a 1% agarose gel. PCR results were considered to be negative if no amplicon was obtained after two attempts with independent DNA preparations.

Detection and quantification of ␤-exotoxin I. ␤-exotoxin I was extracted from the culture supernatant by solvent extraction and quantified by HPLC [START_REF] Gohar | Sample preparation for beta-Exotoxin determination in Bacillus thuringiensis cultures by reversed-phase high-performance liquid chromatography[END_REF]. Briefly, for solvent extraction, acetone was added to the exotoxin in 0.2 mL of the culture supernatant to a final concentration of 90%, and the mixture was centrifuged. The pellet was resuspended in 0.2 mL of double distilled (dd) water. Acetonitrile was added to a final concentration of 40% and the mixture was then centrifuged. The pellet was discarded and the acetonitrile concentration of the supernatant brought up to 90%. The precipitate was collected by centrifugation and the pellet was resuspended in 100 L of 50 mM potassium phosphate buffer, pH 2.5. For HPLC, we injected 25 L of the sample into a Lichrospher (Merck) C18 endcapped 4-ϫ 250-mm column. A gradient of 5% to 15% methanol in 50 mM potassium phosphate buffer, pH 2.5 was developed over 10 min. The flow rate was 1 mL/min, and UV absorption was monitored at 260 nm. ␤-exotoxin I was eluted at 5.5 min. The detection limit of this method for ␤-exotoxin I was 2 g/mL. A standard sample (70% purity) of ␤-exotoxin I was kindly provided by I. Thiery from the Laboratoire des Bacte ´ries Entomopathoge `nes (Institut Pasteur, Paris, France).

Statistical analyses. Frequencies of virulence factors.

Total frequencies of occurrence for each virulence factor were calculated from the experimental data. Assuming that there was no linkage between the virulence factors, the expected frequencies for each combination of the virulence factors is the product of the frequencies of the presence or of the absence of each virulence factor. For example, let us consider the four virulence factors {a, b, c, d} present with frequencies of respectively {f a , f b ,

f c , f d }. The expected frequency of the combination (a, b, c) is F abc ϭ f a ⅐ f b ⅐ f c ⅐ (1 Ϫ f d ).
Expected frequencies were compared to the observed ones, using the Pearson chi-square test, the null hypothesis being no linkage between the virulence factors. The Pearson correlation coefficient was calculated for each pair of virulence factors. The significance of all the correlations in the matrix was tested with the Bartlett chi-square test, and the probabilities associated with each correlation were corrected for multiple tests using the Bonferroni method [START_REF] Holm | A simple sequentially rejective multiple test procedure[END_REF]. Relationship between sample or geographic origin and genotypic or toxicity profiles. For each toxicity profile (i.e. nontoxic, toxic only to S. littoralis, or toxic only to A. grandis), the homogeneity of the frequencies observed for the different categories in sample or geographic origin were tested using a two-sided Fisher's exact test on contingency tables. Fisher's exact tests were performed using the STRUC program implemented in GENEPOP 3.2a [START_REF] Raymond | GENEPOP (version 1.2): Population genetics software for exact tests and ecumenicism[END_REF].

Results

Toxicity of stationary-phase culture supernatants of Bt isolates to S. littoralis and A. grandis. We assessed the activity of the culture supernatants of 392 Bt isolates against larvae of S. littoralis and A. grandis (Table 2). Bt isolates were differentiated into four groups: nontoxic, toxic only to S. littoralis, toxic only to A. grandis, or toxic to both insects. We found that 212 culture supernatants (54.1%) were inactive or only slightly toxic to both insect species (0 -30% mortality), 139 were toxic to S. littoralis only, 33 were toxic to both S. littoralis and A. grandis, and the supernatants of eight isolates were toxic only to A. grandis.

Screening of 125

Bt isolates for the presence of cry1I, vip3, and vip1-2 genes and for the production of ␤-exotoxin I. We selected 125 isolates for further analysis and investigated whether these toxicity patterns were correlated with the presence of one or several of the four known secreted insecticidal toxins. These strains were representative of the various toxicity profiles observed among the 392 samples initially studied and originated from different locations and sampling sources.

The number of samples with each of the four toxicity profiles was representative of the frequency of these profiles in the bioassay. For each strain, we determined the amount of ␤-exotoxin I present in the culture supernatant, using an improved quantification method based on reverse-phase HPLC [START_REF] Gohar | Sample preparation for beta-Exotoxin determination in Bacillus thuringiensis cultures by reversed-phase high-performance liquid chromatography[END_REF]. ␤-exotoxin I production (Ͼ2 g/mL) was detected in 47.2% of the strains. We also used PCR to investigate whether these 125 strains contained the cry1I, vip3, and vip1-2 genes. The PCR products were of the expected sizes, 378 bp for cry1I, 676 bp for vip3, and 914 bp for vip1-2 (Table 3). The cry1I gene was detected in 57.6%, vip3 in 52.8%, and vip1-2 in 34.4% of the 125 strains studied. Fifteen percent of the strains generated no PCR product for the three insecticidal toxin genes examined and produced less than 2 g/mL ␤-exotoxin I. The 45 isolates with supernatants toxic to S. littoralis harbored both the vip3 and cry1I genes, the products of which are known to be active against certain lepidopteran insects but not against coleopterans [8], and produced small amounts of ␤-exotoxin I (Ͻ5 g/mL). Twelve of the 13 isolates toxic to A. grandis produced significantly larger amounts of ␤-exotoxin I (Ͼ5 g/mL) and harbored vip1-2 genes. Three of these isolates were toxic to A. grandis (100% mortality in all cases) but not to S. littoralis (Ͻ15% mortality). We determined precisely the LC 50 of purified ␤-exotoxin I for both insects and found that A. grandis was five times more susceptible than S. littoralis to this toxin (LC 50 s of 6 g/mL and of 30 g/mL, respectively). The three strains toxic to A. grandis only produced intermediate amounts of ␤-exotoxin I (ϳ16 -50 g/mL), consistent with the observed pattern of mortality. Nine of the 10 isolates toxic to both insects produced more than 50 g/mL ␤-exotoxin I, accounting for their activity against both species. We found only one isolate toxic to A. grandis that did not produce ␤-exotoxin I or possess vip1-2 genes. This particular isolate contained a vip3like gene whose expression product was not toxic to A. grandis. This suggests that this strain secretes another, possibly novel, type of insecticidal compound highly active against A. grandis and/or S. littoralis.

Statistical analysis of toxicity and genotype data.

The observed frequencies of the various genotypes among the 125 Bt strains were compared to the expected frequencies of the 16 possible genotypes in the absence of linkage among these four factors (the toxin genes and ␤-exotoxin I production). We observed an overall significant difference (p Ͻ 0.0001) between the observed and expected frequencies, indicating linkage disequilibrium between the virulence factors. The sample only contained representatives of twelve of the sixteen possible genotypes (Fig. 1, black bars) and seven of them accounted for most of the variability. For all but two of these categories (vip3 ϩ cry1I ϩ ␤-exotoxin I and vip3 ϩ cry1I ϩ vip1-2), the deviations of the observed frequencies from the expected frequencies was significant (p Ͻ 0.01). These linkages-i.e., vip3 ϩ cry1I (23 strains), none of the parameters examined (19 isolates), ␤-exotoxin I (Ͼ2 g/mL) (17 isolates), all four toxins (13 isolates), and vip1-2 and ␤-exotoxin I (12 isolates)resulted in the frequencies of the most represented of the secreted virulence factors (i.e., Cry1I and Vip3), present singly in a strain, being much lower (p Ͻ 0.05 and p Ͻ 0.005 for cry1I and vip3, respectively) than expected.

Similarly, the frequency of ␤-exotoxin I production in strains with vip3 or cry1I genes was also lower than expected (p Ͻ 0.005). The other combinations were also underrepresented in the population. We used pairwise correlation analysis to investigate this differentiation further (Fig. 2). We found that there was a positive correlation between the vip3 and cry1I genes (r ϭ 0.81, p Ͻ 0.001) and between the presence of vip1-2 and ␤-exotoxin I production (r ϭ 0.3, p ϭ 0.002). This is consistent with our finding that toxicity to S. littoralis was almost always associated with the presence of vip3 and cry1I and that toxicity to A. grandis was associated with the presence of vip1-2 and ␤-exotoxin I production. Conversely, we found that ␤-exotoxin I production was negatively correlated with the presence of vip3 (r ϭ Ϫ0.24, p Ͻ 0.05) and with the presence of cry1I in a given strain (r ϭ Ϫ0.28, p Ͻ 0.01). We also analyzed the frequency of the various toxicity profiles as a function of the geographic origin or nature of the samples. The toxicity profile was not dependent (p Ͼ 0.05) on the geographic origin (Europe, Asia, Africa, or America) or the nature of the sample (invertebrates, dust, water, soil, or plants). However, strains producing more than 5 g of ␤-exotoxin I were preferentially isolated from invertebrates (this difference being marginally significant: p ϭ 0.05).

Ϫ Ϫ Ϫ 0 Ͻ g Ͻ 2 2 0 Ϫ Ϫ Ϫ 2 Ͻ g Ͻ 5 1 4 Ϫ Ϫ Ϫ 7 g 1 ϩ ϩ Ϫ 0 Ͻ g Ͻ 2 8 ϩ ϩ Ϫ 2 Ͻ g Ͻ 5 3 ϩ ϩ ϩ 0 Ͻ g Ͻ 2 4 ϩ ϩ ϩ 2 Ͻ g Ͻ 5 2 ϩ Ϫ ϩ 0 Ͻ g Ͻ 2 5 ϩ Ϫ Ϫ 0 Ͻ g Ͻ 2 4 Ϫ ϩ Ϫ 0 Ͻ g Ͻ 2 1 Ϫ Ϫ ϩ 0 Ͻ g Ͻ 2 2 Ϫ Ϫ ϩ 2 Ͻ g Ͻ 5 2 Ϫ ϩ ϩ 0 Ͻ g Ͻ 2 1 S. littoralis only 45 (36%) ϩ ϩ Ϫ 0 Ͻ g Ͻ 2 1 5 ϩ ϩ Ϫ 2 Ͻ g Ͻ 5 1 3 ϩ ϩ Ϫ 10 g 1 ϩ ϩ ϩ 0 Ͻ g Ͻ 2 7 ϩ ϩ ϩ 2 Ͻ g Ͻ 5 9 A. grandis only 3 (2.4%) Ϫ Ϫ ϩ 15 Ͻ g Ͻ 50 3 S. littoralis and 10 (8%) Ϫ Ϫ Ϫ Ͼ 50 g 1 A. grandis ϩ ϩ ϩ Ͼ 50 g 1 Ϫ ϩ Ϫ 0 Ͻ g Ͻ 2 1 Ϫ Ϫ ϩ Ͼ 50

Discussion

In this study, we investigated the prevalence and distribution of various insecticidal compounds, secreted by Bt, to obtain an insight into the biological diversity of natural isolates. Analysis of the toxicity spectra of the culture supernatants of the 392 Bt isolates revealed that over half the strains were not toxic to S. littoralis and A. grandis. This suggests that a large number of isolates do not secrete insecticidal compounds (at the concentrations tested) or that they produce toxic compounds active against other insect species. Among the toxic isolates those with activity against S. littoralis were the most frequent in this screening, consistent with previous observations that active Bt strains are generally more prevalent among lepidopteran species [START_REF] Chaufaux | Research on natural strains of Bacillus thuringiensis in different biotopes throughout the world[END_REF][START_REF] Iriarte | Environmental distribution and diversity of Bacillus thuringiensis in Spain[END_REF], the larvae of which are commonly found on leaves of many plants [START_REF] Smith | The phylloplane as a source of Bacillus thuringiensis variants[END_REF]. Analysis of ␤-exotoxin I production and search for the genes encoding the secreted Vip1-2, Vip3, and Cry1I toxins in 125 of these isolates indicated that these four toxins are widely distributed among natural Bt populations. The vip3 genes were detected in more than 50% of the strains tested; this figure was much higher than the 23% reported for the 125 colonies examined by Rice [START_REF] Rice | Specific primers for the detection of vip3A insecticidal gene within a Bacillus thuringiensis collection[END_REF], and this difference was highly significant (twotailed Fisher's exact test, p Ͻ 10 Ϫ5 ). Similarly, we found that 57% of the strains contained a copy of cry1I, whereas Gleave et al. [START_REF] Gleave | Screening by polymerase chain-reaction of Bacillus thuringiensis serotypes for the presence of cryV-like insecticidal protein genes and characterization of a cryV gene cloned from Bacillus thuringiensis subsp kurstaki[END_REF] found that the cry1I gene was present in 7 of 21 (31%) Bt serotypes, a frequency also significantly lower than that obtained in this study (twotailed Fisher's exact test, p ϭ 0.003). The reasons for these discrepancies are unclear, but may reflect natural variation in bacterial population structures from different locations. On the other hand, the frequency of ␤-exotoxin I production reported here, in 47.2% of the strains, although lower than the 59% (17 of 29 isolates) reported by Perani et al. [START_REF] Perani | Prevalence of beta-exotoxin, diarrhoeal toxin and specific delta-endotoxin in natural isolates of Bacillus thuringiensis[END_REF], is not significantly different (twotailed Fisher's exact test, p ϭ 0.307). Finally, we found that 34.5% of the isolates contained vip1-2 genes indicating that Vip1-2 toxins are also common among Bt strains. Given the large number of natural isolates producing one or several of these toxins, it is likely that these compounds make a significant contribution towards determining the insecticidal host range of individual Bt strains. Our results also indicated that the toxicity pattern of the isolates was independent of the nature of the sample (soil, insect, or other sources) and of its geographic origin. This is not very surprising since the ecology and environmental distribution of spore forming bacteria is probably largely due to aerial dispersion. However, we found a fairly good correlation between the genotypic profiles of the isolates and their insecticidal properties. It is therefore possible that some combinations of toxin genes of which the host can take a particular advantage may have become fixed. None of the 67 nontoxic strains produced significant amounts of ␤-exotoxin I, and the majority of them neither contained vip3, cry1I, or vip1-2. Similarly, all of the isolates toxic to A. grandis produced more than 5 g of ␤-exotoxin I, and the production of this compound in most of these strains was associated with the presence of vip1-2 genes. This suggests that toxicity to A. grandis is associated mainly with the production of ␤-exotoxin I but that both toxins may play an important role for virulence against A. grandis. We also found that all of the isolates toxic only to S. littoralis, harbored both vip3 and cry1I. It is not known whether the vip genes are present on the chromosome or on plasmids, but some of the genetic determinants involved in ␤-exotoxin I production have been shown to be present on the same plasmids as certain cry genes [START_REF] Gonzalez | Transfer of Bacillus thuringiensis plasmids coding for delta-endotoxin among strains of B. thuringiensis and B. cereus[END_REF][START_REF] Levinson | Identification of beta-exotoxin production, plasmids encoding beta-exotoxin and a new exotoxin in Bacillus thuringiensis by using high-performance liquid-chromatography[END_REF]. We found a strong linkage between vip1-2 and the genetic determinants of ␤-exotoxin I, on one hand, and between the cry1I and vip3 genes, on the other hand. This provides indirect evidence that the genes encoding these four insecticidal toxins are located on plasmids. We therefore suggest that vip1-2 and the genetic determinants of ␤-exotoxin I are generally present on the same plasmid and that cry1I and vip3 are generally present on another plasmid. The evolution and dissemination of adaptative traits such as toxin production in bacterial populations is frequently mediated by DNA elements such as plasmids [START_REF] Arber | Genetic variation: Molecular mechanisms and impact on microbial evolution[END_REF]. One possible driving force and mechanism responsible for the genetic shaping and divergence of natural Bt populations could be the transfer and fixation of advantageous plasmid/ toxin genes combinations that allow a given strain to become established as a new clonal population and colonize a special niche. However, further research is needed to establish whether there is indeed a correlation between the genetic differentiation of Bt populations and the presence of certain plasmid types or specific plasmidborne toxin genes combinations.

Fig. 1 .

 1 Fig. 1. Chi-square analysis of the linkage between various virulence factors in natural Bt isolates (*p Ͻ 0.05; **p Ͻ 0.01; ***p Ͻ 0.005; NS: not significantly different; exo, ␤-exotoxin I; v12, vip1-2; v3, vip3; cry, cry1I; fobs, frequency observed; ftheor, theoretical frequency).

* Ͼ 2 g

 2 of ␤-exotoxin I/mL.

Fig. 2 .

 2 Fig. 2. Pairwise correlations between virulence factors in natural Bt isolates (*p Ͻ 0.05; **p Ͻ 0.01; ***p Ͻ 0.005; NS, not significantly different; B, ␤-exotoxin I; v12, vip1-2; v3, vip3; cry, cry1I).

Table 1 .

 1 Oligonucleotides used in this study

	Primer	Sequence	Description, accession number	Product size (bp)
	cry1I-fw	5Ј aacactcagtatatgaat 3Ј	position 1059-1374, L49391	378
	cry1I-rev	5Ј cacatgtgatgctgaaat 3Ј	position 1437-1419, L49391	378
	vip2-fw	5Ј gggaaagaaaaagaaaaagagtgg 3Ј	position 231-255, A64136	914
	vip2-rev	5Ј catatatcctttgtcttcttt 3Ј	position 1145-1124, A64136	914
	vip3-fw	5Ј acatcctccctacactttctaatac 3Ј	position 1700-1725, L48811	678
	vip3-rev	5Ј tcttctatggacccgttctctac 3Ј	position 2378-2353, L48811	678

Table 2 .

 2 Insecticidal activity of the culture supernatants of natural Bacillus thuringiensis isolates against larvae of Spodoptera littoralis and Anthonomus grandis

				Number of strains with supernatants toxic to:
			Both
	Number of	Mortality	S. littoralis	S. littoralis	A. grandis
	supernatants	rate (%)	and A. grandis	only	only	Neither species
	Toxic supernatants	71-100%	29	107	5	-
	n ϭ 180	31-70%	4	32	3	-
	Nontoxic supernatants	6-30%	-	-	-	28
	n ϭ 212	0-5%	-	-	-	184
	Total n ϭ 392	-	33	139	8	212

n ϭ number of Bt strains.

Table 3 .

 3 Detection of cry1I, vip3, and vip1-2 genes by PCR analysis and quantification of ␤-exotoxin I in 125 isolates of Bacillus thuringiensis

		Total number	Genotype screening by PCR	Determination of ␤-
	Toxicity profile of	of isolates				exotoxin I g per mL	Number
	the supernatants	examined	cry1I	vip3	vip1-2	of culture supernatant	of isolates
	Nontoxic	67 (53.6%)

ACKNOWLEDGMENTS

We thank Didier Lereclus, in whose laboratory this work was conducted. We also thank Alex Edelman and Associates for editing the English manuscript. This work was supported by the Institut Pasteur and the Institut National de la Recherche Agronomique (INRA). We would like to thank Sophie Fleurde ´pine for technical assistance with PCR. Jeroen Van Rie is acknowledged for the gift of the primers used for PCR. Sylvain Espinasse was supported by a grant from Aventis Crop Science, Ghent, Belgium.