Biodiversity decline is transpiring in rates unimaginable, yet biodiversity monitoring is hampered by a plethora of factors including incomplete inventory of baseline diversity and dwindling population of skilled taxonomists. Bridging this impediment in species inventory in groups regarded as “dark taxa”, or hyperdiverse groups which are often understudied in taxonomy, is particularly crucial in the tropics which face insurmountable environmental pressures. While morphology-based alpha-taxonomy remains the gold standard in species identification, DNA-based methods have been showing tremendous potential in accelerating species delineation and discovery in various magnitudes.  Here, we describe a concoction pipeline for generating molecular operational taxonomic units (MOTUs) involving eight molecular species delimitation approaches implemented on highly inconspicuous riparian and aquatic beetles from the Philippines. Prior the MOTU generation, specimens were collected using microhabitat approach, and standard DNA isolation, cytochrome c oxidase subunit I (COI) amplification, and Sanger sequencing protocols were performed. Given the COI barcodes, the pipeline utilizes a combination of four threshold-based molecular species delimitation approaches, namely TaxonDNA, K2P, ABGD, and ASAP, which use sequence alignment as input data. Additionally, the pipeline proceeds with four coalescent-based approaches, namely PTP-ML, PTP-BI, mPTP, and GMYC, which require a tree file as input data. Some approaches entail the use of web servers, while others utilize local softwares. With the selection of algorithmic methods presented, MOTUs can be identified whether by consensus or by majority. Remarkably, this pipeline further delineated beetle species from the Philippines which are almost morphologically indistinct even those documented from the same or adjacent localities.

We thank the Bureau of Fisheries and Aquatic Resources (BFAR) and the concerned local government units for granting us prior informed consent to perform biodiversity sampling. We thank Dr Hendrik Freitag for guidance on taxonomic work. We also thank the reviewers and editors for their helpful, constructive comments. Lastly, we express our utmost gratitude to the University Research Council Big Project Grant (URC 2022-09) and the Open Access Publication Grant of the Ateneo de Manila University, the SC Johnson Student Research and Development Fund for Environmental Leadership of the Ateneo Research Institute of Science and Engineering, and the Linn&#233;Sys: Systematics Research Fund of the Linnean Society of London and the Systematics Association for funding this study.

The starting materials for this concoction pipeline are DNA sequences, which can either be downloaded from online repositories (e.g., BOLD, GenBank) or generated from Sanger sequencing71. Table 1 and Table 2 list the accession numbers for the datasets whose representative results are presented in this paper.
1. Preparatory step for threshold-based approaches: Sequence alignment
<ol>
	<li>Load sequences in MEGA X72. 
	NOTE: MEGA X can be downloaded from https://www.megasoftware.net/
	<ol>
		<li>Click on Align to open a dropdown menu. Click on Edit/Build Alignment and select Create a new alignment. Click on OK to confirm this selection. Select DNA as the datatype for alignment.</li>
		<li>Mouse over the Edit tab, and select Insert Sequence from File. Navigate to the directory with the sequences and select them to be loaded into MEGA.</li>
	</ol>
	</li>
	<li>Align sequences by clicking on Alignment, and then Align by ClustalW73. Continue with default settings and click on OK.</li>
	<li>Manually edit the sequences by trimming both ends and by cleaning any insertions or deletions.
	<ol>
		<li>To delete insertions, click on the inserted bases and/or positions, and press Delete on the keyboard.</li>
		<li>To correct deletions, click on the position supposed to have been deleted as denoted by "-". Delete the "-" and type the intended base.</li>
		<li>Find the earliest position where all sequences have a character. Click on the blank box on the row header of the position to the left of the determined position, and drag until the excess positions for all sequences until the starting end are selected. Press Delete to trim.</li>
		<li>Repeat step 1.3.3 by looking for the last position where all sequences have a character. Click on the blank box on the header row of the position to the right of the determined position, drag until the excess positions for all sequences until the terminal end are selected. Press Delete. 
		NOTE: Figure 2A shows a good alignment of DNA sequences. The final dataset is a matrix of sequences of equal length; thereby, sequence length is not a factor that affects the outcomes of the study. The reference sequences serve as a basis for the final length of the alignment matrix.</li>
	</ol>
	</li>
	<li>Translate the sequence to check for stop codons, given that COI is a protein-coding sequence.
	<ol>
		<li>Select all sequences, then click on the Translated Protein Sequences tab.</li>
		<li>When prompted, verify the genetic code as Invertebrate Mitochondrial. If the genetic code is different, click No, and a menu will appear allowing one to tick the box for Invertebrate Mitochondrial genetic code.</li>
		<li>If stop codons, as denoted by asterisks in the alignment body, are present for an entire column, click on DNA Sequences and delete the first position for all sequences.</li>
		<li>Repeat the previous step if stop codons are still present in any sequence.</li>
		<li>If stop codons are still present, restore the previously deleted positions. Select all sequences, then mouse over the Data tab and click on Reverse.</li>
		<li>Repeat step 1.4.1. If stop codons are still present, repeat steps 1.4.3 and 1.4.4. 
		NOTE: Figure 2B shows a good alignment of codons.</li>
	</ol>
	</li>
	<li>Click on DNA Sequences, and save the sequence alignment to .mas/x. Export the sequence into other pertinent file types, such as .meg file and .fasta file. 
	NOTE: Certain molecular species delimitation approaches work better without the outgroup. Thus, a possible consideration is generating a copy of the alignment file without the outgroup.</li>
