Which frog’s legs do froggies eat? The use of DNA barcoding for identification of deep frozen frog legs (Dicroglossidae, Amphibia) commercialized in France

. Several millions frogs captured in the wild in Indonesia are sold for food yearly in French supermarkets, as deep frozen frog legs. They are commercialized as Rana macrodon , but up to 15 look-alike species might also be concerned by this trade. From December 2012 to May 2013, we bought 209 specimens of deep frozen frog legs, and identified them through a barcoding approach based on the 16S gene. Our results show that 206 out of the 209 specimens belong to Fejervarya cancrivora , two to Limnonectes macrodon and one to F. moodiei . Thus only 0.96 % of the frogs were correctly identified. Unless misclassification was intentional, it seems that Indonesian frog leg exporters are not able to discriminate between the species. The quasi absence of L. macrodon in our samples might be an indication of its rarity, confirming that its natural populations are declining rapidly, in agreement with its “vulnerable” status according to the IUCN Red List. Our results show that the genetic and morphological diversity of the frogs in trade is much higher than the genetic and morphological diversity measured so far by scientific studies. These results underline the need for large scale studies to assess the status of wild populations.


Introduction
The international traffic of wild animals is considered an important threat to many animal species and is subject to international regulations.From 1998 to 2007 (10 years), 35 million animals and plants registered on the CITES lists have been exported from South-East Asian countries (Nijman 2010).For the same period, between 180 million and 1 billion specimens of frogs were collected every year in the wild in Indonesia, and one eighth of these frogs were then exported to Europe (Kusrini 2005;Kusrini & Alfold 2006).For the ten-year period from 2000 to 2009, 46 400 tonnes of frog legs were imported into Europe, which corresponds to approximately 928 million to 2.3 billion frogs (Altherr et al. 2011).
The main European countries that import frog legs from Indonesia are France and Belgium (Kusrini & Alfold 2006;Altherr et al. 2011).The consumption of frogs is a tradition in French cuisine, causing the French to be called "froggies" by their neighbours, and in former days private people and restaurants would collect local frog specimens at certain periods of the year.However, in the 1970s deep freezing technology was developed, which allowed long term storage and large scale transport of frog legs, and large numbers of frogs were collected by commercial enterprises (Dubois 1983;Neveu 2004).This modified the conditions of traditional local sustainable collection and drove some local frog populations to extinction.A wildlife protection law was voted in 1979 to protect the French frog species (Le Serrec 1988).The market turned to tropical countries, mainly Bangladesh, India and Indonesia, but India and Bangladesh quickly stopped frog exports: the main species concerned (Hoplobatrachus tigerinus (Daudin, 1802)) was put on CITES lists and very protective laws for the local anuran species were enacted.Thus, in 1987 Indonesia became the main exporter of frogs and more than 80% of European frog leg imports now come from this country (Altherr et al. 2011). From 1973to 1987, 830 to 2659 tonnes of frog legs were imported to France from Indonesia (Le Serrec 1988).This import was 5600 tonnes in 1992 (Kusrini & Alford 2006), 4600 tonnes per year in 2000-2009 (Altherr et al. 2011), and 2906 to 3275 tonnes per year in 2009-2011(French Customs 2012).Frog legs from Indonesia are sold in supermarkets, particularly in those specialised in deep frozen food, or in Asian food markets.Kusirini & Alford (2006) showed that three species are predominantly collected in Indonesia for consumption: the giant Javan frog Limnonectes macrodon (Duméril & Bibron, 1841), the crab-eating frog Fejervarya cancrivora (Gravenhorst, 1829) and the common grass frog Fejervarya limnocharis (Gravenhorst, 1829), the latter probably representing a mixture of F. limnocharis and F. iskandari Veith, Kosuch, Ohler & Dubois, 2001(Djong 2007).About fifteen look-alike species of large size might also be involved in collection for international trade (Altherr et al. 2011;Ohler unpubl. report), and they all belong to the genera Fejervarya Bolkay, 1915 andLimnonectes Fitzinger, 1843.Some of these species are common in the wild (e.g., F. cancrivora) while others are uncommon and listed as vulnerable in the IUCN Red List (e.g., L. macrodon, see Iskandar et al. 2004).Nonetheless, the European regulation no.2065/2001 (Anonymous 2001) asks that the labels of fishery products show the common and scientific names on pre-packed, non-transformed fresh products in order to assure traceability of specimens.In France, deep frozen frog legs are usually sold in the supermarkets in plastic bags of 500 g or 1 kg and bear labels with a Latin species name and their country of origin.
Monitoring and management of harvested populations require accurate species identification; yet, high error rates are likely (Warkentin et al. 2009).Frog species can often be distinguished on various morphological characters such as body proportions, foot morphology and coloration pattern.Because exported frogs are skinned and the bodies cut off before packing, the characters allowing a morphological identification cannot be observed on deep frozen legs and thus identification is difficult (Warkentin et al. 2009).In this paper we explore the use of DNA barcoding as an identification tool to identify deep frozen frog legs commercialized in France.DNA barcoding is a technique that uses a short DNA sequence from a standard locus as a species identification tool (Hebert et al. 2003).
From December 2012 to May 2013, we bought 209 frogs in different supermarkets in France that were labeled as Rana macrodon (a variant of the valid name Limnonectes macrodon).These frogs were sold under three different brands but all come from the same importer.The samples were identified through a barcoding approach based on the 16S gene.The 16S gene was previously recognized as one of the most effective genes for molecular amphibian species identification (Vences et al. 2005;Grosjean et al. 2015).Moreover, 16S reference sequences were available in the GenBank database for the 15 species potentially collected in Indonesia for commercialization.
A diminution in size of adult frogs due to overharvesting has been suggested (Iskandar, cited in Anonymous 2007).As snout-vent length and tibia length are highly correlated in frogs (Emerson 1978;Vidal-García et al. 2014), we used tibia length to obtain an estimate for the body size of commercialized frogs.

