Barcoding Analysis of Paraguayan Squamata

: Paraguay is a key spot in the central region of South America where several ecoregions converge. Its fauna (and speciﬁcally its herpetofauna) is getting better studied than years before, but still there is a lack of information regarding molecular genetics, and barcoding analyses have proven to be an excellent tool in this matter. Here, we present results of a barcoding analysis based on 16S rRNA gene sequences, providing valuable data for the scientiﬁc community in the region. We based our ﬁeldwork in several areas of Paraguay. We analyzed 249 samples (142 sequenced by us) with a ﬁnal alignment of 615 bp length. We identiﬁed some taxonomic incongruences that can be addressed based on our results. Furthermore, we identify groups, where collecting e ﬀ orts and research activities should be reinforced. Even though we have some blanks in the geographical coverage of our analysis—and there is still a lot to do towards a better understanding of the taxonomy of the Paraguayan herpetofauna—here, we present the largest genetic dataset for the mitochondrial DNA gene 16S of reptiles (particularly, Squamata) from Paraguay, which can be used to solve taxonomic problems in the region.


Introduction
Ideally, profound knowledge of biodiversity is the first step before any conservation action, sustainable management or biological study is carried out in a given area [1,2]. The scientific community is not only aware that the world is facing a major extinction event [3], but also that many lineages are disappearing even before becoming known to science [4,5]. In the last decades, the use of molecular data has helped to improve our knowledge about biological diversity and the use of genetics as a tool for species recognition is now routine. Molecular data are more often used every day and are applied to the species identification [6][7][8] and species delimitation [9][10][11][12] of all kinds of living organisms, with applications even in food control quality [13].
In this context DNA barcoding analysis is the examination and comparison of short and stable fragments of DNA (DNA barcodes), usually mitochondrial, that represent genetic identifiers for a species [14], and have proven to be a reliable technique for taxonomy [15][16][17]. Nevertheless, researchers have to be cautious because despite being a power tool, barcoding analysis potentially can lead to misinterpretations if sequences used for comparison were generated from misidentified

Data Collection
Even though there are some blanks in the areas sampled, the coverage of collecting sites in this study is rather vast ( Figure 2). Nevertheless, there are two ecoregions in the Occidental Region (Cerrado Chaqueño and Médanos del Chaco), from where we have no samples. The methods used The topography of the Chaco region is a flat savanna with few small isolate hills in the center/north and east. Bad drainage creates vast flooded areas especially south and east, north-west has some dunes formation and is extremely dry. The oriental region is undulated with hills and the highest point in the center is well irrigated by tributaries of the Paraguayan and Parana rivers basins.
The climatic conditions vary in a northwestern-southeastern gradient, being more humid and cooler in the southeast. The mean annual temperature in the whole country is about 23 • C, being 24.5 • C in the western region and 22.5 • C in the eastern region. It is important to note that there are two seasons, the wet season, in which it rains frequently, coincide with the warm period from October to April and the dry season from May to September, where rain is less frequent and coincides with the coldest period.
There is a big difference with respect to the variation in temperature, given that the mean maximum is 25 • C, but the absolute maximum temperature could reach around 50 • C, especially in the northwestern portion of the country in January or February. The coldest month is July, and the absolute minimum temperature can be −6 • C in the south. Nevertheless, in Paraguay the "true" winter usually does not last longer than 16 days each year. Thus, Paraguay has a warm/hot climate during most parts of the year.