</ol>
<img alt="Genetic sequence alignment, MAFFT software, DNA to protein analysis, comparative genomics." class="xfigimg" src="/files/ftp_upload/68323/68323fig2.jpg" data-altupdatedat="1788304762000" data-alt="Figure 2"> 
Figure 2: Sequence alignment in MEGA X. (A) Manually cleaned final DNA sequence alignment and (B) translated protein sequence alignment of the Byrrhinus dataset showing no insertions, deletions, and gaps. <a href="https://www.jove.com/files/ftp_upload/68323/68323fig2large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
2. Delimitation by Taxon DNA module in Species Identifier 1.8
NOTE: Species Identifier can be downloaded from https://taxondna.sourceforge.net/
<ol>
	<li>Open TaxonDNA. Click on Import, then click on FASTA. Upload the alignment in .fasta format.</li>
	<li>Click on Modules, and click on Cluster. Set the threshold to 3%. Then, check on Generate information about individual cluster (Figure 3).</li>
	<li>Click on Make clusters now!. Save the results by taking a screenshot.</li>
</ol>
<img alt="DNA sequence clustering interface for taxonomy analysis; shows sequence clustering results and statistics." class="xfigimg" src="/files/ftp_upload/68323/68323fig3.jpg" data-altupdatedat="1788304762000" data-alt="Figure 3"> 
Figure 3: TaxonDNA in Species Identifier 1.8 window showing the 3% user-identified threshold for delimitation. <a href="https://www.jove.com/files/ftp_upload/68323/68323fig3large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
3. Delimitation by Kimura 2-parameter (K2P) in MEGA X
<ol>
	<li>In MEGA, open the .meg file by clicking on File, then Open A File/Session. Click on Distance, then click on Compute Pairwise Distances. Confirm the .meg file for delimitation.</li>
	<li>Mouse over the yellow box beside Model/Method, then click the arrow that appears on the right of the box.</li>
	<li>Select Kimura 2-parameter model from the dropdown menu. Click on Ok to run the program. Open the data output window (Figure 4).</li>
	<li>Click on File, then Export/Print Distances.</li>
	<li>Open the dropdown menu beside Output Format and select XL: Microsoft Excel workbook. For Decimal Places, choose 4.</li>
	<li>Save the resulting spreadsheet in Microsoft Excel. Use Conditional Formatting to highlight those with a distance of less than 0.03 (or 3%). 
	​NOTE: Unlike the other approaches discussed in this article, K2P does not immediately provide the clustering because it is dependent on the user-identified threshold. This pipeline is working on the assumption of the threshold being 3%, which means that any pair of sequences whose genetic distance is less than 3% belongs to the same molecular cluster.</li>
</ol>
<img alt="Phylogenetic distance matrix, DNA comparison, table, species genetic differences." class="xfigimg" src="/files/ftp_upload/68323/68323fig4.jpg" data-altupdatedat="1788304762000" data-alt="Figure 4"> 
Figure 4: Resulting window in MEGA X showing the outcomes of Kimura-2-parameter. <a href="https://www.jove.com/files/ftp_upload/68323/68323fig4large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
4. Delimitation by Assemble Species by Automatic Partitioning (ASAP)
<ol>
	<li>Go to the ASAP web server: https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html</li>
	<li>Upload the alignment by clicking on the orange box labeled Choose a file… and selecting the .fasta file. Alternatively, drag the .fasta file to the orange box. Scroll down, and click on Go.</li>
	<li>To download the clustering, click on list for the row with the lowest ASAP-score and highest p-val rank (Figure 5). 
	​NOTE: The result output is a .txt file, which presents subsets of the molecular clusters.</li>