Biological samples, DNA extractions, PCR conditions, DNA sequencing
In this study, a total of 209 biological samples of commercialized frogs were bought in different French supermarkets from December 2012 to May 2013 (Table 1, Appendix).All these samples were sold as frozen frog legs in bags of 500 or 1000 g.They came from Indonesia and were sold as Rana macrodon.All these samples are stored in the collections of the National Museum of Natural History, Paris, France (Appendix).
DNA extracts were prepared from muscle tissue using the NucleoSpin Tissue Core kit (Macherey Nagel).16S amplification and sequencing were obtained using the primer pair 16SA-L (5'-CGCCTGTTTATCAAAAACAT-3') and 16SB-H (5'-CCGGTCTGAACTCAGATCACGT-3') (Palumbi et al. 1991).The thermal profile consisted of 38 cycles at 94°C for 45 s, 55°C for 45 s and 72°C for 1 min.The amplicons obtained were about 550 bp long.PCR products were purified and sequenced using Abi technology at the Genoscope (Evry, France).Sequences were checked by eye in CodonCode v. 5.1.The sequences obtained were deposited in the GenBank database under accession numbers KX055940−KX055957.

Molecular species identification of samples
To identify all our samples we first determined the number of haplotypes in our dataset using TCS (Clement et al. 2000).Each haplotype was then assigned to a species through a nucleotide BLAST approach conducted in GenBank (http://blast.ncbi.nlm.nih.gov) and through a phylogenetic method (Bayesian tree).
BLAST species identification was accomplished following three different criteria (Meier et al. 2006 identifications to the closest match regardless of genetic distance.The BCM criterion is similar to BM, but the query is identified by the closest match with a distance below a defined threshold.Finally, the ASB criterion is similar to the BCM by applying a threshold but it returns all the sequences within it.A query is identified when all the matching sequences below the threshold are conspecific.For BCM and ASB, a query may provide ambiguous results if sequence divergences of different species are below the threshold (ASB) or sequences from different species are the closest match below the threshold (BCM).
A tree-based approach of species delimitation was also used.The tree-based criterion of reciprocal monophyly was used to define species boundaries.In our analysis we included: 1) all the haplotypes recovered in our deep frozen frog legs, 2) all species of the genera Fejervarya and Limnonectes known from Indonesia (except Limnonectes dammermani (Mertens, 1929), L. kenepaiensis (Inger, 1966), L. rhacodus (Inger, Boedi & Taufik, 1996), L. sinuatodorsalis Matsui, 2015, L. timorensis (Smith, 1927) and Fejervarya schlueteri (Werner, 1893) for which no 16S gene sequences are currently available in GenBank), and 3) one sequence of Hoplobatrachus rugulosus (Wiegmann, 1834) and one sequence of Occidozyga laevis (Günther, 1858) included as outgroups.Depending on 16S gene sequence availability in GenBank, we included one to four individuals per species in our analysis.Sequences were aligned with ClustalW (Thompson et al. 1994).In order to minimize the number of missing data we kept 553 sites in our final analyses.Evolutionary relationships among sequences were estimated by conducting Bayesian Markov chain Monte Carlo phylogenetic analyses (MCMC) with MrBayes v. 3.1.2(Huelsenbeck & Ronquist 2001).MrModeltest v. 3.04 (Nylander 2004) was used to evaluate the fit of 24 nested models of nucleotide substitution to the data.According to the Akaike information criterion, MrModeltest recommended the GTR + I + G model.This model was used in the Bayesian analysis.
Three heated chains and one single cold chain were employed, and runs were initiated with random trees.Two independent MCMC runs were conducted with five million generations per run; trees (and parameters) were sampled every 100 generations.Stationarity was assessed by examining the average standard deviation of split frequencies and the Potential Scale Reduction Factor (Ronquist et al. 2005).
For each run, the first 25% of sampled trees were discarded as burn-in.

