The genus Sylvirana includes 12 species widely distributed in South China and Southeast Asia. The phylogenetic relationships and species diversity for Sylvirana and allied genera remain unresolved and controversial due to insufficient data and incomplete taxon sampling. Using a combined dataset of mitochondrial genes (16S and COI) and 101 nuclear genes obtained through the amplicon sequence capture approach, we generated the most comprehensive phylogenetic analysis for the genus Sylvirana to date, inferring diversity, phylogenetic relationships, and historical biogeography with unprecedented levels of taxon and geographic sampling. Our results conservatively reveal six undescribed species, mostly distributed in peninsular Indochina. Phylogenetic analyses strongly support the non-monophyly of Sylvirana with respect to Pterorana. Additionally, phylogenetic results place Sylvirana guentheri and Pelophylax lateralis into genus Humerana, supporting the inclusion of Hylarana latouchii, Papurana milleti, and Hylarana attigua within Pterorana + Sylvirana. The long-disputed species of Hylarana bannanica (previously Sylvirana) cluster with genus Papurana. Because the results of multiple non-monophyletic genera create taxonomic confusion, we suggest relegating all genera to subgenus rank of Hylarana. Sylvirana is a junior synonym of the Pterorana. Biogeographically, we trace the origin of Pterorana to Southeast Asia during the early Miocene, with subsequent dispersal thereafter. Our study shows that climatic changes may have profoundly influenced the diversification of Pterorana during the Miocene.

Library preparation, bait preparation, hybridization, and sequencing were conducted at Sun Yat-sen University. The protocol followed that of Zhang et al. (2019). In brief, we mixed 20 ng of genomic DNA of each of the 38 samples to make a DNA pool. The pooled DNA was used as template to amplify 101 nuclear protein-coding (NPC) markers using the primers and protocol developed by Shen et al. (2013). The 101 PCR products were then mixed together and purified. The purified PCR mixture was subjected to end-repair and A-tailing and then ligated with a biotinylated adapter. The biotinylated amplicons were subsequently immobilized on Dynabeads MyOne streptavidin magnetic beads (Life Technologies) to obtain bait-coated beads.

For each hybridization reaction, 500 ng of pooled libraries and 2.5 μl of bait-coated beads (containing 100 ng of biotinylated amplicons) were used. A touch-down hybridization program was adopted. After denaturation, the hybridization started from 65°C, decreased by 5°C every 6 hours and ended at 45°C, for a total duration of 30 hours. The captured libraries were amplified, pooled in equal concentrations, and sequenced on an Illumina HiSeqX10 sequencer using the paired-end 150-bp mode.

Data processing followed Zhang et al. (2019). In brief, the Illumina paired-end reads were quality-controlled using Trimmomatic (Bolger et al., 2014), FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) and FastUniq (Xu et al., 2012). Clean reads of each sample were assembled into contigs using the SPAdes v.3.8.1 genome assembler (Bankevich et al., 2012). The sequencing depths for the filtered contigs were calculated by SAMtools v.1.4.1 (Li et al., 2009). Only contigs with an average sequencing depth ≥ 5× were retained for further analysis. The 101 NPC marker sequences were used as the reference to call the orthologous sequences from the contigs set of each sample by tBLASTN. Then, a reversed BLASTN was performed with the orthologous contigs obtained from the preceding steps against the reference sequences to identify potential chimeras. Following this process, the orthologous contigs set for each sample contained no more than one sequence for each marker.