</ol>
<img alt="Statistical analysis results; dataset table, histogram, ranked distances graph; data distribution analysis." class="xfigimg" src="/files/ftp_upload/68323/68323fig5.jpg" data-altupdatedat="1788304762000" data-alt="Figure 5"> 
Figure 5: Resulting window in ASAP showing the possible delimitations. The first outcome (blue box) is accepted because of the lowest ASAP-score (red box) and highest P-val rank (yellow box). <a href="https://www.jove.com/files/ftp_upload/68323/68323fig5large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
5. Preparatory step for coalescent-based approaches: Tree estimation
<ol>
	<li>Open the .meg file by clicking on Data, then Open A File/Session. Test the best-fit substitution model.
	<ol>
		<li>Click on Models, then find Best DNA/Protein Models (ML). Confirm the .meg file for analysis.</li>
		<li>Click on Ok on the Analysis Preferences menu (Figure 6A). Take note of the best model using the lowest BIC and AICc value (Figure 6B). 
		NOTE: This is usually the first in the list of models in the resulting window. Refer to the complete names of abbreviated model names at the bottom part of the resulting window.</li>
	</ol>
	</li>
	<li>Generate the tree.
	<ol>
		<li>Click on Phylogeny and then Construct/Test Maximum Likelihood Tree.</li>
		<li>Use the best-fit nucleotide substitution model. Do this by clicking on the yellow box beside Model/Method. Select the model from the dropdown menu.</li>
		<li>If the best fit nucleotide substitution model has '+G' or '+I' or both, click on the box beside Rates among Sites and select the parameters (G, I, G + I) from the dropdown menu (Figure 6C).</li>
		<li>Under Phylogeny Test, click on the box beside Test of Phylogeny and select Bootstrap method from the dropdown menu. Click on the box beside "No. of Bootstrap Replications" and type 1000 (Figure 6C).Click on Ok to run the analysis.</li>
	</ol>
	</li>
	<li>Once a window opens containing the resulting tree (Figure 6D), save the tree session as .mts/x. Save the output as .nwk file. Save also the tree as .png picture file.</li>
</ol>
<img alt="Phylogenetic analysis method settings, likelihood table, and tree diagram; maximum likelihood, bootstrap." class="xfigimg" src="/files/ftp_upload/68323/68323fig6.jpg" data-altupdatedat="1788304762000" data-alt="Figure 6"> 
Figure 6: Tree estimation in MEGA X as a preparatory step for coalescent-based approaches. (A) Setting to determine the best-fit nucleotide substitution model. (B) For the Byrrhinus dataset, the best-fit model is GTR + I due to the lowest BIC and AICc values (red box). (C) In running the maximum likelihood tree, GTR (yellow box) + I (blue box) were set as parameters. (D) The resulting raw ML tree to be used in the subsequent analyses. <a href="https://www.jove.com/files/ftp_upload/68323/68323fig6large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
6. Delimitation by Poisson Tree Processes (PTP)
NOTE: The steps below will yield the outcomes of PTP-ML and PTP-BI.
<ol>
	<li>Go to the PTP web server: https://species.h-its.org/ptp/</li>
	<li>Upload a newick file by clicking on Choose file. Under My tree is, select Rooted (Figure 7A).</li>
	<li>Input the parameters for the analysis. Under No. MCMC generations, input 100000. Under Thinning, input 100. Under Burn-in, input 0.1. Under Seed, input 123.</li>
	<li>Input the outgroup by typing the name of the tips in the box under Outgroup taxa names (if any). Instructions for formatting input are to the right of the input box.</li>
	<li>Input email address in the corresponding field, and click Submit. In the PTP species delimitation results window, click on here to open the result window.</li>
	<li>Under the Maximum likelihood solution, download the results of PTP-ML by clicking on Download delimitation results (Figure 7B).</li>
	<li>Under the Bayesian-supported solution, download the results of PTP-BI by clicking on Download delimitation results. 
	​NOTE: The results would indicate groups of sequences in partitions. These will be labeled as numbered species followed by a support value, wherein 1.000 is ideal.</li>
</ol>
<img alt="bPTP species delimitation results, phylogenetic tree diagram, MCMC settings, maximum likelihood analysis." class="xfigimg" src="/files/ftp_upload/68323/68323fig7.jpg" data-altupdatedat="1788304762000" data-alt="Figure 7"> 
Figure 7: Molecular species delimitation by Poisson tree processes. (A) Windows showing the input parameters. (B) Resulting window in PTP web server, with link containing delimitation results (red box). <a href="https://www.jove.com/files/ftp_upload/68323/68323fig7large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
7. Delimitation by Multi-rate Poisson Tree Processes (mPTP)
<ol>
	<li>Go to the mPTP web server: https://mptp.h-its.org/#/tree</li>
	<li>Upload the newick file by dragging it onto the gray rectangle or clicking on the rectangle. Once the data loads on the same page, click on Proceed to outgroup selection (Figure 8A).</li>
	<li>On the Outgroup specification page, select the outgroup by clicking on the box beside the taxa labels of outgroup specimens. Click on the Crop Outgroup option at the bottom of the sequence. Click on Model selection.</li>
	<li>Click on MPTP, then click on Visualization Options. On the Visualization options page, accept the default settings. 
	NOTE: Alternatively, the total image size may be adjusted by inputting a value (in px) under SVG width. Font and spacing of text can be adjusted by inputting desired sizes (in px) under SVG font size, Space between tips, Left Margin, Right Margin, Top Margin, and Bottom margin.</li>
	<li>Click on Contact details. Input the email address for future reference of results (Optional).</li>
	<li>Click on Overview. Make sure to review the number of taxa input and the outgroups selected.</li>
	<li>Click on Submit. Right-click the files under Downloadable Files (Figure 8B), and select Save as… to save the results. 
	NOTE: The output text file will indicate the total number of delimited species at the beginning of the results. It will then present which sequences fall under each of the clusters. These clusters and corresponding sequences are arranged similarly to the generated tree in the .svg file. The mPTP SVG file will display the clustering of the sequences. Clusters are presented by groups of sequences connected by red branches along the tree. Each isolated set of red branches indicates a singular MOTU.</li>