Geographical origin of commercialized frogs
For the species Fejervarya cancrivora, 16S sequences of specimens from different geographical regions are available in GenBank (Table 2).We conducted a phylogeographic analysis to see if there is any phylogeographic signal within this species and whether the geographical origin of the commercialized frogs could thus be determined based on 16S gene sequences.To this aim, 16S gene sequences of F. cancrivora for which the geographical origin was known were retrieved from GenBank.To minimize the number of missing data, 257 sequences of 494 bp were retained in final analyses (206 sequences of commercialized frogs + 51 specimens from GenBank).Relationships among haplotypes were inferred by constructing a network using the median-joining method available in Network v. 5.0.0.0 (Bandelt et al. 1999).

Morphometry
Snout vent length (SVL) and tibia length (TL) of Fejervarya cancrivora were obtained from Boulenger (1920) and from measurements of adult specimens from Borneo, Java and Sumatra stored at the Field Museum of Natural History, Chicago (FMNH 256671, 256688, 256690, 256692-97, 256709-14).Tibia length of frog legs was measured with a slide caliper (Appendix).The Pearson correlation coefficient showed a highly significant correlation between snout vent length and tibia length (r = 0.973; p < 0.001***).For all specimens measured we calculated the ratio r = TL/SVL.We used R = mean of all r to estimate body size depending on TL.We obtained the value R = 0.48177; n = 28; standard deviation = 0.197.Thus, we could use the following formula to calculate an estimate for snout vent length for the frog legs: SVL = TL / 0.48177.
Our nucleotide BLAST analysis showed that haplotypes h01 to h16 correspond to the species Fejervarya cancrivora (best close match: 99 to 100% of sequence identity; Table 1).Three distinct species have been recognized in "Fejervarya cancrivora" (designated as large, mangrove and Sulawesi types; Hasan et al. 2012).The large type of F. cancrivora was designated as the nominal F. cancrivora (Kotaki et al. 2010), while the mangrove and Sulawesi types were designated as F. moodiei (Taylor, 1920) and an undescribed species, respectively (Kurniawan et al. 2011) In our phylogenetic tree (Fig. 1), all Indonesian species of the genera Fejervarya and Limnonectes are reciprocally monophyletic and Bayesian posterior probabilities are high (pp > 0.98) for all species.This tree clearly shows that haplotypes h01 to h16 correspond to the species Fejervarya cancrivora, haplotype 17 corresponds to the species F. moodiei and haplotype h18 corresponds to the species Limnonectes macrodon.

Phylogeography of Fejervarya cancrivora
The results of our network analysis show that there is some genetic variability within F. cancrivora (Fig. 2): 21 haplotypes differing by 1 to 5 mutations were found.Thirteen out of the 16 haplotypes identified in the commercialized frogs (representing 95 individuals) have not been recorded in previous analyses.Haplotype h02, found in 54 commercialized frogs, was recovered in Taiwan and several Indonesian islands (Kalimantan, Sumatra, Bali and Java).Haplotype h11 (49 commercialized frogs) was recovered in Malaysia and several Indonesian islands (Bangka, Sumatra, Java).Haplotype h10 (8 commercialized frogs) is present in Taiwan.Due to the lack of phylogeographic structure within previously sequenced F. cancrivora specimens and to the high number of new haplotypes detected in commercialized frogs, it is not possible to infer the geographical origin of the frogs sold in France.

