Eric Pante
email: pante.eric@gmail.com

Jawad Abdelkrim

Amélia Viricel

Delphine Gey

Scott France

Marie-Catherine Boisselier

Sarah Samadi

Boisselier Mc

Use of RAD sequencing for delimiting species

Keywords: (3-6): phylogenomics, octocoral, Stacks, PyRAD, SNP, species delimitation

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The advent of next-generation sequencing tools has permitted significant advances in our understanding of evolutionary processes such as speciation (e.g. [START_REF] Ekblom | Applications of next generation sequencing in molecular ecology of non-model organisms[END_REF], but some other practical applications of genomic data have been less explored, including phylogenomics and species delimitation. Among genomic approaches that are applicable to these fields, the usefulness of restriction-siteassociated DNA tag (RAD-tag; [START_REF] Baird | Rapid SNP discovery and genetic mapping using sequenced RAD markers[END_REF] sequencing has been investigated in few studies to date. This methodology typically provides short sequences (~ 100-150 bp) flanking the cut sites of a restriction enzyme (or several enzymes), generally yielding thousands of loci distributed throughout the genome. This approach does not require a reference genome, and can therefore be applied to non-model organisms.

However, some technical difficulties remain for groups where very little genomic knowledge is available (see [START_REF] Davey | Genomewide genetic marker discovery and genotyping using next-generation sequencing[END_REF]. For instance, the choice of restriction enzyme(s) and methodology (single-digest versus double-digest RAD) is key to estimating the number of expected cut sites and coverage, but relies on prior knowledge of genome size and GC content. Despite these difficulties, RAD-tag sequencing constitutes one of the reduced genomic approaches that are suitable for investigating inter-specific evolutionary questions. Published RAD-tag sequencing research beyond the species level includes in silico studies (Drosophila, mammals, and yeasts in [START_REF] Rubin | Inferring phylogenies from RAD sequence data[END_REF] Drosophila in [START_REF] Cariou | Is RAD-Seq suitable for phylogenetic inference? an in silico assessment and optimization[END_REF] and empirical work (e.g. Restionaceae flowering plants in [START_REF] Lexer | Next generation' biogeography: towards understanding the drivers of species diversification and persistence[END_REF]cetaceans in Viricel et al., 2014), which both suggest this approach is promising for taxa having diverged up to 60 million years ago. For instance, RAD-tag sequencing has proven useful in species delimitation and phylogenies within recently and rapidly diverged groups (e.g. Orobanchaceae flowering plants in [START_REF] Eaton | Inferring phylogeny and introgression using RADseq data: an example from flowering plants (Pedicularis: Orobanchaceae)[END_REF]swordtails in Jones et al., 2013; Heliconius butterflies in [START_REF] Nadeau | Genome-wide patterns of divergence and gene flow across a butterfly radiation[END_REF]cichlids in Wagner et al., 2013;geckos in Leaché et al., 2014). Comparatively, reconstructing the phylogeny of more distantly related taxa has been the topic of a single study (Carabus beetles, [START_REF] Cruaud | Empirical Assessment of RAD sequencing for interspecific phylogeny[END_REF], to the best of our knowledge. Herein we use this approach on a group of deep-sea octocorals for which little genomic data are available. Thus, our contribution constitutes one of the first studies investigating the use of RAD-tag sequencing for practical species delimitation within a taxonomic group composed of divergent species (up to 16 million years ago).

Deep-sea octocorals are one of the groups for which RAD-tag sequencing can significantly advance our understanding of evolutionary patterns. As for shallowwater octocorals, deep-water octocorals present significant challenges for taxonomists, with few morphological characters being available for species delimitation (e.g., [START_REF] Mcfadden | Molecular phylogenetic insights into the evolution of Octocorallia: A review[END_REF]. In addition, several studies have shown conflicting patterns of morphological and molecular data [START_REF] France | Genetic analysis of bamboo corals (Cnidaria: Octocorallia: Isididae): does lack of colony branching distinguish Lepidisis from Keratoisis?[END_REF][START_REF] Dueñas | Character lability in deep-sea bamboo corals (Octocorallia, Isididae, Keratoisidinae)[END_REF][START_REF] Pante | Pseudochrysogorgia bellona n. gen. n. sp.: a new genus and species of chrysogorgiid octocoral (Coelenterata: Anthozoa) from the Coral Sea[END_REF], suggesting that an integrative approach to taxon delimitation must be applied in this group (e.g. [START_REF] Schlick-Steiner | Integrative taxonomy: a multisource approach to exploring biodiversity[END_REF].

Octocorals, as with other anthozoans (e.g. scleractinians and sea anemones), are also plagued with remarkably low levels of mitochondrial genome evolution that renders the use of classical barcoding gene regions such as cox1 of limited use [START_REF] Mcfadden | Limitations of mitochondrial gene barcoding in Octocorallia[END_REF]. Comparatively, a few studies have successfully used nuclear markers within octocoral species (e.g. [START_REF] Concepcion | An alternative to ITS, a hypervariable, single-copy nuclear intron in corals, and its use in detecting cryptic species within the octocoral genus Carijoa[END_REF][START_REF] Mokhtar-Jamaï | From global to local genetic structuring in the red gorgonian Paramuricea clavata: the interplay between oceanographic conditions and limited larval dispersal[END_REF], but these are either not widely useable across octocorals (e.g. SRP54; France and Pante unpublished observations), or not informative at multiple phylogenetic scales (e.g. microsatellites). Multi-copy markers have been employed (e.g. [START_REF] Herrera | Molecular systematics of the bubblegum coral genera (Paragorgiidae, Octocorallia) and description of a new deep-sea species[END_REF], however their use implies that lack of concerted evolution within and across genomes will not blur the phylogenetic signal [START_REF] Vollmer | Testing the utility of internally transcribed spacer sequences in coral phylogenetics[END_REF][START_REF] Calderón | Evaluation of the utility of COI and ITS markers as tools for population genetic studies of temperate gorgonians[END_REF].

In this group, RAD-tag genotyping may therefore offer a panel of markers to help describe patterns of population structure, delimit species, and investigate phylogenetic relationships. This technique may however be difficult to implement in this group. Indeed, the composition of the deep-sea octocoral genome is unknown (size, GC content, prevalence of cut sites for restriction enzymes, etc.); the size of known cnidarian genomes, for instance, varies between 224 Mb and 1.8 Tb (Animal Genome Size Database; [START_REF] Gregory | Animal genome size database[END_REF]. In addition, sampling of deep-sea animals can be associated with a loss of quality of genomic DNA samples, particularly when sampling in tropical waters using trawls or dredges.

The genus Chrysogorgia (Calcaxonia: Chrysogorgiidae) is a noteworthy model for testing the utility of RAD sequencing for delimiting octocoral species, as it is diverse (62 nominal species described, 93% of which were based solely on morphology), widely distributed, and can be locally abundant [START_REF] Watling | Biology of deep-water octocorals[END_REF].

The large geographic, bathymetric, and ecological distributions of some Chrysogorgia species [START_REF] Pante | Deep-sea origin and in-situ diversification of chrysogorgiid octocorals[END_REF] question whether taxa are appropriately delimited, and whether cryptic diversity is important in the group. In the northwestern Atlantic, congruence exists between morphological and genetic data, suggesting that a relatively short fragment of the mitochondrial mtMutS gene can be used to formulate "Primary Species Hypotheses" [START_REF] Pante | Chrysogorgia from the New England and Corner Seamounts: Atlantic -Pacific connections[END_REF]. It is suspected that little to no intra-specific variation exists for this marker within the group [START_REF] Mcfadden | Limitations of mitochondrial gene barcoding in Octocorallia[END_REF], but the null hypothesis that single mutations at mtMutS are diagnostic of species limits must be evaluated using genetic data from markers informative within and above the species level. RAD loci allow to test whether lineages that putatively belong to different species do not exchange genes.

In this communication we test the utility of RAD-tag genotyping for delimiting species in Chrysogorgia using the genealogical criterion defined by Taylor et al., (2001). More specifically, we test whether single mutations on the mitochondrial mtMutS gene can be used as a criterion for grouping Chrysogorgia colonies into separate, putative species (or, more specifically, "Primary Species Delimitation hypotheses" as in [START_REF] Puillandre | Large-scale species delimitation method for hyperdiverse groups[END_REF]. We compare the results from two analysis pipelines, Stacks [START_REF] Catchen | Stacks: an analysis tool set for population genomics[END_REF] and PyRAD [START_REF] Eaton | PyRAD: assembly of de novo RADseq loci for phylogenetic analyses[END_REF], which significantly differ in the method employed for detecting homologous loci.

Material and methods

Specimen collection and mtDNA typing

Chrysogorgia specimens were collected from the SE slope of New Caledonia (NC) and adjacent seamounts of the Norfolk Ridge (82 colonies; Terrasses cruise, 2008), from Papua New Guinea (PNG; 8 colonies; BioPapua cruise, 2010), and from the northwestern Atlantic (1 colony, Extreme Coral 2010 cruise; Tables 1 andS1). Pacific specimens were retrieved from dredges and trawls (details on cruises of the Tropical Deep Sea Benthos research program: [START_REF] Bouchet | A quarter-century of deep-sea malacological exploration in the South and West Pacific: Where do we stand? How far to go?[END_REF]; details on the BioPapua cruise: Pante et al., 2012a); the Atlantic specimen was collected using the Jason II ROV (Woods Hole Oceanographic Institution). Specimens were fixed in 80% ethanol as soon as possible after collection. Genomic DNA was extracted using a CTAB protocol according to [START_REF] France | DNA sequence variation of mitochondrial large-subunit rRNA provides support for a two-subclass organization of the Anthozoa (Cnidaria)[END_REF]. A 700-bp fragment of the mitochondrial mtMutS gene (identified as more informative than cox1 or 18S in chrysogorgiids, [START_REF] Pante | Deep-sea origin and in-situ diversification of chrysogorgiid octocorals[END_REF] was amplified using the ND4L2475F -MUT3458R primer pair and sequenced using an ABI PRISM (R) 3100 or 3130xl Genetic Analyzer (primer information, PCR and sequencing conditions: [START_REF] Pante | Deep-sea origin and in-situ diversification of chrysogorgiid octocorals[END_REF]. Sequences were checked for quality and edited in Sequencher (TM) 4.7 (Gene Codes), aligned by eye (a single, 3 bp indel was present in the alignment), and haplotypes were submitted to GenBank (Table S1).

Divergence times among putative species were estimated using the molecular clock from [START_REF] Lepard | Analysis of variation in the mitochondrial encoded msh1 in the genus Lepidogorgia (Cnidaria: Octocorallia) and implications for population and systematics studies[END_REF], which was calculated for the shallow-water octocoral genus Lepogorgia based on mtMutS genetic distances for clades located on either sides of the Isthmus of Panama (0.14-0.25%/million years).

Library construction, RAD sequencing, and quality control

Genomic DNA quality was evaluated by 1% agarose gel electrophoresis, and quantified using a Thermo Scientific Nanodrop ND-1 000 spectrophotometer. DNA was sent to Eurofins Genomics (Ebersberg, Germany) for RAD-tag library preparation and sequencing. Libraries were constructed from 1-2 µg of DNA per colony using the SbfI restriction enzyme. This enzyme was chosen because it was successfully used in RADseq experiments with marine invertebrates (sea-anemones, [START_REF] Reitzel | Going where traditional markers have not gone before: Utility of and promise for RAD sequencing in marine invertebrate phylogeography and population genomics[END_REF]abalone, Gruenthal et al., 2014), and was expected to allow an acceptable compromise between prevalence of cut sites and depth of coverage, based on RADcounter (the University of Edinburgh, https://www.wiki.ed.ac.uk/display/RADSequencing/Home).

As the genome size and GC content of Chrysogorgia (or other octocorals, to the best of our knowledge) are not known, we estimated the prevalence of SbfI cut sites based on a range of genome sizes and GC content, based on information from the Animal Genome Size Database (see Introduction) and with a GC content of 40% (e.g. Soza-Ried et al., 2009). Barcodes 6-9 nucleotides long and differing by at least 2 nucleotides were used to differentiate multiplexed samples (Table S1). Sequencing was performed on two lanes of the Illumina (R) HiSeq (TM) 2 000 instrument (Illumina Inc., San Diego CA, USA) using the single read, 100 nucleotide configuration. Raw HiSeq output was processed using the CASAVA v1.8.2 software pipeline (Illumina Inc., San Diego CA, USA), and de-multiplexed and quality filtered using the process_radtags.pl module (default quality settings) of the Stacks v.0.99994 pipeline [START_REF] Catchen | Stacks: an analysis tool set for population genomics[END_REF]. A single sequencing error was tolerated in the barcode. Reads were truncated to 91 nucleotides. Quality (as measured by phred scores and percentage of sequence overrepresentation) was checked before and after treatment by process_radtags using FastQC v.0.10.1 (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/).

Exploration of the divergence parameter space

Two main pipelines specifically designed for analysis of RADseq data are currently available. The most used to date is the Stacks pipeline. It constructs a catalog of loci for a set of samples mainly based on three parameters: the minimum stack depth parameter m (i.e. the minimum number of reads allowed per allele), the intraindividual divergence parameter M (i.e. the maximum number of mutations that can be observed between stacks within a sample), and the inter-individual divergence parameter n (i.e. the maximum number of mutations that can be observed between loci across samples).

PyRAD [START_REF] Eaton | PyRAD: assembly of de novo RADseq loci for phylogenetic analyses[END_REF]) is a more recently developed pipeline and differs from Stacks in several ways, the most important one being that it allows the presence of indels, since the clustering process of reads into loci uses alignment tools. This is anticipated to be an advantage compared to the first pipeline when considering more phylogenetically distant species. PyRAD relies on a large number of parameters used at different steps of the process. Most of them are related to reads quality control, detection of homology and filtering of paralogs. Two main parameters are of particular importance: the minimum depth coverage Mindepth (minimum depth necessary to make a statistical base call at each position of a cluster) and the similarity threshold Wclust (similarity value to be used for the alignment during both the within and across-sample clustering).

For both pipelines, these parameter settings are expected to influence greatly the number of markers available for intra-and inter-specific comparisons and it is necessary to explore which parameter combinations maximize the number of orthologous loci [START_REF] Viricel | Applicability of RAD-tag genotyping for interfamilial comparisons: empirical data from two cetaceans[END_REF]. To explore the effect of these parameters at different phylogenetic depths, we randomly selected pairs of specimens that (1) were separated by 0 to 16 mutations at mtMutS (representing different levels of phylogenetic divergence), and (2) were characterized by 1 to 1.5 million reads (to alleviate potential effects of depth of coverage on the number of assembled loci). For each level of divergence, we used three replicate pairs of specimens. We refer to specimens with mtMutS haplotypes differing by few mutations as pairs of closelyrelated colonies, and those with haplotypes differing by many mutations as distantlyrelated colonies.

In Stacks, m was kept to 3 (the default value); M was incremented from 1 to 10 in two cases (specimens separated by 0 and 12 mutations at mtMutS), and from 1 to 7 in all other cases. Similarly, n was incremented from 1 to 10 (0 and 12 mutations cases) and from 1 to 8 (all cases). All combinations of M and n were not tested: only similar values of M and n were used together (two settings were used: M=n and M+1=n), as to (1) keep maximum levels of intra-and inter-individual divergence levels close, and (2) keep the number of Stacks analyses to a reasonable number. A total of 408 Stacks catalog construction tests were therefore performed using the denovo_map.pl script available in Stacks. Catalogs were parsed with the populations.pl script, where each sample was considered as a separate population, no missing data were allowed, and a minimum of 10 reads per SNP was set.

In PyRAD v. 2.0, combinations of two values for Mindepth (3 and6) and3 values for Wclust (0.89, 0.93 and 0.96) were tested, resulting in 156 analyses. For these analyses, the maximum number of sites per read with a quality < 20 (NQual) was set to 4, the minimum number of samples in a final locus (MinCov) was set to 1 and the maximum proportion of shared polymorphic sites in a locus (MaxSH) was set to 10%. For this last parameter, which aims at detecting paralogs, preliminary tests showed that in our case, changing this value did not drastically affect the number of loci and SNPs detected. Finally, optional parameters were kept to default values.

Comparison of Stacks and PyRAD

To evaluate what proportion of loci was detected by both PyRAD and Stacks, a custom BLASTN search was performed (BLAST toolkit v. 2.2.25; [START_REF] Zhang | 93% divergence Intersect (%) N. loci (pyRAD) N. loci[END_REF]. Local BLAST databases were constructed using PyRAD sequences (locus file containing consensus sequences for each individual; PyRAD parameters m=6 and Wclust=93% and 89%) for three groups of specimens with different numbers of reads (Table 2).

Stacks loci for these specimens (based on the locus file produced by the populations script, for which a single allele was retained per locus; denovo_map parameters m=3, M=4, n=4, and m=3, M=10, n=12) were then compared to the PyRAD database using BLASTN (percent identity set to 93% and 89%, word size 80 and 84 nt, ungapped alignments). The XML output of BLASTN searches was then parsed in bash using grep.

Phylogenetic reconstruction and species delimitation

RAxML v. 8.0.9 [START_REF] Stamatakis | RAxML-VI-HPC: Maximum likelihood-based phylogenetic analyses with thousands of taxa and mixed models[END_REF][START_REF] Stamatakis | A rapid bootstrap algorithm for the RAxML web servers[END_REF] was used on the CIPRES Portal [START_REF] Miller | Creating the CIPRES Science Gateway for inference of large phylogenetic trees[END_REF] to infer phylogenetic relationships among Chrysogorgia colonies, based on mitochondrial and nuclear sequences, using the GTRCATI model and automating boot-stopping. The mitochondrial phylogeny was inferred from the first 700 nt of the mtMutS gene (see above); the nuclear phylogeny was inferred using concatenated RAD loci obtained based on two parameter sets in Stacks, and one parameter set in PyRAD. The first Stacks set ("m3M4n4", denovo_map parameters m=3, M=4, n=4; populations script parameters m=6, p=2, r=0.5) corresponds to parameters that maximize the total number of loci detected while minimizing divergence parameters (see "Exploration of the divergence parameter space" section above). For this analysis, each mtMutS haplotype was considered as a separate population. The Stacks populations script parameters that were used signify that 50% missing data were allowed within each population, a locus had to be present in at least two populations to be included in the output and a minimum of 6 reads per SNP was required. The second Stacks set ("m3M10n12", Stacks script denovo_map parameters m=3, M=10, n=12; populations script parameters m=6, p=2, r=0.5) allowed more divergence between loci. The PyRAD dataset ("m6s93") was constructed with m=6 and Wclust=93% (details above). In all analyses, the Atlantic colony JAC1018 was used as the outgroup.

Once clades were delimited with RAxML, a Discriminant Analysis on Principal Components (DAPC, [START_REF] Jombart | Discriminant analysis of principal components: a new method for the analysis of genetically structured populations[END_REF] was used to explore genetic structure within three clades represented by 18 to 31 colonies (see below). This method takes into account the multilocus genotype of each individual and forms clusters based on genetic similarity without considering a model of evolution. We also used TESS [START_REF] Durand | Spatial inference of admixture proportions and secondary contact zones[END_REF] to investigate population structure using the conditional auto-correlative (CAR) admixture model with a spatially explicit, Bayesian framework. In TESS, the Deviance Information Criterion (DIC) was used to compare population structure in the presence of different numbers of clusters (the maximum number of cluster K was set to the total number of individual in the tested clade; for example, K was set from 2 to 18 for clade 1). Five replicate runs were used per K, with 1 200 MCMC steps and a 200-step burnin. The best K was determined by minimizing DIC and its variance; once the best K determined, a longer analysis with 12 000 steps and a 2 000-step burnin was run to obtain reliable individual assignments. The populations script in Stacks was re-run to keep only one SNP per locus, in order to minimize the probability of co-analyzing linked markers. The Stacks m3M4n4 dataset was chosen for these analyses for two reasons: (1) the DAPC and TESS analyses are run within clades at shallow phylogenetic depths, and (2) as only one SNP / locus is retained, divergence level should be kept minimal to prevent the inclusion of nonhomologous loci. The DAPC analysis was run using adegenet in R [START_REF] Jombart | adegenet: a R package for the multivariate analysis of genetic markers[END_REF]R Development Core Team 2014).

Results

Mitochondrial typing and RAD-tag sequencing

A total of 12 mtMutS haplotypes were detected among the 91 colonies investigated, 10 of which were from NC, 3 from PNG, 1 from the northwestern Atlantic, and 2 being shared between NC and PNG. The biogeography of these mitochondrial haplotypes at these locations is further discussed in Pante et al. (2012ab). A total of 236 million raw reads, corresponding to 35 463 Mbp were produced on two HiSeq2000 lanes. The number of quality-filtered reads (in millions) per colony varied between 0.04 (TER11108) and 5.82 (TER2044), with a median of 1.6. There was a significant correlation between the number of quality-filtered reads per colony and haplotypes (Kruskal-Wallis chi-squared = 25.11, df = 13, p-value = 0.02), haplotypes 6 and 10, for instance, yielded fewer reads than other haplotypes (haplotype 10 colonies were sampled from depths down to 880 m, and haplotype 6 colonies had remarkably small polyps that may have been particularly sensitive to prolonged times to preservation).

Loci, SNPs, and indel cataloguing using Stacks and PyRAD

Results from both pipelines (Stacks and PyRAD) show variations in the number of loci and SNPs depending on the set of parameters used (Figure 1a-e, 1g-k), as well as the mitochondrial genetic distance between samples (Figure 1f). For Stacks, as the mitochondrial genetic distance among included samples decreases, both the total number of loci and the number of polymorphic loci increases (Figure 1ab). The former ranges from a few loci to more than 2 000, whereas the latter ranges from a few loci to ~1 000, depending on the set of parameters used. When related to time of divergence (in MY, based on mtDNA), the total number of loci obtained decreases exponentially (Figure 1f). Inversely, the percentage of polymorphic loci is lower for more closelyrelated colonies (~40%) than for distantly related-colonies (~90%; Figure 1c). These three measures (number of loci, number of polymorphic loci and percentage of polymorphic loci) show the same response to an increase in divergence parameters M and n, namely a rapid increase followed by a plateau. This plateau is reached for the m3M4n4 set of parameters. Conversely, the number of SNPs increases drastically without reaching a plateau, from a few SNPs for the most stringent set of parameters and the most distantly-related colonies to around 3 000 for the most closely-related colonies and the most relaxed set of parameters (Figure 1d). Thus, the effect of increasing mitochondrial genetic distance among samples or decreasing stringency of parameters is to increase SNPs densities, from one SNP every 250 bp to one SNP every 20 bp (Figure 1e).

Results of the PyRAD analyses follow the general trends observed for the Stacks pipeline. These trends are an increase in total number of loci and polymorphic loci (Figure 1gh) for more relaxed parameters sets, as well as for more closely-related colonies. As for Stacks, more distantly-related specimen pairs have fewer loci than for closely-related ones, but a larger proportion of those is polymorphic (Figure 1i). While the percentage of polymorphic loci shows similar ranges of values for Stacks and PyRAD, the total number of loci as well as the number of polymorphic loci are almost doubled (from 2 000 to almost 4 000 and from 1 000 to almost 2 000, respectively). The same pattern is observed for the number of SNPs and SNP densities (Figure 1jk): PyRAD output differs from Stacks output by a factor of almost two, resulting in SNPs densities twice as high (from one SNP every 130 bp to one SNP every 20bp). Finally, unlike Stacks, PyRAD allows for indels within loci. The percentage of loci containing indels increases with less stringent sets of parameters (Figure 1l). Depending on the pair of samples considered, this measure varies from a few percent to almost 40 %.

For PyRAD, the number of catalogued loci decreased rapidly with the number of specimens included in the analysis (with significant drops corresponding to the number of individuals in the haplotype clades revealed by the phylogenetic reconstruction, see below) (Figure 2). Most loci bore <3 SNPs even when 10 polymorphisms were allowed on a single RAD locus (Figure 2).

We measured the proportion of loci catalogued by Stacks that was also detected by PyRAD using custom BLASTN database searches. Overall, 0.6 to 42.7% of loci detected by Stacks were present in the PyRAD catalog. This pattern is partly explained by the proportion of PyRAD loci with indels (see above), but might also be influenced by the differential detection of repeated regions (i.e. deleveraging algorithm in Stacks), or the number of reads per individual (the proportion of loci in common between Stacks and PyRAD was lower for individuals with fewer reads; Table 2).

Phylogenetic reconstruction and species delimitation

The automatic boot-stopping method implemented in RAxML yielded 1 000 bootstrap replicates for the mitochondrial phylogeny (91 taxa x 700 nt), 500 replicates for the Stacks RAD phylogenies (91 taxa x 1 080 352 nt, 11 872 loci for the first dataset, and 1 146 054 nt, 12 594 loci for the second dataset), and 200 replicates for the PyRAD phylogeny (91 taxa x 6 120 523 nt, 69 851 loci). The proportion of gaps and undetermined characters ranged between 83 and 84% for Stacks and was 92% for PyRAD. The three RAD phylogenies were similar but not identical, the second Stacks dataset being better resolved than the first, and the PyRAD dataset being better resolved than the Stacks sets (nodes with bootstrap >70%: 19% for m3M4n4, 29% for m3M10n12, 40% for m6s93; Figure 3). Divergence levels were much higher in the RAD phylogenies compared to the mitochondrial phylogeny. For instance, the groups composed of haplotypes 9 and 10 were separated by a distance of 0.001 substitution/site on the mtMutS tree, while these clades were separated by 0.27 and 0.25 substitutions/sites on the m3M4n4 and m3M10n12 RAD phylogenies, respectively (Figure 3).

Out of nine mitochondrial haplotypes represented by more than one individual, six formed well-supported monophyletic groups on the RAD phylogenies, for all datasets. One of these clades (corresponding to haplotype 10) contained specimens from both NC and PNG. The group formed by mitochondrial haplotype 7 was polyphyletic on the RAD phylogenies, with specimens grouping in two wellsupported clades on the PyRAD phylogeny: one composed of five closely-related NC specimens and one composed of three more divergent PNG colonies (this clade was split in two on the Stacks phylogenies). Specimens characterized by mtMutS haplotype 7 may therefore belong to at least three distinct species. On the other hand, specimens characterized by three distinct mitochondrial haplotypes (2, 8, 13) clustered into a single, well-supported clade (with the exception of one individual, TER13034, haplotype 8, which clusters well outside this clade). These three haplotypes, which form a paraphyletic group on the mitochondrial phylogeny and are one to two mutations different from each other, would therefore be considered as one evolutionary unit based on the RAD phylogenies (and population clustering analyses with DAPC and TESS failed to detect structure within this clade; see below). Finally, out of three singleton haplotypes (J, 13, 14), two (J, 14) sit on long branches and are clearly differentiated from other haplotypes using RAD-tag data.

We ran a DAPC on the three clades that contained the most colonies (clade 1: 18 colonies of haplotype 9; clade 2: 20 colonies of haplotype 4; clade 3: 31 colonies of haplotypes 2, 8, 13). Within these clades, 3 685, 1 470 and 8 201 loci were retained (with 25, 42 and 55% missing data, respectively). In all three cases, DAPC failed to detect intra-clade genetic structure, as the most likely number of group (based on BIC, discounting the scenario in which each sample belongs to its own group), in each case, was one (Figure S1). The spatially-explicit admixture model implemented in TESS also failed to detect genetic structure within clades 1 and 3, but suggested the presence of three clusters in clade 2, these clusters being composed of colonies sampled (1) on the slope of New Caledonia, (2) Munida Seamount (Norfolk Ridge), and (3) Jumeaux Ouest Seamount (Norfolk Ridge; Figure S1). The population genetics of Chrysogorgia will be further discussed in a separate study.

Detection of environmental contaminants

As octocoral DNA was extracted from whole polyps rather than dissected, internal tissue, some loci may come from environmental contaminants such as bacteria. To evaluate the prevalence of such loci, we blasted all the loci that were catalogued for the m3M4n4 Stacks dataset from individual JAC1018 (n = 1 202). The BLASTN algorithm [START_REF] Altschul | Gapped BLAST and PSI-BLAST: a new generation of protein database search[END_REF] was used to match RAD loci to the nonredundant NCBI nucleotide database, using 10 -3 as a statistical significance threshold (e-value). Most sequences (92.6%) could not be assigned to a match in the nucleotide database and 4.5% of loci were similar to bacterial sequences (78-100% similarity between match and query). A single locus matched human mitochondrial DNA (84% similarity); other matches (n = 34) included other invertebrates and plant sequences.

Given (1) the small prevalence of potential contaminants, (2) our inability to determine whether these loci really belong to contaminant DNA or correspond to coral sequences which closest matches are non-cnidarian taxa, and (3) the large number of Stacks analyses performed (>400), we decided to run our analyses without trying to filter loci from exogenous DNA sources.

Discussion

A critical decision in RAD analyses is the way the sequencing data are filtered to get to the final SNP dataset. This process goes through several steps to ensure that the final loci will correspond to homologous sequences. The main filters involve several quality filters (sequencing quality, sequencing depth) as well as several similarity thresholds aimed at identifying the different allelic states of homologous loci. Finally, for each sample, an algorithm is used to tell apart sequencing errors from real mutations in order to conduct the final SNP calling. Even though the overall process is quite similar for Stacks and PyRAD analyses pipelines, a strict comparison of their results is not straightforward since they use sets of parameters that differ to some extent. A main difference between these two pipelines is in the assessment of similarity of loci: Stacks uses a strict similarity criterion (maximum number of mutations) in order to cluster reads into loci, whereas PyRAD uses an overall similarity criterion, after an alignment step, allowing for the presence of indels within clusters. This should be a critical difference when comparing genetically more-distant samples as indels are more likely to occur, and would thus result in sequences being assigned to different loci using Stacks (which will then be excluded from the final catalog since not present in both individuals) while PyRAD would theoretically allow these reads to be considered as homologous loci.

Our results show that more loci are recovered using the PyRAD pipeline. Despite these differences, general trends are similar using both pipelines. First, fewer loci and SNPs are recovered when comparing more genetically distant samples. This result is expected and has been anticipated through simulation [START_REF] Cariou | Is RAD-Seq suitable for phylogenetic inference? an in silico assessment and optimization[END_REF] and observed empirically [START_REF] Cruaud | Empirical Assessment of RAD sequencing for interspecific phylogeny[END_REF]. Our data show an exponential decay of the number of loci recovered as a function of divergence time of samples. Second, the stringency of the filtering process has a significant effect on the number of loci and SNPs identified. Indeed, higher minimum depth of sequencing thresholds and higher similarity threshold lead to fewer loci being identified. This trend is observed regardless of the level of genetic divergence between samples, but it seems to be accentuated when samples are more closely related.

Despite the similarities in general trends, quantitative and qualitative differences are observed in the outputs of each pipeline. Indeed, whatever the set of parameters used, almost twice as many loci are identified using PyRAD compared to Stacks. This difference cannot be solely attributed to the management of indels since our results show that the percentage of loci containing indels is usually around 5-20% and never reaches 40% whatever the genetic distance between samples and the parameters set. Another interesting result is that PyRAD is not simply adding extra loci to the total loci identified by Stacks: only half of the loci identified using Stacks are also present in the PyRAD loci catalogs. It is thus necessary to invoke other filtering processes and differences in algorithm to explain these differences in output. More thorough analyses would be needed in order to identify precisely what are the main sources of divergence in the processing of raw data, in addition to the treatment of indels.

One major result is the remarkable loss of homologous loci with increasing divergence among specimens with different mitochondrial haplotypes. For instance, compared to specimens sharing the same haplotype, specimens two mutations apart at mtMutS (estimated divergence of 1-2 My) had on average 70% fewer homologous loci (Stacks analysis at m3M7n8). Within the genus, specimens from mitochondrial clades 16 mutations apart (i.e. the highest divergence level included in our study, estimated between 9 and 16 My) share 97% fewer loci. This rate of loss of homologous RAD-tags is far greater than what has been observed in cetaceans [START_REF] Viricel | Applicability of RAD-tag genotyping for interfamilial comparisons: empirical data from two cetaceans[END_REF], for which 66% of homologous loci were retained at the inter-familial level (short-beaked common dolphins, Delphinus delphis, vs. harbour porpoise, Phocoena phocoena; estimated divergence of 14-19 My) compared to the intra-specific level (within Delphinus delphis). Comparisons within cetaceans were performed using the same custom pipeline as used in the present study, using Stacks parameters m3M3n3 (the results for corals were similar when comparing m3m3n3 to m3M7n8).

The differences observed between our study and that of [START_REF] Viricel | Applicability of RAD-tag genotyping for interfamilial comparisons: empirical data from two cetaceans[END_REF] may be explained by various factors. For example, the choice of restriction enzyme was different (Sbf1 here, Not1 for Viricel et al.), and differences in genome composition (most importantly GC content and size) are unknown. While both studies were conducted with two lanes of Illumina HiSeq2000 sequencing (conducted by Eurofins Genomics in both cases), throughput may have been influenced by the quality of genomic DNA (trawled deep-sea samples here, stranded animals for Viricel et al.).

These various factors may have significantly influenced the number of cut sites. Our comparisons might also be significantly affected by the precision of the molecular clocks available. Divergence times between cetacean families were inferred based on fossil evidence (see references in [START_REF] Viricel | Applicability of RAD-tag genotyping for interfamilial comparisons: empirical data from two cetaceans[END_REF], while no such fossil-calibrated molecular clock exists, to the best of our knowledge, for octocorals. The mtMutS divergence rates estimated by [START_REF] Lepard | Analysis of variation in the mitochondrial encoded msh1 in the genus Lepidogorgia (Cnidaria: Octocorallia) and implications for population and systematics studies[END_REF] are based on a group of shallow-water octocorals that may evolve faster than the deep-sea Chrysogorgia (a long standing question in deep-sea biology is whether evolutionary process take longer in deeper water, compared to shallower waters; e.g. [START_REF] Wilson | Speciation in the Deep Sea[END_REF], and rely on a geological event (rising of the Isthmus of Panama), which can introduce further bias.

The exploration of divergence parameter space, as outlined above, was made using pairs of specimens, and not allowing any missing data. Stacks and PyRAD can build catalogs with loci shared by a set proportion of individuals within pre-defined groups. Hence, our phylogenetic matrix based on over 12K loci (Stacks parameters m3M10n12) resolved most deeper nodes of the tree despite 83 to 84% of missing data. Similarly, [START_REF] Cruaud | Empirical Assessment of RAD sequencing for interspecific phylogeny[END_REF] constructed a phylogeny of 18 species of the beetle genus Carabus, and found that the deepest node of the tree (17 My divergence between species) was characterized by 67% of missing data but strong statistical support. [START_REF] Jones | The evolutionary history of Xiphophorus fish and their sexually selected sword: a genome-wide approach using restriction site-associated DNA sequencing[END_REF] reconstructed phylogenetic relationships among congeneric species of swordtail and platyfish (Xiphophorus sp.) that diverged <3 My, and estimated up to 70% missing data (ingroup data; their Table S2). They noted, however, that missing data had little effect on tree topology and branch support. The rate of loss of homologous loci observed in swordtail and platyfish is more on par with what we observed for Chrysogorgia than what was reported for cetaceans and Carabus beetles, and further emphasizes that (1) the utility of RAD sequencing for phylogenetic reconstruction may be taxon-dependent, and (2) molecular clocks must be critically interpreted. It must be underlined, however, that notable differences in tree topologies were observed between the three inferred RAD phylogenies, such as deep but well-supported nodes (e.g. relative positions of clade 3 and haplotypes 6, 7 and 8).

RAD-tag sequencing has also proven very useful in testing the criterion used for our primary species delimitation hypotheses, namely that single mitochondrial mtMutS haplotypes discriminate species that fit within the General Lineage Concept of species as defined by de Queiroz (1998). Indeed, a large numbers of variable loci could be catalogued within and among closely-related colonies (sharing the same mtMutS haplotype, and therefore putatively belonging to the same species) and more distantly-related colonies (separated by 1-16 mutations at mtMutS, putatively belonging to different species), allowing us (1) to plot our primary delimitation hypotheses onto well-supported phylogenies, and (2) to explore the spatial structure of populations. Three patterns were evidenced from the data: (1) in the majority of cases we noted a complete congruence between mtMutS haplotypes and RAD clades (6/9 non-singleton haplotypes and 2/3 singleton haplotypes); (2) in one case incomplete congruence was noted (with PyRAD, haplotype 7 corresponding to two RAD clades (one NC, one PNG) that did not form a monophyletic group; (3) in one case a single RAD clade included specimens with different (but closely-related) haplotypes. This result is significant for octocoral taxonomy and systematics, as mtMutS has been widely used to assist species delimitation across a large number of families (e.g. review of [START_REF] Mcfadden | Molecular phylogenetic insights into the evolution of Octocorallia: A review[END_REF]. While morphological, mitochondrial [START_REF] Pante | Chrysogorgia from the New England and Corner Seamounts: Atlantic -Pacific connections[END_REF] and genomic data (this study) all point to the utility of mtMutS for delimiting Chrysogorgia species, its resolution should be interpreted in two ways.

First, as we did not find 100% congruence between RAD clades and mtMutS haplotypes, and tested only a restricted set of putative species, mtMutS should still be considered as one of the first steps in an integrative taxonomic loop incorporating more variable markers (e.g. [START_REF] Schlick-Steiner | Integrative taxonomy: a multisource approach to exploring biodiversity[END_REF][START_REF] Kekkonen | DNA barcode-based delineation of putative species: efficient start for taxonomic workflows[END_REF].

Second, the evolutionary speed of mtMutS may well vary among octocorals, and its resolving power may therefore vary from one group to another (e.g. [START_REF] Baco | Comparing molecular variation to morphological species designations in the deep-sea coral Narella reveals new insights into seamount coral ranges[END_REF]. Nevertheless, combining mitochondrial markers such as mtMutS and RAD-tag data will without doubt be of tremendous value for testing the large number of outdated species hypotheses within the Octocorallia. Chrysogorgia specimens. Wclust: percent divergence permitted between loci within and across specimens; in addition to the 93% Wclust level used to infer the Chrysogorgia phylogeny, the 89% Wclust level was tested here. sequences. The number of loci detected within nine individuals (with high, medium and low read numbers) is presented for the two analyses performed on the entire set of 91 specimens. The number of quality-filtered reads is given in million.

Titles and Legends to Figures

Table S1. Excel table with information on collection (location, date, coordinates, depth), mitochondrial haplotypes (haplotype number and GenBank accession number), and number of quality-filtered reads for the 91 Chrysogorgia specimens used in this study. The 6-9 nucleotide barcodes used to distinguish specimens after Illumina sequencing are also included.

Figure 1 .

 1 Figure 1. Comparison of locus detection for Stacks (a-f) and PyRAD (g-l). The number of loci, SNPs and indels detected for specimens separated by 0-16 mutations at the mitochondrial mtMutS gene are shown for the different read coverage (m parameter) and divergence levels (M and n parameters, see text). In PyRAD analyses, "s" corresponds to the "Wclust" parameter.

Figure 2 .

 2 Figure 2. Information content of the locus catalog built by PyRAD for all 91

Figure 3 .

 3 Figure 3. Maximum likelihood phylogenetic trees inferred using RAxML for the mitochondrial mtMutS data (a), and RAD loci (b-d). Bootstrap node support (1000 replicates for a, 500 replicates for b-c, 200 for d) is presented only for nodes with ≥ 70% support. At the tips, colored dots, which represent mtMutS haplotype membership (each color represents a unique haplotype), are followed by specimen identifiers and haplotype numbers. Each tree was rooted to the Atlantic specimen (JAC1018, haplotype J). Genetic structure within clades 1, 2 and 3 were further investigated using a DAPC and TESS (see text and Figure S1). Scale bars: substitution / site.

Figure S1 .

 S1 Figure S1. Population genetic structure within three clades of the phylogenetic analysis. a-c: Bayesian Information Criterium (BIC) values for each tested number of DAPC cluster. For each clade, the maximum number of clusters was set as the

Table 1 .

 1 Summary table of haplotype information (sample size, geographical spread, depth range, habitat (seamounts vs. slopes) and mtMutS vs. RAD delimitation. NC:

	New Caledonia, PNG: Papua New Guinea

Table 2 .

 2 Results of the BLASTN alignments performed between Stacks and PyRAD

divergence Specimen Haplotype read.category N. reads (M) N. loci (pyRAD) N. loci (Stacks)

	Haplotype N. colonies Geography Habitat J 1 Atlantic slope 2 11 NC slope 4 20 NC slope & seamoun 150 -Depth range (m) 627 -390 -Delimitation mtMutS / RAD congruence mtMutS / RAD incongruence mtMutS / RAD congruence 89% TER2044 11 high 5.82 6580
	6 mtMutS / RAD congruence 2 NC JAC1018 J high	seamount 5.49	270 -3305
	7 8 9 10 11 13 mtMutS / RAD incongruence 8 NC-PNG 20 NC 18 NC 3 NC 3 NC 1 NC TER7092 7 high mtMutS / RAD incongruence TER130424 9 median mtMutS / RAD congruence TER13064 8 median mtMutS / RAD congruence TER13087 9 median mtMutS / RAD congruence mtMutS / RAD incongruence TER11101 4 low	slope & seamoun 300 -slope 390 -slope 390 -slope & seamoun 458 -seamount 750 -slope 460 -4.04 6867 1.61 6151 1.61 6876 1.60 5959 0.09 1046
	14 30 mtMutS / RAD congruence 1 NC 3 PNG TER13047 9 low mtMutS / RAD congruence TER11108 4 low	slope slope 0.08 0.04	400 -220 -1145
		TER12047, 8 TER12047, 8
		TER120414, 8
			Clade 1
			Clade 3
			Clade 2

Acknowledgments

The authors thank the participants of the national research group "Génomique Environnementale" for stimulating discussions on the use of RAD-tags for inferring phylogenies, in particular R. Debruyne and A. Cruaud. Samples used in this study were collected during the Terrasses and BioPapua cruises (PIs S. Samadi and L. Corbari) as part of the MNHN/IRD Tropical Deep-Sea Benthos programme (Bouchet et al., 2008), and during the Extreme Corals 2010 cruise in the northwestern Atlantic (PIs SW Ross and SD Brooke; funding from NOAA's Deep Sea Coral Research and Technology Program). We warmly thank the participants and crew members on these cruises for they indispensable help at sea. This study was funded by the "Institut Écologie et Environnement" of the CNRS ("Appel { Projets en Génomique Environnementale," organizers Dominique Joly and Denis Faure), the French Muséum national d'Histoire naturelle (Action Thématique du Muséum "Taxonomie moléculaire: DNA Barcode et gestion durable des collections") and by the Agence National de la Recherche (ANR 12-ISV7-0005-01 French-Taiwanese project TF-DeepEvo). Parts of the analysis were run on the YMIR super-computer (partly funded by the European Union (contract 31031-2008, European Regional Development Fund) of the University of La Rochelle (many thanks to Mikael Guichard, Marc-Henri Boisis-Delavaud and Frédéric Bret). We also thank Julio Pedraza (UMS 2700, MNHN) for his help with bioinformatics. Salary for EP was covered by a grant to the Poitou-Charentes region (Contrat de Projet État-Région 2007-2013). Finally, the authors thank the editor and two anonymous reviewers for their constructive comments.

Conflict of Interest:

The authors declare no conflict of interest.

Data Archiving: Mitochondrial haplotypes were deposited on GenBank (Table S1).

Phylogenetic data were deposited on Dryad: doi:xxxxx.


Skeletal Quads: Human Action Recognition Using Joint Quadruples Georgios Evangelidis, Gurkirt Singh and Radu Horaud INRIA Grenoble Rhône-Alpes 655, avenue de l'Europe 38330 Montbonnot Saint-Martin, FRANCE email: firstname.lastname@inria.fr

Abstract-Recent advances on human motion analysis have made the extraction of human skeleton structure feasible, even from single depth images. This structure has been proven quite informative for discriminating actions in a recognition scenario. In this context, we propose a local skeleton descriptor that encodes the relative position of joint quadruples. Such a coding implies a similarity normalisation transform that leads to a compact (6D) view-invariant skeletal feature, referred to as skeletal quad. Further, the use of a Fisher kernel representation is suggested to describe the skeletal quads contained in a (sub)action. A Gaussian mixture model is learnt from training data, so that the generation of any set of quads is encoded by its Fisher vector. Finally, a multi-level representation of Fisher vectors leads to an action description that roughly carries the order of subaction within each action sequence. Efficient classification is here achieved by linear SVMs. The proposed action representation is tested on widely used datasets, MSRAction3D and HDM05. The experimental evaluation shows that the proposed method outperforms state-of-the-art algorithms that rely only on joints, while it competes with methods that combine joints with extra cues.

I. INTRODUCTION

Action recognition is an active topic in computer vision and pattern recognition, with many potential applications in human-computer interaction. Despite the advances in recent years, however, recognising human actions remains a challenging problem, mainly because of the articulated nature of human motion. Therefore, the discrimination of human postures and actions benefits from the segmentation of the body into parts. While this kind of segmentation remains a quite difficult task using monocular visual sensors, the release of depth sensors (e.g. Kinect) has simplified the pose estimation by means of 3D body joints [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF].

Typically, an action recognition method employs an L-class classifier that performs a 1-of-L assignment to input vectors. By putting aside the classifier itself, what mostly differentiates the majority of the methods is the way of building the input of the classifier, i.e., the vector representation of a video segment. Commonly, multiple descriptors of raw data are summarised into a vector, e.g. through a Bag-of-Words (BoW) paradigm, in order to encode an action sequence.

Thanks to recent achievements [START_REF] Shotton | Real-time human pose recognition in parts from single depth images[END_REF], the detection of the human pose by means of skeleton joints is feasible even from Support from the European Research Council (ERC) through the Advanced Grant VHIA (#340113) is greatly acknowledged. A Gaussian Mixture Model (GMM), learnt on training data, is supposed to generate skeletal quads. Based on the GMM parameters, the skeletal quads of any action example are encoded into a Fisher vector, thus building the action descriptor which is led to a multi-class linear SVM. single depth images. This implies a powerful representation for action recognition, since actions can be seen as a set of poses. However, when joints are pooled into global features, e.g., using all pairwise joint differences [START_REF] Yang | Eigenjoints-based action recognition using naivebayes-nearest-neighbor[END_REF], the articulated nature of the human pose is not well encoded in the action descriptor. Therefore, local skeleton features are more meaningful.

In this paper, we propose a compact yet effective local skeleton descriptor that makes the pose representation invariant to any similarity transformation, hence view-invariant. The novel skeletal feature, referred here to as skeletal quad, locally encodes the relation of joint quadruples, so that 3D similarity invariance is guaranteed. Unlike the common BoW paradigm, we adopt a Fisher kernel representation [START_REF] Jaakola | Exploiting generative models in discriminative classifiers[END_REF]. Inspired by [START_REF] Perronnin | Improving the fisher kernel for large-scale image classification[END_REF], we consider a Gaussian mixture model that generates the skeletal quads, while we enable a power normalisation for the Fisher vectors. Further, a multi-level splitting of sequences into segments is invoked to integrate the performing order of subactions into the vector representation. Such vectors constitute the input of a multi-class linear SVM. Fig. 1 illustrates our action recognition pipeline.

The remainder of the paper is organised as follows. We summarise the related work in Sec. II and we propose our video representation in Sec. III. Sec. IV in short discusses the used classifier, while our method is tested on public datasets in Sec. V. Finally, Sec. VI concludes this work. II. RELATED WORK Shotton et al. presented in [15] an articulated pose recognition algorithm that makes the extraction of skeleton structure from single depth images possible. This result inspired many researchers either to rely on skeleton information only or to combine joints with other depth and/or color cues in order to recognise actions. In what follows, we first describe methods that only use skeleton data, hence more related to our approach, and we proceed with the state-of-the-art methods that use multiple features. For a recent detailed survey on human motion analysis from depth data, we refer the reader to [START_REF] Aggarwal | Human activity recognition from 3d data: A review[END_REF], [START_REF] Ye | A survey on Human Motion Analysis from Depth Data. Time-of-Flight and Depth Imaging[END_REF].

Xia et al. [START_REF] Xia | View invariant human action recognition using histograms of 3d joints[END_REF] suggest a compact posture representation through a histogram of 3D joint locations. Linear discriminant analysis along with a BoW model translates each action into a series of symbols (quantized postures). Then, a generative classifier (discrete HMM), that deals with the temporal nature of data, classifies each input. Yang and Tian [START_REF] Yang | Eigenjoints-based action recognition using naivebayes-nearest-neighbor[END_REF] rely on position differences of in-frame and cross-frame joints. The resulting long vectors are compressed by PCA so that each frame is described by an EigenJoint descriptor. As for the classification, a naive Bayes nearest-neighbor classifier assigns to unknown inputs the label with the minimum video-to-class distance. Ofli et al. [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] represent an action instance as a sequence of the most informative joints (SMIJ) per action. This selection is based on joint related measures such as the moments of the joint angles. Several encoding methods are suggested for the vector representation of SMIJ, while two classifiers are tested, nearest-neighbor classifier (NNC) and SMVs. Chaundry et al. [START_REF] Chaudhry | Bioinspired dynamic 3d discriminative skeletal features for human action recognition[END_REF] are inspired by a neural encoding method and describe skeletal features based on a medial-axis template. Global linear dynamical systems (LDS) model the temporal dynamics of such features, while the LDS parameters describe the sequences and define the input of the classifier (either SVMs or NNC).

Regardless of the skeleton structure, the raw depth map itself provides a source for extracting discriminative features. Li et al. [START_REF] Li | Action recognition based on a bag of 3d points[END_REF] propose a BoW model to describe the points close to the silhouette, after their projection to the three Cartesian planes. Local occupancy patterns (LOP) of depth sequences have been also used as features. Random and spatio-temporal LOPs are proposed in [START_REF] Wang | Robust 3d action recognition with random occupancy patterns[END_REF] and [START_REF] Vieira | On the improvement of human action recognition from depth map sequences using spacetime occupancy patterns[END_REF] respectively. Wang et al. in [START_REF] Wang | Mining actionlet ensemble for action recognition with depth cameras[END_REF] combine LOPs around joints with joint differences to build joint-based time series, whose Fourier coefficients are used to describe the action. Moreover, a mining step provides a pool of informative actionlets (subset of joints) per action, that are taken into account by the classifier being used. Instead of the raw depth data, the 4D normals are used by Oreifej et al. [START_REF] Oreifej | Hon4d: Histogram of oriented 4d normals for activity recognition from depth sequences[END_REF], while their distribution is encoded in a histogram whose bins are irregularly spaced in 4D; this binning results from a learning process. Xia and Aggarwal [START_REF] Xia | Spatio-temporal depth cuboid similarity feature for activity recognition using depth camera[END_REF] recently presented a depth-based spatio-termporal detector along with a self-similarity feature that describes local depth areas of adaptive size. Finally, Ohn-Bar and Trivedi [START_REF] Ohn-Bar | Joint angles similiarities and hog 2 for action recognition[END_REF] track the joint angles and build a descriptor based on similarities between angle trajectories. This feature is further combined with a double-HOG descriptor that accounts for the spatio-temporal distribution of depth values around the joints. Both [START_REF] Xia | Spatio-temporal depth cuboid similarity feature for activity recognition using depth camera[END_REF] and [START_REF] Ohn-Bar | Joint angles similiarities and hog 2 for action recognition[END_REF] illustrate the benefits of combining different features.

III. SKELETON-BASED VIDEO REPRESENTATION

We propose new representation based on skeleton joints, namely a local skeletal descriptor and an associated representation that. This skeletal descriptor, which is referred to as skeletal quad, was inspired by a geometry hashing method that describes the positions of nearby stars in night sky images [START_REF] Lang | Astrometry.net: Blind astrometric calibration of arbitrary astronomical images[END_REF]. It was also used to describe keypoint constellations in video frames [START_REF] Evangelidis | Efficient subframe video alignment using short descriptors[END_REF]. Our action descriptor is a Fisher kernel representation which encodes the Jacobian of the probability function that generates the skeletal quads contained in a depthimage sequence.

A. Skeletal quads

Let x ∈ R 3 denote the coordinates of a skeleton joint in some world centered frame. Suppose also that we are given a quadruple of (nearby) joints X = [x 1 x 2 x 3 x 4 ], where (x 1 , x 2 ) is the most widely separated pair of points within the quadruple. We consider a local coordinate system, such that x 1 becomes the origin and x 2 is mapped onto [1, 1, 1] ⊤ . This constrains a similarity transformation (a 3×3 rotation matrix, a translation vector and a scale factor) whose parameters can be easily computed from x 1 and x 2 . Once these parameters are estimated, the quadruple is mapped onto its new coordinates:

S(x i ) = sR[x i -x 1 ], i = 1 . . . 4, (1)
with

S(x 1 ) = [0, 0, 0] ⊤ and S(x 2 ) = [1, 1, 1] ⊤ .
Hence the quadruple is encoded by six parameters, namely q = [S(x 3 ); S(x 4 )], where the notation [•; •] denotes vertical vector concatenation. We refer to this descriptor as skeletal quad. Fig. 2 shows an upper-body skeleton example and how the skeletal quad is formed for a joint quadruple.

Note that R should also normalize any possible rotation around the axis x 1 x 2 , e.g., by further mapping x 3 to a specific plane, that leads to a 5D descriptor. However, when there is no strong cross-subject variation in body orientation within each action category, it is sufficient enough to just map x 2 to [1, 1, 1] ⊤ , thus increasing the discriminability. 1Unlike representations of skeleton joints that rely on joint differences (e.g. [START_REF] Wang | Mining actionlet ensemble for action recognition with depth cameras[END_REF]), our descriptor is scale, viewpoint and body-orientation invariant. Moreover, this coding scheme leads to well distributed points in R 6 [START_REF] Lang | Astrometry.net: Blind astrometric calibration of arbitrary astronomical images[END_REF]. Note that there is a kind of symmetry between descriptors owing to the different order in the pairs x 1 , x 2 , and x 3 , x 4 , which can be easily broken or taken into account. Here, we consider both orders for the first pair, while the third point of X is always the closest to the local origin between the remaining points.

B. A Fisher Kernel representation

The superiority of Fisher vectors (FV) against the popular BoW representation has been analyzed in the image classification context [START_REF] Perronnin | Improving the fisher kernel for large-scale image classification[END_REF]. We follow a similar approach in order to describe an action sequence. It is important to note that the low dimension of the proposed descriptor compensates for the large inherent dimensionality associated with Fisher vectors.

p(Q|θ) = M i=1 K k=1 w k N (q i |µ k , σ k ), (2)
where θ = {w k , µ k , σ k }, k = 1, ..., K is the set of the mixture parameters with mixing parameters w k , means µ k ∈ R 6 and diagonal covariance matrices (represented here as vectors) σ k ∈ R 6 . These parameters can be easily estimated via the standard EM algorithm based on a training dataset. Once the GMM parameters are estimated, any set Q may be described by its Fisher score [START_REF] Jaakola | Exploiting generative models in discriminative classifiers[END_REF], namely the Jacobian of the logprobability with respect to the GMM parameters:

J Q θ = ∇ θ log p(Q|θ) . (3)
The Fisher kernel, however, relates any two such vectors through a bilinear form based on the inverse of the Fisher information matrix. Since the decomposition of this matrix is possible, one can write the Fisher kernel as a linear kernel of the so called Fisher vectors, denoted here as J . The reader is referred to [START_REF] Jaakola | Exploiting generative models in discriminative classifiers[END_REF] for a detailed analysis.

As in [START_REF] Perronnin | Improving the fisher kernel for large-scale image classification[END_REF], we consider the Jacobians with respect to non-scalar parameters only, so that the FV consists of the concatenation of two vectors J Q µ k and J Q σ k . One can easily show (see [START_REF] Sánchez | Image classification with the fisher vector: Theory and practice[END_REF]) that the ((k -1)6 + j)-th element of the above vectors (1 ≤ j ≤ 6, 1 ≤ k ≤ K), i.e. the j-th entry for the k-th mixture component, is given by:

J Q µ k (j) = 1 M √ π k M i=1 γ k,i q j i -µ j k σ j k , J Q σ k (j) = 1 M √ 2π k M i=1 γ k,i q j i -µ j k σ j k 2 -1 , (4)
where γ k,i is the posterior probability that q i belongs to kth cluster conditioned by Q. The normalization by M is added to avoid dependence on the Q's cardinality. Since quads live in R 6 , our Fisher vectors are reasonably long, i.e., of dimension 12K. Typically, d-dimensional descriptors imply 2Kd-dimensional Fisher vectors in a similar framework, a strong disadvantage when long descriptors are employed.

Once the two vectors of (4) are computed and appended in a single vector, a power-normalisation step is applied, i.e., (

) 5
This normalisation step eliminates the sparseness of the Fisher vector, thus increasing its discriminability. Note that the FVs tend to be sparse since the majority of the quads are assigned with high posterior probability to a couple of components.

The impact of such a normalization is evident in Fig. 3 which depicts the distribution of a Fisher vector's elements before and after the power normalisation. We refer the reader to [START_REF] Perronnin | Improving the fisher kernel for large-scale image classification[END_REF], [START_REF] Sánchez | Image classification with the fisher vector: Theory and practice[END_REF] for a detailed discussion about power normalization.

Since any permutation of Q's elements would lead to the same FV, the order of different postures within an action is not taken into account. Therefore, and similar to spatial pyramids in scene recognition [START_REF] Lazebnik | Beyond bags of features: Spatial pyramid matching for recognizing natural scene categories[END_REF], we further enable a multilevel splitting of sequences, so that n video non-overlapping segments are present in the n-th level. The concatenation of FVs of all segments somewhat makes the representation temporally-and order-dependent as long as the order of subactions is encoded in the action descriptor. Three pyramid levels are used in this paper as shown in Fig. 4.

IV. ACTION CLASSIFICATION

It is beyond the scope of this paper to focus on the classification step. We simply employ linear SVMs trained in an one-versus-all fashion in order to build a multi-class classifier, while a cross-validation on training sets provides the best offset per classifier. Notice that a Fisher kernel classifier is equivalent to a linear classifier on FVs [START_REF] Jaakola | Exploiting generative models in discriminative classifiers[END_REF]. Moreover, linear SVMs easily deal with the high-dimensional representations that result from our mutli-level FVs. However, a more in-depth analysis on the classification step will possibly lead to higher performance. To implement the classifier, we make use of the LIBSVM library [START_REF] Chang | Libsvm: a library for support vector machines[END_REF].

V. EXPERIMENTS

In this section, we test our action recognition method on widely used datasets and compare it with the state-of-the-art. Two publicly available datasets are used: MSR-Action3D [START_REF] Li | Action recognition based on a bag of 3d points[END_REF] and HDM05 [START_REF] Müller | Documentation of mocap database hdm05[END_REF]. The datasets are captured from different modalities and regard different human activities.

A. MSR-Action3D Dataset

MSR-Action3D dataset [START_REF] Li | Action recognition based on a bag of 3d points[END_REF] is a set of depth videos, captured by a Kinect device, that contains 20 actions: high arm wave, horizontal arm wave, hammer, hand catch, forward punch, high throw, draw x, draw tick, draw circle, hand clap, two hand wave, sideboxing, bend, forward kick, side kick, jogging, tennis swing, tennis serve, golf swing, pick up & throw. Ten actors perform each action three or two times. Commonly, this dataset is divided into three action subsets as shown in Table I, while the average performance over these subsets is reported as the recognition accuracy [START_REF] Li | Action recognition based on a bag of 3d points[END_REF], [START_REF] Xia | View invariant human action recognition using histograms of 3d joints[END_REF]. Notice that AS1 and AS2 were intended to group similar actions, while AS3 defines a subset with high cross-action variability.

We consider here the challenging case of cross-subject splitting into training and testing sets [START_REF] Xia | View invariant human action recognition using histograms of 3d joints[END_REF], [START_REF] Wang | Mining actionlet ensemble for action recognition with depth cameras[END_REF], [START_REF] Xia | Spatio-temporal depth cuboid similarity feature for activity recognition using depth camera[END_REF], [START_REF] Li | Action recognition based on a bag of 3d points[END_REF], i.e. half of the subjects are used for training and the rest for testing. In particular, the common splitting with odd and even subject-IDs defining the training and testing sets respectively is adopted [START_REF] Li | Action recognition based on a bag of 3d points[END_REF], [START_REF] Wang | Mining actionlet ensemble for action recognition with depth cameras[END_REF]. As in [START_REF] Wang | Mining actionlet ensemble for action recognition with depth cameras[END_REF], we exclude ten sequences where the skeleton data is missing or corrupted, namely, 557 sequences are used in total. A 20-joint skeleton is provided that implies 4845 skeletal quads per frame. The number of GMM components is K = 128, thus leading to 1560-element FVs that are power normalised with a = 0.3. As mentioned, three levels are considered, hence (1 + 2 + 3 =)6 FVs are obtained per sequence.

A fair comparison suggests using competitors that exploit the skeleton structure only. For the sake of completeness, however, we present in a separate table the recognition accuracy of methods that use multiple features and combine joints with other depth cues. Table II shows the recognition accuracy per action subset along with the corresponding results of methods that rely on skeleton joints. Notice that SMIJ [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] and LDS [START_REF] Chaudhry | Bioinspired dynamic 3d discriminative skeletal features for human action recognition[END_REF] methods use a reduced dataset of 17 actions. Unlike most of the competitors, our method provides equally good results in all action subsets, while it provides the best overall accuracy. Fig. 5 depicts the confusion matrices we obtain per action subset. Actions with similar poses, such as the drawing actions in AS2, are more easily confused owing the similarity of the resulting quads. Instead, actions in AS3 are better discriminated, as expected. For comparison reasons, we also show the confusion matrices obtained by the EigenJointbased method [START_REF] Yang | Eigenjoints-based action recognition using naivebayes-nearest-neighbor[END_REF].

Table III shows the recognition accuracy of methods that use the full-body depth map or combine skeleton information with other features. To the best of our knowledge, the method by Ohn-Bar and Trivedi [START_REF] Ohn-Bar | Joint angles similiarities and hog 2 for action recognition[END_REF] performs best when it combines the joints angles with two other features. However, when only joint information is employed, our method performs better (see Table II). It is important to note that our algorithm competes with the majority of these methods, despite the fact that they employ multiple features. As a consequence, the use of skeletal quads in conjunction with other features suggests an even promising approach.

B. Mocap databse HDM05

In this subsection, we present results when applying our method on skeleton sequences of a Motion Capture dataset, the HDM05 database [START_REF] Müller | Documentation of mocap database hdm05[END_REF]. Unlike MSRAction3D, this dataset is captured by motion-capture sensors that acquire more precise data. Moreover, the frame rate is much higher (120 fps), while 31 joints are provided per pose instance. However, we consider here only 15 joints: root, L/Rhip, L/Rknee, L/Rankle, neck, head, L/Rshoulder, L/Relbow, L/Rwrist. As a result, each skeleton pose implies 1365 quads.

We adopt the experimental setup of [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] that suggests a set of 11 actions: deposit floor, elbow to knee, grab high, hop both legs, jog, kick forward, lie down floor, rotate both arms backward, sneak,squat, and throw basketball. The actions are performed by 5 subjects, while each subject performs each action a couple of times (not fixed); this suggests a set of 249 sequences. As with [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF], we use a cross-subject splitting with Table IV shows the recognition accuracy of our method along with several counterparts of the algorithms presented in [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] and [START_REF] Chaudhry | Bioinspired dynamic 3d discriminative skeletal features for human action recognition[END_REF]. Notice that all methods in [START_REF] Chaudhry | Bioinspired dynamic 3d discriminative skeletal features for human action recognition[END_REF] regard global linear dynamical systems (LDS) that describe the dynamics of several joint features across the whole sequence. Then, the LDS parameters encode each action sequence and are classified by SVMs. Moreover, both methods enable a mining step that provides the most discriminative features. Even so, our simple method is more robust and provides higher recognition accuracy. The confusion matrix of our results is shown in Fig. 6. 2 Despite the more accurate joints in this dataset, there is still confusion between some actions, such as deposit floor and kick forward. Note that deposit floor means "deposit an item on the floor with knees bent" [START_REF] Müller | Documentation of mocap database hdm05[END_REF]. As a consequence, there are similar leg poses between these two actions, hence similar quads.

VI. CONCLUSIONS

A local skeletal representation was proposed in this paper. This representation implies a short, view-invariant descriptor of joint quadruples. Furthermore, a Fisher kernel representation was devised that encodes the generation of such a representation from a Gaussian mixture model. The final action descriptor results from a multi-level representation of Fisher 2 confusion matrices are not provided by [START_REF] Ofli | Sequence of the most informative joints (smij): A new representation for human skeletal action recognition[END_REF] and [START_REF] Chaudhry | Bioinspired dynamic 3d discriminative skeletal features for human action recognition[END_REF] Future work regards the combination of skeletal quads with other cues. As well, the dimensionality reduction of the final action descriptor will lead to more efficient classification and will allow the use of non-linear SVMs.

Fig. 1 .

 1 Fig. 1.A Gaussian Mixture Model (GMM), learnt on training data, is supposed to generate skeletal quads. Based on the GMM parameters, the skeletal quads of any action example are encoded into a Fisher vector, thus building the action descriptor which is led to a multi-class linear SVM.

Fig. 2 .

 2 Fig. 2. An example of coding the marked joint quadruple {Neck, Shoulder, Elbow, Hand}. The Neck and Hand joints are identified with the points (0, 0, 0) and (1, 1, 1) in the new local coordinate system. The local 3D coordinates of the Shoulder and Elbow joints describe the structure of the quadruple and the quad descriptor is q = [Sx, Sy, Sz, Ex, Ey, Ez] T .

Fig. 3 .

 3 Fig. 3. The distribution of entries (left) before and (right) after the power normalisation (α = 0.5) of Fisher vector.

Fig. 4 .

 4 Fig. 4. A 3-level representation of an action example. The complete action descriptor is the concatenation of the six FVs.

3 and 2

 2 subjects in training and testing sets respectively, thus having 139 training and 110 testing examples at our disposal. The same parameters with MSR-Action3D are used.

Fig. 5 .

 5 Fig. 5. The confusion matrices obtained by (top) [21] and (bottom) our method for each action subset of MSRAction3D dataset (left: AS1, middle: AS2, right: AS3).

  angles + LDS (Ofli et al. [10]) 76.15% Joint angles + HMW * (Ofli et al. [10]) 78.90% Joint angles + HMIJ * * (Ofli et al. [10]) 82.57% Joint angles + SMIJ (Ofli et al. [10]) 84.47% Joint Shape + LDS (Chaudhry et al. [3]) 82.57% Joint Tangents + LDS (Chaudhry et al. [3]) 88.07% Joint positions + LDS (Chaudhry et al. [3]) 91.74% FV of Skeletal Quads 93.89% * histogram of motion words, * * histogram of most informative joints vectors that encodes the temporal nature of action examples.Experimental validation of the proposed method verified its state-of-the-art performance in human action recognition from depth data.

Fig. 6 .

 6 Fig.6. The confusion matrix of our method on HDM05 dataset.

TABLE I

 I

	.	THE THREE ACTION SUBSETS (AS) OF MSRACTION3D
		DATASET AS DEFINED IN [8]
		AS1	AS2	AS2
	Horizontal arm wave	High arm wave	High throw
		Hammer	Hand catch	Forward kick
		Forward punch	Draw x	Side kick
		High throw	Draw tick	Jogging
		Hand clap	Draw circle	Tennis swing
		Bend	Two hand wave	Tennis serve
		Tennis serve	Forward kick	Golf swing
		Pickup & throw	side boxing	Pickup & throw

TABLE III

 III

	.	RECOGNITION ACCURACY ON MSRACTION3D DATASETS.
	THE METHODS EITHER COMBINE SKELETON JOINTS WITH OTHER DEPTH
		CUES OR USE VARIOUS DEPTH FEATURES
	Methods		Accuracy
	Bag of 3D points (Li et al. [8])	74.70%
	Radnom Occupancy Patterns (Wang et al. [17])	86.50%
	Space-time Occupancy Patterns (Vieira et al. [16])	87.50%
	Joints + Actionlets (Wang et al. [18])	88.20%
	HON4D (Oreifej and Liu [12])	88.89% *
	Joints + Depth Cuboids (Xia and Aggarwal [19])	89.30%
	Joint Angles + MaxMin + HOG 2 (Ohn-Bar and Trivedi [11])	94.84% *
	*performance obtained with different cross-subject splitting (first five actors for training,
	last five actors for testing).

TABLE II .

 II RECOGNITION ACCURACY ON MSRACTION3D DATASET USING SKELETON JOINTS.

		AS1	AS2	AS3	Overall
	Histogram of 3D Joints (Xia et al. [20])	87.98%	85.48%	63.46%	78.97%
	EigenJoints (Yang & Tian [21])	74.50%	76.10%	96.40%	82.33%
	Joint Angles (Ohn-Bar & Trivedi [11])	N/A	N/A	N/A	83.53%
	Joint angles + SIMJ * (Ofli et al. [10])	N/A	N/A	N/A	47.06%
	Joint angles + LDS * (Chaundry et al. [3])	N/A	N/A	N/A	83.89%
	FV of Skeletal Quads	88.39%	86.61%	94.59%	89.86%

* results on a reduced dataset with less actions

The code to extract skeletal quads for both cases of normalization is provided by https://team.inria.fr/perception/research/icpr2014/