</ol>
<img alt="Phylogenetic analysis interface for tree upload, outgroup selection, and file parsing in research setup." class="xfigimg" src="/files/ftp_upload/68323/68323fig8.jpg" data-altupdatedat="1788304762000" data-alt="Figure 8"> 
Figure 8: Molecular species delimitation by multi-rate Poisson tree processes. (A) Windows showing the mPTP-read information regarding the uploaded tree. (B) Resulting window in mPTP web server, with link containing delimitation results (red box). <a href="https://www.jove.com/files/ftp_upload/68323/68323fig8large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
8. MOTU generation
<ol>
	<li>Open the tree using a photo editor program or using PowerPoint.</li>
	<li>Create a bar to represent the results of each molecular species delimitation approach (see Figure 9 and Figure 10 for samples). Ensure that any single bar aligns with the sequences it includes in its molecular cluster.</li>
	<li>Determine the MOTU by one of these two approaches74
	<ol>
		<li>If all the approaches generate identical results for a given molecular cluster, erect the molecular cluster as a MOTU by consensus.</li>
		<li>If most approaches generate identical results for a given molecular cluster, erect the molecular cluster as a MOTU by majority.</li>
	</ol>
	</li>
</ol>

We have nothing to disclose.

While DNA sequences, especially COI49, may be compared against reference sequences22,53,78 in online public repositories54,55,79,80 for a bioidentification system46, such an approach is limited by the available deposited sequences. On the contrary, the use of molecular species delimitation can set species limits and form molecular clustering without the need for a sequence repository as a prerequisite56. Given the two classes of molecular species delimitation approaches, there are two critical steps in the concoction pipeline, which correspond to the preparation of the input file, namely the assembly of the sequence alignment and the generation of the gene tree.
As demonstrated here, the concoction pipeline, when applied to both Byrrhinus and Anacaena datasets, generated MOTUs that can provide clues on species identification. The concoction pipeline expanded the support for the delineation of known species. First, the pipeline provided additional support for B. negrosensis and B. villarini, which were previously delineated using only one approach, namely K2P35, by having concordant outcomes from five other approaches. This is particularly striking because freshwater beetles are usually assumed not to coexist in the same sampling locality with their congeners, yet outcomes of the concoction pipeline provide stronger support for syntopy or congeners being collected from the same two rivers35. This trend agrees with observations on other endemic beetle species51,81,82. Instances like syntopy punctuate the importance of proper molecular species delimitation to challenge said assumption. There are instances where beetle richness was previously underestimated due to poor delineation approaches, which led to some challenges in conservation8. Second, the concoction pipeline expanded the support for delineation of A. angatbuhay and A. auxilium, which were previously delineated only by morphology-based alpha-taxonomy70. The results (Figure 10) demonstrate that the pipeline validates the morphological species limits.
The delimitation worked for both the 723bp-long 3&#39;-end and 658bp-long 5&#39;-end of the COI gene as shown in the Byrrhinus (Table 1) and Anacaena (Table 2) datasets, respectively. Some molecular taxonomic studies used much shorter COI segments of &#62;200 bp83,84 owing to justified amplification capacity and sound species-specific signals. Aside from COI, mitochondrial genes COB8,85, and 16S38,86 were previously used in the exploration of species delimitation. This concoction pipeline, however, maintains the use of the full length of COI-3&#39; and COI-5&#39; Folmer region following the common lengths of the reference sequences for the aquatic and riparian beetle taxa of interest in this study35,76,77.
Moreover, this article assembled a concoction of independent molecular species delimitation approaches, such that there is stronger support for delineation. Each method in this concoction pipeline is chosen because of the diversity of assumptions they have. Thus, by using a concoction of different approaches, the clustering in the final MOTU is regarded as robust. Three approaches used in this pipeline are threshold-based methods, which use cut-off values to determine divergence between and within species47,58. For one, TaxonDNA through its Species Identifier module60 functions by computing pairwise uncorrected distance between sequences and clustering them according to the threshold set by the user46,49. Similarly, K2P works on pairwise comparison of genetic distance between the sequences61. While thresholds are identified by the user in Species Identifier and K2P, ASAP presents different clusters by generating various intraspecific diversity p-values for the user and ranks the different outcomes62. Moreover, there are other molecular species delimitation approaches beyond those featured in this concoction pipeline. Other threshold-based approaches include Automatic Barcode Gap Discovery (ABGD)57, refined single linkage (RESL) method as implemented in the Barcode of Life Data (BOLD) Systems55, Structure87, Structurama88, and Gaussian mixed clustering89.
Moving away from comparing nucleotide bases as mere diagnostic sequences, such as in threshold-based approaches, coalescent-based species delimitation approaches identify species boundaries and sort independent lineages by determining change in lineage branching rate of the gene or phylogenetic tree63,64,65,66,90. The approaches implemented in this study are closely related with&#160;the difference that PTP assumes the fixed identical evolutionary rate throughout the tree, while multi-rate PTP assumes that the evolutionary rate differs in different branches68,69. Beyond these, other examples of such approaches are generalized mixed Yule-coalescent (GMYC)56,64,91, Bayesian multispecies coalescent92, and spedeSTEM93.
As presented earlier (Figure 9 and Figure 10), in the context of this study, which focused on aquatic and riparian beetles, there is concordance across the outcomes of the six approaches, as well as between delimitations using morphology-based alpha-taxonomy and molecular species delimitation approaches. As such, the accuracy and precision of the concoction pipeline are parallel to alpha-taxonomy, which is regarded as the gold standard in beetle taxonomy. While the outcomes in this study yielded concordant results, which led to the formulation of MOTUs by consensus74, discordance among approaches is possible27,32,63. In instances where different molecular species delimitation approaches yielded different clustering, MOTUs can still be erected, given that the majority of the approaches generated identical clustering32,74.
For example, hypothetically, four approaches clustered eight sequences as one molecular cluster, while two approaches split the said sequences into two molecular clusters. Should this happen, the eight sequences form one MOTU, following MOTU the majority. The two discordant outcomes may later be examined as this may hint at ongoing or upcoming species divergence94. In any case, given that this concoction pipeline is intentionally designed to feature six approaches, the majority is four. This number guarantees that MOTUs erected have identical clustering for at least one threshold-based and at least one coalescent-based approach, which suggests that such molecular clustering is robust and can withstand the different and variable assumptions explained earlier74,86,95.
The number of investigations on the different molecular species delimitation approaches has increased tremendously in the past decade, especially on simulations optimizing the algorithm. Studies include comparison of approaches designed for single-locus applied to multilocus datasets67, comparison of multilocus analyses7, identifying optimal settings for species-poor data96 and species-rich data97, managing conflicting gene trees for tree-based approaches98, and comparison of sensitivity to sampling effort and non-monophyly74 and to singletons99. While robust outcomes have consistently been extracted in datasets from varying taxa, studies have identified that proper delimitation is affected by a confluence of factors, including mutation rate of marker, number of loci, speciation rate, species richness, presence of singletons, effective population size, population size to divergence time ratio, sample size per species, and dispersal ability57,67,96,97. Fortunately, the effects of these items were not noted in this study.
There are a lot of permutations on how molecular species delimitation approaches can be operationalized in taxonomy. While some studies tend to validate the high congruence of threshold-based approaches to morphological data27,74,100, many studies have heavy reliance on coalescent-based approaches38,67,101. For one, up until recently, threshold-based approaches have been widely ignored in simulation-based comparison of species delimitation methods97.
Adding to the multitude of permutations, the placement of molecular approaches in integrative taxonomic studies on different taxa is highly variable. Traditional taxonomic studies de facto use morphological data as primary clustering to be validated by molecular species delimitation85,102. With molecular methods becoming more accessible, this concoction pipeline supports the call of numerous studies to utilize genetic data for primary clustering, which is further delineated using morphological data52,103,104,105,106. Moreover, this concoction pipeline is similar to the protocol of automatic barcode gap discovery (ABGD)57, molecular biodiversity assessment52, and large-scale integrative taxonomy (LIT)27 in using threshold-based delimitation approaches to generate molecular clustering of COI barcodes. In LIT, the potentially incongruent clusters with high genetic distance are then primarily verified by examining a second line of evidence, such as morphological data. While the three methods compared the efficiency of their clustering with other algorithms, including coalescent-based approaches, the concoction pipeline described here necessitates the use of both threshold-based and coalescent-based approaches for a more robust delimitation. This is similar to the MOTU estimation phase of the recently proposed dark taxonomy pipeline106.
Even with this DNA data-dependent concoction pipeline, the gold standard in erecting new beetle species remains to be morphology-based alpha-taxonomy. Given the resources and specialized alpha-taxonomy skills6,107,108 necessary to document the multitude of undiscovered and undescribed species109,110,111, it was predicted that more than $260 billion and 360 years would be needed to completely elucidate the diversity of life on earth using alpha-taxonomy alone112. With these challenges associated with species discovery, the majority of biodiversity studies are focused on megafauna and iconic species113,114, which further creates taxonomic bias. Thus, the use of the concoction pipeline based on DNA data to formulate preliminary clustering, which can then be examined for morphological similarities and differences, is relevant because it can accelerate species discovery. As shown in this article, such an approach can save taxonomists from the &#39;shotgun&#39; approach of having to sort a multitude of morphospecies of megadiverse and highly indistinguishable aquatic and riparian beetles, and can proceed in species inventory and discovery in a more efficient manner.