Estimation of snout-vent length of frog legs in trade in France
For 192 frog legs, the estimated snout-vent length ranged from 59.4 to 111.7 mm (mean value: 79.0; standard deviation: 11.37).The histogram (Fig. 3, Appendix) of these measurements shows a bimodal distribution which does not correspond to the size distribution of Fejervarya cancrivora specimens collected in natural populations without sampling bias.This is also reflected by the relatively low standard deviation.Measurements of specimens from Borneo, Java and Sumatra published in Boulenger (1920) vary from 54 to 88 mm (mean value: 72.8; standard deviation: 13.34), whereas specimens collected in the early 1990s range from 49.7 to 101.6 mm (mean value: 78.9; standard deviation: 14.98).

Misidentification of frogs commercialized in France
Results based on the phylogentic tree and the BLAST approach (BM, BCM and ASB) are congruent and show (threshold of 5%; Kurabayashi et al. 2005) that 206 out of the 209 analyzed specimens belong to the species Fejervarya cancrivora, two to the species Limnonectes macrodon and one to the species F. moodiei.Thus only 0.96% of the frogs were correctly identified.
In a previous study on frog legs imported to Belgium as F. cancrivora, L. macrodon, L. limnocharis and Rana catesbeiana Shaw, 1802, all samples proved to be Fejervarya cancrivora.In this case, 34.5% of the identifications were correct (Veith et al. 2000).In that study 36.8% of the frog legs were sold as Limnonectes macrodon.Although all our samples were labelled as Rana macrodon, only 2 frog legs (0.96%) could be identified as this species.As proposed by several authors (Kusrini 2005;Kusrini & Alford 2006;Veith et al. 2000), it seems likely that misclassification is not intentional but due to the fact that frog leg exporters are not able to discriminate between the species in trade.Managers of export companies stipulate that they should be supplied only with frog legs of L. macrodon (Kusrini & Alford 2006).On the other hand, in the local markets L. macrodon and F. cancrivora are correctly identified and sold at different prices, because the meat of L. macrodon is considered of better taste (Kusrini & Alford 2006;A. Ohler, pers. obs.).