Data Collection
Even though there are some blanks in the areas sampled, the coverage of collecting sites in this study is rather vast (Figure 2). Nevertheless, there are two ecoregions in the Occidental Region (Cerrado Chaqueño and Médanos del Chaco), from where we have no samples. The methods used in the field were the traditional techniques for herpetology: Active searching at different times of the day and night, examining potential shelters (e.g., barks, logs, caves, mud, leaf litter, etc.) [46] (Figure 3). Fast moving lizards (e.g., Ameiva and Teius) were collected using compressed air rifles [47]. Additionally, some habitats, such as ant nests and swamps, were dug looking for hypogeal organisms [48], and floating vegetation was sampled using a trawl net (Figure 3). In total, 147 days of fieldwork were accounted for this project, and about 400 specimens collected.
It is important to highlight that the exotic lizard Hemidactylus mabouia is currently widely distributed in the country and now is part of the Paraguayan herpetofauna, and thus also included in the study. These genetic data may help in future studies about the colonization of the species, which in Paraguay has been recorded in the Concepción, San Pedro, Central, Alto Paraná, and Itapúa departments [49].
Reptiles that were captured alive were euthanized with a pericardial injection of a solution of embutramide, mebezonium iodide, and tetracaine hydrochloride (T-61 ® , Intervet International GmbH, Unterschleissheim, Germany) or Sodium Thiopental (Tiopental Sódico ® , Biosano, Chile in the field were the traditional techniques for herpetology: Active searching at different times of the day and night, examining potential shelters (e.g., barks, logs, caves, mud, leaf litter, etc.) [46] ( Figure  3). Fast moving lizards (e.g.: Ameiva and Teius) were collected using compressed air rifles [47]. Additionally, some habitats, such as ant nests and swamps, were dug looking for hypogeal organisms [48], and floating vegetation was sampled using a trawl net ( Figure 3). In total, 147 days of fieldwork were accounted for this project, and about 400 specimens collected.
It is important to highlight that the exotic lizard Hemidactylus mabouia is currently widely distributed in the country and now is part of the Paraguayan herpetofauna, and thus also included in the study. These genetic data may help in future studies about the colonization of the species, which in Paraguay has been recorded in the Concepción, San Pedro, Central, Alto Paraná, and Itapúa departments [49].    Reptiles that were captured alive were euthanized with a pericardial injection of a solution of embutramide, mebezonium iodide, and tetracaine hydrochloride (T-61® , Intervet International GmbH, Unterschleissheim, Germany) or Sodium Thiopental (Tiopental Sódico® , Biosano, Chile). The Secretaría del Ambiente from Paraguay (Currently "Ministerio del Ambiente y Desarrollo Sustentable") authorized the collecting of specimens though permits SEAM [Secretaría del Ambiente] N° 004/11 and 009/2014. Exportation permits for tissues and specimens were also issued by the same authority through the permits SEAM N° 002/14, 016/2016, and 084/2016. After euthanasia, tissue samples were taken either from the muscle of the thigh, tongue, finger clips, tail (when regenerated), or liver. Tissues were preserved in vials containing 98% non-denatured ethanol, and stored at −20 °C as soon as possible.
Hemipenes of Squamata were everted after euthanasia, with an injection of 70% ethanol after manually everting the organs. All specimens were fixed with a solution of 36% formalin and 96% ethanol in the proportion of 5:1000 (e.g., 5 mL formalin in 1 l ethanol), injected in the body cavity, thighs, and thickest part of the tail. Following fixation, the specimens were maintained in 70% ethanol.

Molecular Protocols
We used two different methods of DNA extraction. For sets containing few samples (usually eight or fewer), we used the DNeasy® Blood & Tissue Kit of Qiagen® (Hilden, Germany), whereas for sets of 96 samples we used the fiberglass plate [50]. Both methods are detailed below. The DNA was isolated from tissues whenever possible, or taken from preserved specimens that had been stored for a considerable time in 70% ethanol at room temperature in some cases.
For the DNeasy® method, we used tissue fragments of ~2 mm 2 . When buffers formed precipitates, they were warmed up at 56 °C before use. All reagents for this protocol are included in the kit. Tissues were digested adding 180 μL (all values are for individual samples) of ATL Buffer and 20 μL of proteinase K. Samples in that mix were incubated in a rocking platform at 56 °C for 4 to 12 h until the tissue was completely lysed. Sampling methods during fieldwork included diverse techniques to search in different environments.
After euthanasia, tissue samples were taken either from the muscle of the thigh, tongue, finger clips, tail (when regenerated), or liver. Tissues were preserved in vials containing 98% non-denatured ethanol, and stored at −20 • C as soon as possible.
Hemipenes of Squamata were everted after euthanasia, with an injection of 70% ethanol after manually everting the organs. All specimens were fixed with a solution of 36% formalin and 96% ethanol in the proportion of 5:1000 (e.g., 5 mL formalin in 1 L ethanol), injected in the body cavity, thighs, and thickest part of the tail. Following fixation, the specimens were maintained in 70% ethanol.

Molecular Protocols
We used two different methods of DNA extraction. For sets containing few samples (usually eight or fewer), we used the DNeasy ® Blood & Tissue Kit of Qiagen ® (Hilden, Germany), whereas for sets of 96 samples we used the fiberglass plate [50]. Both methods are detailed below. The DNA was isolated from tissues whenever possible, or taken from preserved specimens that had been stored for a considerable time in 70% ethanol at room temperature in some cases.
For the DNeasy ® method, we used tissue fragments of~2 mm 2 . When buffers formed precipitates, they were warmed up at 56 • C before use. All reagents for this protocol are included in the kit. Tissues were digested adding 180 µL (all values are for individual samples) of ATL Buffer and 20 µL of proteinase K. Samples in that mix were incubated in a rocking platform at 56 • C for 4 to 12 h until the tissue was completely lysed.
Following digestion, 200 µL of AL lysis Buffer + 200 µL of ethanol (98%) were added. This mix was centrifuged (8000 rpm) in DNeasy ® Mini spin columns, discarding all the flow-through. Then, 500 µL of AW1 washing Buffer was added and centrifuged (8000 rpm) discarding the flow-through. Finally, 500 µL of AW2 washing Buffer was added and centrifuged (14,000 rpm) discarding the flow-through. The final elution was made with 200 µL of AE Buffer, after an incubation of one minute, followed by centrifugation (8000 rpm).
For the fiberglass extraction method, we used tissue fragments of~1 mm 2 . Specifications of reagents used in this protocol, are detailed in Table S1. Initially, the samples were washed with 50 µL (values per sample) of a solution of 1× Tris-Ethylenediaminetetraacetic acid (TE) Buffer to remove the remaining ethanol, for~15 h. Following, the samples were digested with 50 µL of a solution of Vertebrate lysis Buffer and proteinase K (10:1), and incubated in a rocking platform at 56 • C for 12-24 h.
Once the samples were digested, the DNA extraction was made adding 100 µL of Binding Buffer and centrifuging at 2800 rpm. These products were transferred to a Pall ® (Cortland, NY, USA) AcroPrep ® filter plate, where the plate was vacuumed for 2 min. Then it was added 180 µL of Washing Buffer 1 and vacuumed again for 2 min. Posteriorly, it was added 750 µL of the Washing Buffer 2 and vacuumed for 2. Then TE Buffer was used to elute the DNA, adding 50 µL and incubating it for 2 min at 56 • C.
We amplified fragments of the mtDNA 16S gene with forward and reverse reactions using the following primers respectively: F: L2510 (5 -CGCCTGTTTAACAAAAACAT-3 ) and R: H3056 Additionally to our own samples, we included non-Paraguayan data from sequences downloaded from GenBank, selecting preferably those sequences associated with museum vouchers, to avoid common problems of misidentifications in that repository [52][53][54]. In most cases, we downloaded only sequences from species represented in our samples, except for Bothrops and Amphisbaena. In these two cases, we downloaded samples from all the species present in Paraguay, because of the difficulty of these taxa for morphological identification. In the case of Micrurus, there is only one sequence of a species present in Paraguay available in GenBank (JQ627286). A particular case was the only available sample of Vanzosaura rubricauda (AF420716) uploaded in the framework of a lizards' phylogeny [55]. That specimen (MRT 05059) from Vacaria, Estado de Bahia, Brazil, actually is V. multiscutata [35]. Nevertheless, it was included in the analysis to evaluate the clustering with the genetic sample from Paraguay.
Codes of sequences downloaded from GenBank, plus accession numbers of sequences generated in this work, are available in the Table S2.

Data Analysis
Chromatograms of forward and reverse sequences were assessed, and a consensus sequence for each sample generated in SeqTrace 0.9.0 [56]. For sequences alignment, we employed MAFFT2 [57,58] through the webserver [59], which includes a special search strategy (Q-INS-i) for the secondary structure of the rRNA 16S [60]. No later manual edition was introduced. Results of MAFFT2 were visualized in MSA Viewer [61] and exported as fasta files.
The best scheme for substitution model was explored in PartitionFinder 2.1.1 [62], using linked branch lengths (supported by most of the phylogenetic programs) using a PhyML 3.0 analysis [63]. Given that the analysis is based on Squamata, which can show highly divergent clades, we estimated the relative quality of the statistical models using the Bayesian Information Criterion (BIC) [64] since it penalizes more the number of parameters in the model and then is better for a large degree of heterogeneity [65].
Given that it is not recommended to use both +I (significant proportion of invariable sites) plus +G (rate of variation among sites follows a gamma distribution) together in the same substitution model, we chose the best suggested model using +I or +G in the partition schemes, but never both together [66].
The phylogenetic hypothesis was performed under a Maximum Likelihood (ML) approach, using IQ-Tree [67] through its webserver [68], setting 10,000 non-parametric bootstrap replicates plus 10,000 replicates of Shimodaira-Hasegawa approximate likelihood ratio (SH-aLRT) [69] and 10,000 ultrafast bootstrap (UFBoot) approximation replicates [70]. We used a sequence of Sphenodon punctatus to root the tree [71], which has been proposed as the sister clade to Squamata [72]. Here, it is important to highlight that the phylogenetic hypothesis is used to sorting groups, and we do not seek for a comprehensive evolutionary reconstruction.
For visualization and edition (branch arrangement, colors, font sizes, etc.) of the tree generated through ML analysis, we used FigTree 1. 4.3 [73]. The final alignment plus the ML tree are stored in TreeBASE repository (https://treebase.org/) under the submission number 24616. To do this, we first managed the alignments and trees in nexus format and combined them in a single file (containing one alignment and the ML tree) using Mesquite 3.31 [74].

Results
We generated a total of 142 sequences of 64 species of Squamata from Paraguay, including one exotic species: Hemidactylus mabouia. In the Table S2, we present a list of specimens used for genetic analyses based on the field work for this project. For comparison we added 107 sequences from GenBank (Table S2). The final alignment constituted of a dataset of 249 samples of 615 bp length. Sequences are available in GenBank.
The best substitution model for the Barcoding dataset was GTR+G, according to the BIC. The sample of Sphenodon punctatus was retrieved as the sister clade to the Squamata (Figure 4). Deep nodes have low bootstrap values, meaning that the phylogenetic relationships are weakly supported. Nevertheless, the shallowest divergences have higher support values, recovering most of the genera included in the analysis as monophyletic, with the exception of Manciola (Scincidae) and the tribe Xenodontini (Colubridae).
The tribe Xenodontini ( Figure S1) of the Subfamily Dipsadinae (Colubridae) contains the samples of Erythrolamprus aesculapii in a monophyletic clade, whereas E. poecilogyrus appears as paraphyletic. Erythrolamprus reginae clusters sister to the above-mentioned taxa. The genus Xenodon seems to be paraphyletic, given that two samples of Xenodon pulcher are sister to Erythrolamprus, whereas Xenodon merremi is sister to Xenodon pulcher + Erythrolamprus. Finally, in this clade a sample of Erythrolamprus sagittifer is sister to a sample of Lygophis dilepis.
Sister to Xenodintini is a clade composed of Phalotris + Philodryas ( Figure S2). Both genera are monophyletic in the tree. The genera Psomophis and Dipsas are clustered together, and nested as sister to the above-mentioned snakes ( Figure S3). Other genera of Dipsadinae that are rendered as monophyletic are Hydrodynastes, Helicops, and Thamnodynastes ( Figure S4). A clade grouping members of the Colubrinae subfamily is composed of three genera (Chironius, Leptophis, and Palusophis) that also show monophyly ( Figure S5). The Pseudoboini (Dipsadinae: Colubridae) is shown in its own clade ( Figure S6) with four of the five genera used in the analysis being monophyletic, whereas the two species of the genus Phimophis appear in different positions of the gene tree.
The genus Micrurus (Elapidae) is the sister clade of the Colubridae, whereas the genus Bothrops (Viperidae) seems to be the sister to Elapidae + Colubridae ( Figure S7). Located in a most basal position among snakes are the two species of Epicrates (Boidae), with Amerotyphlops (Typhlopidae) as the sister clade of the remaining snakes ( Figure S8). The genus Amphisbaena (Amphisbaenidae) is monophyletic, where A. alba and A. bolivica are in their own clades, and A. mertensii shows also monophyly ( Figure S9). Amphisbaena angustifrons is the sister taxon of the other Amphisbaena, with a sample of Amphisbaena sp. (PCS 314) as the most basal taxon of the clade ( Figure S9). In our analysis, Amphisbaena is sister to Teiidae + Gymnophthalmidae. Gymnophthalmidae appears as a monophyletic clade, and the four genera show monophyly as well ( Figure S10).  The tribe Xenodontini ( Figure S1) of the Subfamily Dipsadinae (Colubridae) contains the samples of Erythrolamprus aesculapii in a monophyletic clade, whereas E. poecilogyrus appears as paraphyletic. Erythrolamprus reginae clusters sister to the above-mentioned taxa. The genus Xenodon seems to be paraphyletic, given that two samples of Xenodon pulcher are sister to Erythrolamprus, whereas Xenodon merremi is sister to Xenodon pulcher + Erythrolamprus. Finally, in this clade a sample of Erythrolamprus sagittifer is sister to a sample of Lygophis dilepis.
Sister to Xenodintini is a clade composed of Phalotris + Philodryas ( Figure S2). Both genera are monophyletic in the tree. The genera Psomophis and Dipsas are clustered together, and nested as sister to the above-mentioned snakes ( Figure S3). Other genera of Dipsadinae that are rendered as monophyletic are Hydrodynastes, Helicops, and Thamnodynastes ( Figure S4). A clade grouping members of the Colubrinae subfamily is composed of three genera (Chironius, Leptophis, and Palusophis) that also show monophyly ( Figure S5). The Pseudoboini (Dipsadinae: Colubridae) is shown in its own clade ( Figure S6) with four of the five genera used in the analysis being monophyletic, whereas the two species of the genus Phimophis appear in different positions of the gene tree.
The genus Micrurus (Elapidae) is the sister clade of the Colubridae, whereas the genus Bothrops (Viperidae) seems to be the sister to Elapidae + Colubridae ( Figure S7). Located in a most basal position among snakes are the two species of Epicrates (Boidae), with Amerotyphlops (Typhlopidae) as the sister clade of the remaining snakes ( Figure S8). The genus Amphisbaena (Amphisbaenidae) is monophyletic, where A. alba and A. bolivica are in their own clades, and A. mertensii shows also monophyly ( Figure S9). Amphisbaena angustifrons is the sister taxon of the other Amphisbaena, with a sample of Amphisbaena sp. (PCS 314) as the most basal taxon of the clade ( Figure S9). In our analysis, Amphisbaena is sister to Teiidae + Gymnophthalmidae. Gymnophthalmidae appears as a monophyletic clade, and the four genera show monophyly as well ( Figure S10).
The Family Teiidae is shown as paraphyletic. The only Tupinambinae in our samples was Salvator, which clusters as sister to Gymnophthalmidae ( Figure S11). Samples of Teiinae are clustered The Family Teiidae is shown as paraphyletic. The only Tupinambinae in our samples was Salvator, which clusters as sister to Gymnophthalmidae ( Figure S11). Samples of Teiinae are clustered together showing monophyly, where Ameivula and Kentropyx are sister clades ( Figure S11), as are Teius and Ameiva ( Figure S12). The Family Scincidae is sister to the all above-mentioned clades ( Figure S13). Samples of Manciola show paraphyly ( Figure S13). The remaining cluster contains members of the Anguidae, Gekkonidae, Phyllodactylidae, Liolaemidae, Polychrotidae, and Tropiduridae families. The clade composed by geckos shows monophyly in the genera, but not in the families given that Phyllopezus and Homonota are currently placed in Phyllodactylidae, whereas Hemidactylus and Lygodactylus are Gekkonidae ( Figure S14). Sister to the Gekkota (Gekkonidae + Phyllodactylidae) is Liolaemus, and Stenocercus is rendered as sister to the Gekkota + Liolaemus.

Discussion
Molecular genetics, and in particular barcoding analyses, proved to be a powerful tool to generate preliminary information about the taxonomic status of problematic taxa [20,75]. We present here the most comprehensive analysis of genetic samples of Squamata from Paraguay. The results obtained here will be useful to help identify questionable specimens and in some cases also to clarify some taxonomic issues of the Squamata fauna from the central region of South America. Thus, the data generated here will have a positive impact in a larger geographic context, beyond Paraguay's borders. As said before, genetics alone will not yield a well-founded taxonomy. Nevertheless, molecular genetics open a path for defining operational taxonomic units (OTUs), identifying potential undescribed species and pointing to taxonomic problems, and thus have to be seen as a first informative step and a complementary evidence line in the framework of the modern integrative taxonomic approach [20,76].
Some taxonomic results of this project were already published. For instance, the samples of Colobosaura exhibit large genetic distances, and then Colobosaura kraepelini was revalidated [41]. The Tropidurus samples show monophyly in the species of the torquatus group (T. catalanensis and T. etheridgei), but indicate several uncertainties within the spinulosus group (formerly T. guarani, T. lagunablanca, T. spinulosus, T. tarara, and T. teyumirim), that resulted in the synonymization of T. guarani with T. spinulosus, and T. tarara and T. teyumirim with T. lagunablanca [43]. Regarding the Family Phyllodactylidae, there is strong evidence for the recognition of two different Homonota species in the Chaco [42,44] and a highly distinctive Phyllopezus clade, separated from populations from Cerrado and Chaco [45].
The samples of Vanzosaura rubricauda from Cerrado (field number "ALA") show a high branch distance compared with Vanzosaura rubricauda from Chaco (GK 3801) which is even larger than the distance from V. multiscutata ( Figure S10). Integrating molecular and morphological data, a new species of Vanzosaura (V. savanicola) was previously described, and Gymnodactylus multiscutatus was transferred to the genus Vanzosaura [35]. Nevertheless, their genetic tree [35] included only a single sample from Paraguay and none from Argentina. In their map, obviously two divergent populations of Vanzosaura rubriauda are recognized: One west of the Paraguay River in the Dry Chaco, and another east of the Paraguay River in the Cerrado. Keeping a conservative approach, the authors maintained V. rubricauda as a single taxonomic unit, but with our additional samples it might be possible to generate new taxonomic hypotheses.
Furthermore, we recommend further studies on Amphisbaenidae, because one of the major and latest revisions of Amphisbaenidae in the Neotropics concluded that A. mertensi and A. cunhai (not recorded in Paraguay) are the most basal lineages of the genus [77]. Our analysis showed that the most basal sample (Amphisbaena sp. PCS 314) seems to be a different species as those within the remaining clade. Additional analyses, including more samples and a detailed morphological revision, are necessary to assess the specific status of that specimen.
The weakest part of this work was the analysis of snakes. These animals are usually the harder ones to sample. Compared to the actual diversity of Colubridae, our dataset had fewer samples of this family and therefore it was not possible to draw detailed taxonomic conclusions. However, the presence of the genus Xenodon in two different clusters suggests that more taxonomic work with this group of snakes is needed. Several taxonomic modifications occurred within the Colubridae in the last decade, where the genera Lystrophis and Waglerophis were synonymized with Xenodon, based on the analysis of gene sequences 12S and 16S for the genus Lystrophis, and Cytb and bdnf for one sample of Waglerophis merremii [78]. In our analysis, we found the samples of X. pulcher (previously Lystrophis pulcher) separated from X. merremii (previously Waglerophis merremii). It is desirable to perform phylogenies in this group using more nuclear data to get more robust relationships in the deep nodes.
In the clade of the genus Thamnodynastes, the two species used in our analysis (T. chaquensis and T. hypoconia) are nested in the same node. A more specific genetic study of the genus is highly recommendable. Furthermore, the revision of the phylogenetic status of Phimophis is advisable, since here it appears polyphyletic. In a former study, two samples of Phimophis were used: P. guerini (GQ457761) and P. iglesiasi (JQ598891) and due to polyphyly the authors described the genus Rodriguesophis to include the latter species [78]. This genus is characterized by the absence of the loreal scale. Both P. guerini and P. vittatus have a loreal scale, so they cannot be assigned to Rodriguesophis. Thus, a deeper integrative (morphological and molecular) analysis is needed to understand their relationships.
The genus Micrurus is scarcely represented in GenBank, and comparisons are not possible. The only sample from GenBank is M. altirostris, which is differentiated from Paraguayan samples by a rather long branch distance. There is a polytomy with three samples (PCS 310, 334, and 337), and the only identified specimen is Micrurus pyrrhocryptus (PCS 310) from Pantanal (northern part of Paraguay). The other samples are from Concepción at the other side of the river, and with a body color different from the pattern of M. pyrrhocryptus. It is important to highlight here that some of the specimens that we have of Micrurus, were decapitated (killed by rural farmers) and therefore, without genetic samples for comparison (in GenBank) and without cephalic data (which contains important diagnostic characters) [79], its specific taxonomic allocation becomes difficult.
Some of our Paraguayan samples of Bothrops of the neuwiedi complex, from distant parts of Paraguay, are clustered with a sample of B. diporus from GenBank (Samples PCS 302, PCS 318, PCS 331, PCS 504, Figure S7). A thorough analysis of this complex of Bothrops is needed to understand the true diversity in the group.
Regarding the Scincidae, a sample from GenBank of Manciola guaporicola (KX364960) from Mbaracayú Reserve (Paraguay) appears out of the clade of the remaining M. guaporicola from Paraguay and Brazil. This is also a topic that should be further investigated.
Regarding conservation, one of the major problems in Paraguay for several years was habitat loss due to extensive soybean crops in the eastern part of Paraguay [80]. Nevertheless, habitat fragmentation is currently also affecting the landscapes of the Occidental Region of Paraguay [81,82]. Thus, currently, the protected areas are the best strategy for conservation of biodiversity in Paraguay, although many conservation units face legal problems (e.g., lack of official measurements, management plans, forest guards, infrastructure, etc.) and then the long-term maintenance of their biodiversity is not guaranteed [80]. There are some reptile species absent from protected areas in Paraguay; therefore, monitoring and conservation efforts should be intensified for these taxa [83]. One species recently revalidated is Colobosaura kraepelini. This lizard is known only from the holotype from the locality of Puerto Max (San Pedro Department), the neotype from Altos and an additional specimen from San Bernardino, both localities in Cordillera Department [41]. There are no protected areas in the Cordillera Department, but there are some in the northern portion of Central Department (border with Cordillera), located less than 10 km from the known localities of C. kraepelini. The presence of this species in a conservation unit should be confirmed, but is possible that it is protected by "Monumento Natural Cerro Chororí" and "Monumento Natural Cerro Kõi". It is important to note that the conservation unit closer to the distribution of C. kraepelini is the "Parque Nacional Lago Ypacaraí", although only the lagoon is protected and not the surroundings. The species Homonota septentrionalis was described from the driest part of Paraguay (northwestern Chaco) and is abundant in the "Parque Nacional Teniente Enciso" [44]. The four species of Tropidurus found in Paraguay [43] are well represented in several protected areas. The last herpetofaunal conservation assessment was published in 2009 [84], and thus a new conservation assessment of Paraguayan reptile fauna including the new taxa, is necessary to provide a sound basis for conservation planning for those species that require special attention.
Finally, given that the sequenced specimens are yet a small portion of the actual diversity of Paraguay ( Figure 5), it will be of the utmost importance to continue and expand these studies that will further improve our taxonomic knowledge. published in 2009 [84], and thus a new conservation assessment of Paraguayan reptile fauna including the new taxa, is necessary to provide a sound basis for conservation planning for those species that require special attention.
Finally, given that the sequenced specimens are yet a small portion of the actual diversity of Paraguay ( Figure 5), it will be of the utmost importance to continue and expand these studies that will further improve our taxonomic knowledge.

Supplementary Materials:
The following are available online at www.mdpi.com/xxx/s1, Table S1: Reagents used for molecular procedures; Table S2: Specimens used in the barcoding analysis; Figures S1-S16: Details of sub-trees sections from the general barcoding analysis.
Funding: This research was developed in the framework of the Ph.D. thesis project "Lizards of Paraguay: an integrative approach to solve taxonomic problems in central South America" presented by P.C. in the Faculty of Biosciences of the Johann Wolfgang Goethe-Universität, which had a strong financial support from the Deutscher Akademischer Austauschdienst (DAAD). Additionally, P.C. received a subsidy from the Consejo Nacional de Ciencia y Tecnología (CONACYT) through the Programa Nacional de Incentivo a los Investigadores (PRONII), field equipment provided by Idea Wild, and a grant to perform fieldwork activities in the Mbaracayú Reserve by funded by the GNB Bank through the Fundación Moisés Bertoni (FMB). Additionally, we received financial support from PRESIDENT ENERGY to pay the APC of this work.
Acknowledgments: This was a major project that had the help of many people. For this, we owe gratitude to Norman Scott, Luciano Avila, Mariana Morando, Tony Gamble, Sebastian Lotzkat, Martin Jansen, Arne Schultze, Christian Printzen, Raúl Maneyro, and Abel Batista, who provided advice and technical support at different stages of the project. This project would not have been achieved without the valuable help of many colleagues and friends who offered their help during field work. Hence, thank you so much Jorge Ayala, Figure 5. Comparison between known diversity of Paraguayan Squamata (blue bars) vs. diversity of sampled taxa for this study (red bars).
Funding: This research was developed in the framework of the Ph.D. thesis project "Lizards of Paraguay: an integrative approach to solve taxonomic problems in central South America" presented by P.C. in the Faculty of Biosciences of the Johann Wolfgang Goethe-Universität, which had a strong financial support from the Deutscher Akademischer Austauschdienst (DAAD). Additionally, P.C. received a subsidy from the Consejo Nacional de Ciencia y Tecnología (CONACYT) through the Programa Nacional de Incentivo a los Investigadores (PRONII), field equipment provided by Idea Wild, and a grant to perform fieldwork activities in the Mbaracayú Reserve by funded by the GNB Bank through the Fundación Moisés Bertoni (FMB). Additionally, we received financial support from PRESIDENT ENERGY to pay the APC of this work.
Dante and Rafa for patience and support. Finally, to the CONACYT for providing important tools for scientific research in the country.