With the unprecedented rates of global species decline across varied taxa1,2,3,4, the impetus to create a comprehensive biodiversity inventory, which includes previously undiscovered and undescribed species, becomes a race against time2,5. For one, biodiversity conservation only makes sense in light of proper taxonomic knowledge6,7. Thus, the identification of biodiversity at the species level is necessary, given that closely related yet different species may have unique fundamental niches8.
In beetle systematics, alpha-taxonomy or traditional taxonomy is the gold standard in species identification similar to most, if not all, animal groups9. This approach relies on morphological characters to classify organisms at higher taxonomy and to provide identification at the species-level10,11. In many instances, family12 or even genus-level identification13 can be done using external morphological features alone, such as body form and specialized adaptations. Meanwhile, species-level identification in beetles is usually done via the comparison of the structural details of the aedeagus or the male genitalia14,15,16.
Despite the wide acceptance of this morphology-based alpha-taxonomy, this approach to species discovery is hampered by numerous 'taxonomic impediments'6,17,18, especially among invertebrates19. In beetle systematics, species identification and/or discovery using alpha-taxonomy faces the challenges of having a limited number of skilled taxonomists20,21, of brief and barely informative earlier descriptions22, and of logistical problems related to access to literature and type specimens23.
This is even exacerbated when the taxon in question is cryptic or highly inconspicuous. Cryptic species refers to the set of species of the same genus, or even different genera, which cannot be easily delineated by comparative morphology, given subtle interspecific phenetic differences24. Some cryptic taxa, regarded as 'dark taxa'25,26,27, additionally suffer from being hyperdiverse yet heavily understudied. Unfortunately, this is not uncommon in riparian and aquatic beetles8,28. Failure to address species delineation among cryptic dark taxa threatens biodiversity as a one-size-fits-all conservation measure is given to organisms of potentially different niches and ecological requirements8,29.
Given these challenges, 'integrative taxonomy'23, or the approach of coupling alpha-taxonomy with another line of evidence, has gained traction as a means to erect new species in the last two decades30,31,32. One such line of evidence used in integrative taxonomic studies on beetles33,34,35,36,37,38,39 and other insect orders40,41,42,43,44,45 is DNA sequence. In particular, the mitochondrial cytochrome c oxidase subunit I (COI) gene is being used to barcode a diverse set of animals due to its moderately conservative nature46,47,48. While its 658bp-long 5′-end (COI-5′), also known as the 'barcoding fragment' or 'Folmer fragment', has been proposed as the 'diagnostic' sequence49, the 723 bp-long COI-3′ is also being used owing to a robust set of primers for the 3′-end50,51,52.
Whereas DNA barcoding using COI was conceived for species-level identification, it is limited by an a priori sequence repository53, such as GenBank54 and BOLD55, molecular species delimitation approaches set species limits using only the provided DNA sequences with no pre-requisite repository56. There are two classes of molecular species delimitation approaches based on the input file. First, using sequence alignment as the input, threshold-based approaches use cut-off values to determine divergence between and within species48,57,58. As distance-based methods, these approaches rely on an observed maximum intraspecific distance of 3% for insects, which is being referred to as the 'barcoding gap'51,59. Some of the threshold-based approaches are SpeciesIdentifier implemented in TaxonDNA60, genetic distance via Kimura 2-parameter (K2P)61, and Assemble Species by Automatic Partitioning (ASAP)62.
Second, using the tree as the input file, coalescent-based approaches identify species boundaries and sort independent lineages by determining change in lineage branching rate63,64,66. Often relying on phylogenetic species concept67, coalescent-based approaches include Poisson tree processes (PTP)68 and their variation, multi-rate PTP69. In any case, the clusters formed by the molecular species delimitation approaches are collectively referred to as molecular operational taxonomic units (MOTUs).
Here, we present methods for a concoction pipeline for generating MOTUs involving three threshold-based approaches and three coalescent-based approaches (Figure 1). The pipeline is tested on two datasets, namely (1) COI-3′ sequences of the riparian Byrrhinus Motschulsky, 1858 beetles, which is composed of two published species35 previously erected using integrative taxonomy and additional sequences from undescribed species, and (2) COI-5′ sequences of Anacaena Thomson, 1859 beetles, which is composed of sequences from two published species70 previously erected using alpha-taxonomy and additional sequences from undescribed species.
<img alt="DNA sequence analysis flowchart; alignment, threshold, coalescent methods, tree estimation, MOTU steps." class="xfigimg" src="/files/ftp_upload/68323/68323fig1.jpg" data-altupdatedat="1788304762000" data-alt="Figure 1"> 
Figure 1: Flowchart of the concoction pipeline for generating MOTUs. <a href="https://www.jove.com/files/ftp_upload/68323/68323fig1large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>