Sustainability and conservation
The species concerned by international trade are of similar, large size but show differences in breeding biology and habitat, although relevant data are scarce.Fejervarya cancrivora inhabits marshes and paddy fields, not avoiding habitat modified by man and thus a large area of potential habitat is available (Inger 1966;Iskandar 1998;Yuan et al. 2004).Virtually nothing is known of the habitat and ecology of F. moodiei (Ohler 2004).Fejervarya cancrivora lays a relatively large number (up to more than 2500) of small-sized eggs in successive clutches (Inger 1956;Alcala 1962), which develop in lotic habitats.Limnonectes macrodon is present in riparian secondary forests (Inger 1966).It can be observed in clearings and secondary growth, and it is very rare in primary forests.A single clutch of about 1000 eggs is laid in side pools of rivers (Iskandar 1998).Overharvesting should thus have a higher impact on L. macrodon, as breeding capacity is smaller in L. macrodon than in F. cancrivora and riverside habitats are scarcer than ponds and paddy fields.The conservation status of F. cancrivora was evaluated as "Least Concern" "in view of its wide distribution, tolerance of a broad range of habitats, presumed large population, and because it is unlikely to be declining to qualify for listing in a more threatened category" (Yuan et al. 2004).Fejervarya moodiei is listed as "Data deficient in view of continuing doubts as to its extent of occurrence, status and ecological requirements" (Ohler 2004).Limnonectes macrodon is listed as "Vulnerable because it depends on streams in lowland forest, and so its Area of Occupancy is probably less than 2000 km 2 , its distribution is severely fragmented, and there is continuing decline in the extent and quality of its habitat, the number of locations, and the number of mature individuals" (Iskandar et al. 2004).These evaluations are not based on population studies that might measure the impact of the heavy harvesting ongoing for decades.However, the quasi absence of L. macrodon in our samples might be an indication of its rarity in the field and the fact that its natural populations are declining rapidly.
In the 1990s, researchers already considered that large-sized frogs had completely disappeared in many parts of Java and Sumatra (Manthey & Grossmann 1997;Inger & Stuebing 1997).In 2006, 18.8% of the frogs captured by harvesters in West Java and East Java belonged to the species L. macrodon (Kuzrini & Alford 2006).In 2012-2013, only 0.96% of the samples from trade sold in France belonged to this species.
A diminution in size of adult frogs due to overharvesting has been suggested (Iskandar, cited in Anonymous 2007).The estimated size of the frogs from deep frozen specimens of F. cancrivora is not smaller than specimens collected 20 or a hundred years ago.The frogs measured in Boulenger (1920) are slightly smaller in size than the frogs collected in the early 1990s and in 2012-2013.However, the body size of the frogs in trade is clearly biased.After collection, frogs are separated into size categories and only large specimens are chosen for export (Kusrini & Alford 2006).In samples studied by Church (1960) and Kusrini & Alford (2006) that had been collected for trade, very large-sized frogs were present (maximum size of 132 mm for 1325 frogs, and of 162 mm for 555 frogs, respectively).In the samples bought in the supermarkets in France imported by the same company, two size classes of frogs are present, but very large frogs have not been observed in our sample.The absence of small-sized frogs is due to the sorting out of larger specimens for export (Kusrini & Alford 2006).Nevertheless, the presence of numerous large frogs means that such large specimens can be found in the field.Based on our genetic results, we could not allocate our samples to precise geographic origins.The absence of an effect of overharvesting on body size can only be studied when the precise geographic origin of specimens is known, because body size varies geographically (Kurniawan et al. 2011).The specimens studied in the present work may well come from newly harvested populations.In fact, the presence of F. moodiei in our samples, a species which has a known range outside Java and Sumatra, a range that includes the Philippines, the northern coasts of the Gulf of Thailand, Bangladesh and Orissa (Kurniawan et al. 2010), indicates that the specimens studied here may have come from other islands than Java and south-eastern Sumatra.
"At present, nothing is known about the impact of extensive frog harvesting on Indonesian frog populations and agricultures."This phrase of Veith et al. (2000) can be cited as such, as there has been no change in the last 15 years.There are still no studies on the effects of frog harvesting on natural populations.Our results show that the genetic and morphological diversity of the frogs in trade is much higher than the genetic and morphological diversity measured so far by scientific studies (Kurniawan et al. 2010(Kurniawan et al. , 2011)).These results underline the need for large-scale studies on the taxonomy, population structure, reproductive data and ecology of the species concerned in international trade.In 2000, Veith et al. concluded that the development of a "quick and cheap test for management authorities" was necessary to monitor the international trade in frog legs.Such a test is already available through DNA sequencing, and the price is decreasing rapidly with the development of next generation sequencing techniques.The lack of any progress in the development of identification tools and assessing the status of wild populations is not due to an absence of convenient methodology but it reflects the fact that conservation is focused on flagship species, such as tigers and pandas, or species that qualify for rarity.
The edible frogs of Indonesia, although billions of individuals are killed annually, are the scarcely visible part of the bulk of species that are disappearing.

Appendix
List of samples of deep frozen frog legs with purchase date, collection number, haplotype number, taxonomic identification, tibia length (TL) and estimated snout vent length (SVL).

Haplotype number
Species identification (Blast)

Fig. 1 .
Fig. 1.Phylogeny of Indonesian species of Fejervarya and Limnonectes recovered by the Bayesian analysis (GTR + I + G model).Hoplobatrachus rugulosus (Wiegmann, 1834) and Occidozyga laevis (Günther, 1858) were used as outgroups.Numbers on nodes represent Bayesian posterior probabilities, * indicates a value higher than 0.98.Only values higher than 0.75 are represented.h01 to h18 indicate the 18 haplotypes from frozen frog legs recovered in this study.

Fig. 2 .
Fig. 2. Minimum spanning network depicting relationships among 16S haplotypes of Fejervarya cancrivora(Gravenhorst, 1829).The size of each circle is proportional to the haplotype frequency and the lengths of the connecting lines are proportional to the number of mutations.Colors refer to distinct regions (Indonesia: Java, Sumatra, Bali, Kalimantan, Bangka; Malaysia; Taiwan) and commercialized frogs of unknown origin are in black.

Fig. 3 .
Fig. 3. Histograms. A. Snout vent length (in mm) in adult Fejervarya cancrivora (Gravenhorst, 1829) from samples collected for scientific purposes (Boulenger 1920) and collection specimens as mentioned in Material and methods.B. Snout vent length estimated from tibia length of genetically identified frog legs from French supermarkets (specimen list, see Appendix).

Table 1 .
Species identification recovered from a Blast in GenBank (16S gene) of the 209 specimens of frogs bought in French supermarkets from December 2012 to May 2013 and identified on the label as "Rana macrodon from Indonesia".

Table 2 .
List of specimens of Fejervarya cancrivora used in our phylogeographic analysis.