<table><tbody><tr><td>ASAP&amp;nbsp;</td><td>Spart explorer</td><td>https://bioinfo.mnhn.fr/abi/public/asap/asapweb.html</td><td></td></tr><tr><td>Computer</td><td></td><td>Any processing unit</td><td></td></tr><tr><td>MEGA X&amp;nbsp;</td><td>MEGA Software</td><td>https://www.megasoftware.net/</td><td></td></tr><tr><td>mPTP web server</td><td></td><td>https://mptp.h-its.org/#/tree</td><td></td></tr><tr><td>PTP web server</td><td>Heidelberg Institute for Theoretical Studies (HITS)</td><td>https://species.h-its.org/ptp/</td><td></td></tr><tr><td>Species Identifier&amp;nbsp;</td><td></td><td>&amp;nbsp;https://taxondna.sourceforge.net/</td><td></td></tr></tbody></table>

Biodiversity decline is transpiring at rates unimaginable, yet biodiversity monitoring is hampered by a plethora of factors, including the incomplete inventory of baseline diversity and the dwindling population of skilled taxonomists. Bridging this impediment in species inventory of &#34;dark taxa&#34;, or hyperdiverse groups often understudied in taxonomy, is particularly crucial in the tropics facing insurmountable environmental pressures. While morphology-based alpha-taxonomy remains the gold standard in species identification, DNA-based methods have been showing tremendous potential in accelerating species delineation and discovery. Here, we describe a concoction pipeline for generating molecular operational taxonomic units (MOTUs) involving six molecular species delimitation approaches implemented on&#160;two highly inconspicuous beetle genera, namely Byrrhinus Motschulsky, 1858 of family Limnichidae and Anacaena Thomson, 1859 of family Hydrophilidae, from the Philippines. As this protocol focuses on processing DNA barcodes and given that these genetic data can be obtained either by sequencing or by downloading from online databases, the starting point for the protocol described here are the DNA sequences. Given the COI barcodes, the pipeline utilizes a combination of three threshold-based molecular species delimitation approaches, namely TaxonDNA, K2P, and ASAP, which use sequence alignment as input data. Additionally, the pipeline proceeds with three coalescent-based approaches, namely PTP-ML, PTP-BI, and mPTP, which require a tree file as input data. Some approaches entail the use of web servers, while others utilize local software. With the selection of algorithmic methods presented, MOTUs can be identified, whether by consensus or by majority. Remarkably, this pipeline delineated almost conspicuous beetle species from the Philippines, even those documented from the same or adjacent localities.

This article implemented a molecular species delimitation concoction pipeline, which includes TaxonDNA, K2P, ASAP, two variations of PTP, and mPTP, to aquatic and riparian beetles. The outcome of each of the six approaches is referred to as 'molecular clusters', while the summary is referred to as 'molecular operational taxonomic units' or MOTUs75.
The concoction pipeline was tested on two datasets. The first dataset (Table 1, Figure 9) is composed of 24 COI-3' sequences of the riparian Byrrhinus Motschulsky, 1858 beetles collected from the Philippines, with two unidentified Limnichus Latreille, 1829 species76 as the outgroup. The Byrrhinus sequences include sequences from two previously published species35 erected using morphology-based alpha-taxonomy with K2P delimitation, as well as sequences from unknown species collected from five different provinces.
<img alt="Phylogenetic tree diagram illustrating genetic relationships; includes taxonomic classification." class="xfigimg" src="/files/ftp_upload/68323/68323fig9.jpg" data-altupdatedat="1788304762000" data-alt="Figure 9"> 
Figure 9: Consensus COI-3′ gene tree of Philippine Byrrhinus following GTR+I model using maximum likelihood. Values on nodes are bootstrap support values, BS &lt; 80 are hidden. Color bars summarize the molecular clusters from the six molecular species delimitation approaches. Yellow bars represent the molecular clusters of threshold-based approaches, namely TaxonDNA, Kimura 2-parameter (K2P), and Assemble Species by Automatic Partitioning (ASAP). Blue bars represent the molecular clusters of coalescent-based approaches, namely two iterations of single Poisson Tree Processes (PTP-BI and PTP-ML) and multi-rate PTP (mPTP). The outcomes of the six approaches are summarized as molecular operational taxonomic units (MOTUs, in black bar). <a href="https://www.jove.com/files/ftp_upload/68323/68323fig9large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
Figure 9 presents the molecular clusters of the threshold-based approaches in shades of yellow and the molecular clusters of the coalescent-based approaches in shades of blue. Consider the first vertical bar under TaxonDNA. The figure can be read to mean that six sequences (EDD102, EDD187, EDD441, EDD044, EDD522, EDD413) were designated to belong to the same molecular cluster by the approach TaxonDNA. Meanwhile, the black bar is the summary of the six molecular species delimitation approach.
With this, all six molecular species delimitation approaches yielded the same delimitations. There are four MOTUs retrieved from the Byrrhinus dataset, namely the two known species from the island of Negros, namely B. negrosensis Delocado &amp; Freitag, 2021 and B. villarini Delocado &amp; Freitag, 2021, and two unknown species labelled as 'Byrrhinus sp. A' and 'Byrrhinus sp. B'. This clustering is standard for all six methods, and the MOTUs can be retrieved by consensus74. Remarkably, the concoction pipeline was able to delineate the sequences of B. negrosensis and B. villarini despite having been collected from the same localities (Table 1). Additionally, the pipeline clustered sequences from four different provinces into one MOTU, namely Byrrhinus sp. A.
<img alt="Phylogenetic tree diagram; DNA barcoding of Anacaena species with clustering analysis results." class="xfigimg" src="/files/ftp_upload/68323/68323fig10.jpg" data-altupdatedat="1788304762000" data-alt="Figure 10"> 
Figure 10: Consensus COI-5′ gene tree of Philippine Anacaena following T92+G+I model using maximum likelihood. Values on nodes are bootstrap support values, BS &lt; 80 are hidden. Color bars summarize the molecular clusters from the six molecular species delimitation approaches. Yellow bars represent the molecular clusters of threshold-based approaches, namely TaxonDNA, Kimura 2-parameter (K2P), and Assemble Species by Automatic Partitioning (ASAP). Blue bars represent the molecular clusters of coalescent-based approaches, namely two iterations of single Poisson Tree Processes (PTP-BI and PTP-ML) and multi-rate PTP (mPTP). The outcomes of the six approaches are summarized as molecular operational taxonomic units (MOTUs, in black bar). <a href="https://www.jove.com/files/ftp_upload/68323/68323fig10large.jpg" target="_blank">Please click here to view a larger version of this figure.</a>
Meanwhile, the second dataset (Table 2, Figure 10) is composed of 23 COI-5' sequences of the aquatic Anacaena Thomson, 1859 beetles collected from the Philippines, with two Agraphydrus Regimbart, 1903 species from the Philippines77 as the outgroup. Four MOTUs were retrieved for the dataset. This includes two known species, namely A. angatbuhay Sanchez, Delocado, &amp; Freitag, 2022 and A. auxilium Sanchez, Delocado, &amp; Freitag, 2022, which were previously erected using alpha-taxonomy only and without any molecular data70. The other two MOTUs erected were from unknown species labelled as 'Anacaena sp. A' and 'Anacaena sp. B'. With the clustering being identical for all molecular species approaches, similar to the first dataset, the four MOTUs were retrieved by consensus74. Notably, there is an overlap in the distribution of A. auxilium with Anacaena sp. A and Anacaena sp. B. This shows that geographic distribution does not strongly influence molecular clustering unless the geographic pattern is supported by genetic structures and lineage branching, which was not the case in the datasets.
Table 1: GenBank and BOLD accession numbers of COI-3′ sequences and locality of Byrrhinus specimens used for the genetic distance analysis. <a href="https://www.jove.com/files/ftp_upload/68323/68323_Table%201.xlsx" target="_blank">Please click here to download this Table.</a>
Table 2: GenBank and BOLD accession numbers of COI-5′ sequences and locality of Anacaena specimens used for the genetic distance analysis. <a href="https://www.jove.com/files/ftp_upload/68323/68323_Table%202.xlsx" target="_blank">Please click here to download this Table.</a>

The protocol presented here aims to guide the users on the computer-based processing of cytochrome c oxidase subunit I (COI) gene sequences generated from beetles, such that the species clustering hypothesis called molecular operational taxonomic units (MOTUs) can be generated from DNA.

a-concoction-pipeline-for-generating-molecular-operational-taxonomic

A Concoction Pipeline for Generating Molecular Operational Taxonomic Units (MOTUs) Among Riparian and Aquatic Beetles

Biodiversity decline is transpiring at rates unimaginable, yet biodiversity monitoring is hampered by a plethora of factors, including the incomplete ...

A Concoction Pipeline for Generating Molecular Operational Taxonomic Units (MOTUs) among Riparian and Aquatic Beetles

Physilia Chua

Effective biodiversity monitoring is critical for understanding environmental changes and implementing appropriate management strategies to reverse biodiversity loss. Advances in DNA-based high-throughput sequencing technologies such as metabarcoding, metagenomics, and target capture have revolutionized the way we carry out biodiversity monitoring. This collection aims to accelerate the use of these DNA-based methods for cataloging our terrestrial biodiversity across the Tree of Life, with possible applications, including species discovery and species interactions. We welcome all research using environmental DNA, specimens, or bulk collected samples to assess or monitor biodiversity.&nbsp;

JoVE Biology welcomes all general biology research methodologies. Content in this section canvases all fields of cell, molecular, and organismal biology, ranging from new applications of standard techniques to novel approaches aimed at understanding the functions of life and living organisms. This diverse section includes, but is not limited to, techniques in physical biology, cellular biochemistry, genetics, physiology, systems biology and a combination of eukaryotic and prokaryotic model systems.

Biology

Advancing Terrestrial Biodiversity Monitoring Through DNA-Based Methods

Effective biodiversity monitoring is critical for understanding environmental changes and implementing appropriate management strategies to reverse biodiversity loss. Curated by Physilia Chua. Explore Biology techniques on JoVE.

Guest Editor

Collection Overview

Articles

A Concoction Pipeline for Generating Molecular Operational Taxonomic Units (MOTUs) Among Riparian and Aquatic Beetles