A novel method of caenophidian snake sex identification using molecular markers based on two gametologous genes

Abstract Sex identification provides important information for ecological and evolutionary studies, as well as benefiting snake conservation management. Traditional methods such as cloacal probing or cloacal popping are counterproductive for sex identification concerning very small species, resulting in difficulties in the management of their breeding programs. In this study, the nucleotide sequences of gametologous genes (CTNNB1 and WAC genes) were used for the development of molecular sexing markers in caenophidian snakes. Two candidate markers were developed with the two primer sets, and successfully amplified by a single band on the agarose gel in male (ZZ) and two bands, differing in fragment sizes, in female (ZW) of 16 caenophidian snakes for CTNNB1 and 12 caenophidian snakes for WAC. Another candidate marker was developed with the primer set to amplify the specific sequence for CTNNB1W homolog, and the PCR products were successfully obtained in a female‐specific 250‐bp DNA bands. The three candidate PCR sexing markers provide a simple sex identification method based on the amplification of gametologous genes, and they can be used to facilitate effective caenophidian snake conservation and management programs.


| INTRODUCTION
Sex identification is important for mating systems and sexual behavior in ecological and conservation research in vertebrates. Intriguingly, vertebrates display considerable diversity in their sex determination systems, especially in squamate reptiles which have both temperaturedependent sex determination (TSD) and genotypic sex determination (GSD) with ZZ/ZW-type, XX/XY-type, or multiple (X 1 X 2 Y and Z 1 Z 2 W) sex chromosomes (Ezaz, Srikulnath, & Graves, 2017). Moreover, sex identification cannot be reliably determined in many squamate reptiles, because males and females have similar morphology, and morphological sexual dimorphism appears in the short period of time before mating (Frýdlová et al., 2011;Garland, 1985;Wangkulangkul, Thirakhupt, & Voris, 2005). Snakes (Serpentes) are a species-rich lineage of extant reptiles that exhibit phenotypically diverse radiation (Castoe, Jiang, Gu, Wang, & Pollock, 2008;Castoe et al., 2009;Secor & Diamond, 1998), with distribution broadly found in arboreal, terrestrial, and aquatic habitats. The snakes are considered as excellent model organisms for biomedical research, and snake venom is extracted for developing antivenoms to treat snakebites or for chemotherapeutical development (Blackburn, 2006;Kerkkamp et al., 2016;Ratanabanangkoon et al., 2016). Two distinct groups of snakes are classified as follows: (i) Scolecophidia including "blind" snakes, and (ii) Alethinophidia comprising Henophidia (pythons, boas, and other "primitive" snakes) and Caenophidia (advanced snakes) (Pyron, Burbrink, & Wiens, 2013;Wiens et al., 2012). Snake ecology, population, and diversity have been widely studied (Laopichienpong et al., , 2017Liu et al., 2015;Supikamolseni et al., 2015). Accurate sex identification is important for snake management and breeding. Improvements in conservation programs are necessary to identify the sex of juveniles before the development of primary and secondary sexual characteristics to reduce the risk of extinction, or when samples are obtained without handling individuals (e.g., noninvasive sampling). This is a very important feature when working with endangered species (Waits & Paetkau, 2005). Sex identification of snakes is commonly conducted by observation of sexually dimorphic characters such as size between sexes, with males usually having larger body, length of tail, and body color as basal practical procedure (Laszlo, 1975). However, males and females have similar morphologies in several snake species such as Acrochordus spp., thereby making it difficult to identify sexes in snakes (Wangkulangkul et al., 2005). This has led to the development of alternative sexing methods based on observation of the sex organ (hemipenes). Several methods have been developed to observe the sex organs, such as cloacal probing or cloacal popping (Laszlo, 1975). However, these methods induce stress with potential injury and require specific skills. There is a crucial need for the development of a rapid, safe, and accurate method for sex identification in snakes.
Molecular and chromosomal processes underlying the sex determination of vertebrates have been extensively investigated, especially in mammals and birds (Ezaz et al., 2017). The homologous gene located in the nonrecombining region of differentiated sex chromosomes is known as a "gametologous gene." In mammals, these are the ZFX and ZFY genes and in birds they are the CHDZ and CHDW, and ATP5A1Z and ATP5A1W genes which can be applied to determine the sex of individuals (Carmichael, Fridolfsson, Halverson, & Ellegren, 2000;Ellegren & Fridolfsson, 1997;García-Moreno & Mindell, 2000;Lawson & Hewitt, 2002). Almost all caenophidian snakes exhibit GSD with ZZ/ZW-type sex chromosomes. Z sex chromosomes which have homology with chicken chromosomes 2 and 27 is the fourth or fifth largest metacentric chromosome in the karyotypes of most snake species (Matsubara et al., 2006(Matsubara et al., , 2012Rovatsos, Vukić, Lymberakis, & Kratochvíl, 2015;Vicoso, Emerson, Zektser, Mahajan, & Bachtrog, 2013). Recently, sequence analysis of two gametologous genes as the catenin (cadherin-associated protein) beta 1 (CTNNB1) and the WW domain containing adaptor with coiled-coil (WAC) has provided information on the evolutionary process of sex chromosome differentiation in snakes (Laopichienpong et al., 2017;Matsubara, Nishida, Matsuda, & Kumazawa, 2016;Matsubara et al., 2006). Size and sequence differences between the Z and W homologs of these genes were found in many caenophidian snakes, leading us to develop molecular sexing markers to identify male and female individuals. In this study, we developed novel molecular sexing markers using three primer pairs to amplify fragments from the Z and W homologs of the genes with a clear size difference of PCR products. This assay provides a rapid and reliable method to identify genetic sex across different caenophidian snake species.

| Specimen and DNA extraction
Twenty-two snake species containing both male and female individuals were examined, and detailed information was presented in Table 1. The sex of each species was identified morphologically and confirmed by mating observations and sexing probes that searched for the male hemipenes. Blood samples were collected from the ventral tail vein using a 25-gauge needle attached to a 1-ml disposable syringe containing 10 mm ethylenediaminetetraacetic acid (EDTA).
Whole genomic DNA was extracted following the standard saltingout protocol as described previously (Supikamolseni et al., 2015) and used as templates for polymerase chain reaction (PCR). DNA quality and concentration were determined using 1% agarose gel electrophoresis and spectrophotometric analysis. Animal care and all experimental procedures were approved by the Animal Experiment Committee, Kasetsart University, Thailand (approval no. ACKU00359) and conducted according to the Regulations on Animal Experiments at Kasetsart University.

| Molecular sexing marker development
For WAC genes, partial DNA fragments of exons 9-10 were amplified using target-specific primers Eq-WAC-int9-F: 5′-CTCAGCCATCTAAT CAGTCCCCAA-3′ and Eq-WAC-int9-R: 5′-GAACGCTGAAGACTTCGA GGAG-3′ (Matsubara et al., 2016). For the CTNNB1 gene, partial DNA fragments were amplified using PCR primers (Eq-CTNNB1-11-F1: 5′-AGAGACGTCCACAATCGGATTG-3′ and Eq-CTNNB1-13-R: 5′-CA GACGTTTCTTATAATCTTGTGG-3′) (Laopichienpong et al., 2017;Matsubara et al., 2016). The female-specific PCR products with different size from the DNA fragments commonly found for both males and females (WACZ and CTNNB1Z) were expected to derive from the genes WACW and CTNNB1W (File S1). In addition, a new female-specific primer, CTNNB1W-F: 5′-GACAAAGAAGCAGCTGAGTCAG-3′, was designed, based on the nucleotide difference between the CTNNB1Z and CTNNB1W sequences at exon 11 identified in our previous study (Laopichienpong et al., 2017), to amplify a female-specific PCR product using CTNNB1W-F and Eq-CTNNB1-13-R (File S1). The BDNF (brain derived neurotrophic factor) gene was used as a positive PCR control marker using primers BDNF-F (5′-GACCATCCTTTTCCTGA CTATGGTTATTTCATACTT-3′) and BDNF-R (5′-CTATCTTCCCCT TTTAATGGTCAGTGTACAAAC-3′) (Leaché & McGuire, 2006). PCR amplification was performed using 20 μl of 1× ExTaq buffer containing 1.5 mm MgCl 2 , 0.2 mm dNTPs, 5.0 μm the primers, and 0.25 U of TaKaRa Ex Taq (TaKaRa Bio, Otsu, Japan), and 25 ng of genomic DNA. PCR conditions was as follows: an initial denaturation at 94°C for 3 min, followed by 35 cycles of 94°C for 30 s, 55°C for 30 s, and 72°C for 1 min, and a final extension at 72°C for 7 min. The nucleotide sequences of the DNA fragments derived from the PCR reaction of the primer sets Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R or Eq-WAC-int9-F and Eq-WAC-int9-R were only determined directly using the DNA sequencing services of First Base Laboratories Sdn Bhd (Seri Kembangan, Selangor, Malaysia). Nucleotide sequences of the Z and W homologs from each species were searched for homologies with the nucleotide sequences in the National Center for Biotechnology Information (NCBI) database to identify DNA fragments of the target gene, using the BLASTx and BLASTn programs. All the sequences were then deposited in the DNA Data Bank of Japan (DDBJ) ( Table 1).

| RESULTS
By agarose gel electrophoresis with the primer Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R, one same-sized DNA fragment was observed in both males and females, and one additional DNA fragment was detected from females in each snake species: Daboia siamensis,  (Table 1). Individual size variation within the same species was also found in the Z-derived fragments in two species: G. oxycephalum and C. flavolineatus (Table 1). In addition to the results from PCR with Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R primers, the primer pairs Eq-WAC-int9-F and Eq-WAC-int9-R yielded one similar-sized DNA fragment for both males and females, and one additional DNA fragment for females in each snake species: D. siamensis, E. enhydris, L. madagascariensis, O. fasciolatus, A. prasina, B. dendrophila, G. oxycephalum, C. flavolineatus, C. radiatus, P. mucosa, and P. guttatus, but no female-specific DNA fragments were observed in C. ruffus, E. maurus, X. unicolor, P. bivittatus, P. regius, A. javanicus, N. kaouthia, N. siamensis, O. hannah, B. candidus, and B. flaviceps. A wide size variation between species was found for both the Zand W-derived fragments. The Z-derived fragments ranged from approximately 600 bp in E. maurus to 3,500 bp in A. javanicus, and the W-derived fragments range from 800 bp in L. madagascariensis, A. prasina, B. dendrophila, C. flavolineatus, C. radiatus, and P. guttatus to 2,200 bp in D. siamensis (Table 1). Individual size variation within the same species was also found in Z-derived fragments in one species: A. javanicus.
Additionally, in PCR with the CTNNB1W-specific primer set, amplification of DNA fragments was found only in females as a single 250-bp DNA band (Figures 1 and 2). To distinguish the absences of PCR products in males from failures of PCR reactions, the BDNF primers were used separately as an internal control under the same PCR condition.
The BDNF primers produced larger fragments than the CTNNB1Wspecific products. Samples were identified as female if two products of 250-bp (CTNNB1W) and 750-bp sized bands (BDNF) were observed and male if only the control band was amplified (Figure 2; File S5). If the control locus failed to amplify, then sex was not assigned.

| DISCUSSION
Sex identification is very important, not only for the basic understanding of the ecology and behavior of endangered or protected animals but also for establishing management and conservation plans (Dubey et al., 2011;Webb, Brook, & Shine, 2002). To contribute to breeding programs, the sexes of the snakes must be precisely identified while avoiding snake injury and stress. Traditional methods of cloacal probing or cloacal popping are counterproductive in very small species, whereas PCR has the advantage of identifying sexually heteromorphic PCR products as genomic DNA from small quantities of tissue such as blood in this study, or applied for skin remnants, slough, or eggshell membranes left behind after hatching as shown in birds (Martín-Gálvez et al., 2011). Based on the hypothesis of sex chromosome differentiation, the cessation of recombination between sex chromosomes leads to the accumulation of gene mutations, resulting in the occurrence of sexually antagonistic alleles or the functional inactivation of genes, followed by the partial deletion of the sex chromosomes (Ezaz et al., 2017). In snake lineages, Boidae and Pythonidae have morphologically homomorphic Z and W chromosomes (Matsubara et al., 2006;Olmo & Signorino, 2005). By contrast, the W chromosomes are highly degenerated in caenophidian snakes that predominantly showed heteromorphic ZW sex chromosomes (Beçak & Beçak, 1969 "−" No DNA band in both male and female. "+" One DNA band in female only. "++" One DNA band in both male and female. "+++" One DNA band in male and two DNA bands in female.  Nazareth, 1964;Matsubara et al., 2006Matsubara et al., , 2016O'Meally et al., 2010;Oguiura, Collares, Furtado, Ferrarezzi, & Suzuki, 2009;Ray-Chaudhuri, Singh, & Sharma, 1971;Singh, 1972;Vicoso et al., 2013).

Accession number
Considering the practically routine technique, cytogenetic approaches to examine the sex chromosomes require large sample volume with long cell culture and chromosome preparation time. Therefore, they are not practical in terms of wildlife ecological studies and conservation programs. A molecular sexing method utilizing sex-specific sequences is, thus, more advantageous than cytogenetic analyses to identify sex chromosome systems. We developed novel PCR-based molecular sexing methods with three primer sets to identify individual caenophidian snake sex, based on the nucleotide sequence differences of two gametologous genes. In molecular sexing with two of the three sets, Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R, and Eq-WAC-int9-F and Eq-WAC-int9-R, the males with the homogametic sex chromosome (ZZ) were characterized by a single DNA fragment band from the two Z homologs, and the females with the heterogametic sex chromosome (ZW) were identified by two bands differing in fragment sizes from the one Z and one W homologs. The two primer sets Eq-CTNNB1-11-F1-Eq-CTNNB1-13-R and Eq-WAC-int9-F-Eq-WAC-int9-R were available for molecular sexing in 16 and 12 caenophidian snakes, respectively. These two markers exhibited co-dominant DNA pattern type. This suggests that the Z and W forms of the CTNNB1 or WAC genes were differentiated by the cessation of recombination in these caenophidian lineages, which led to insertions and deletions of nucleotide sequences on the Z and W chromosomes.
No female-specific PCR products of W homolog were found in A. javanicus for the CTNNB1 gene with the primer set Eq-CTNNB1-11-F1 and Eq-CTNNB1- 13-R and N. kaouthia, N. siamensis, O. hannah, B. candidus, and B. flaviceps for the WAC gene with the primer set Eq-WAC-int9-F and Eq-WAC-int9-R, even though these species belong to Caenophidia. This suggests that the absence of the size differences with female-specific PCR products was caused by the accumulation of mutations at existing primer sites, leading to failure of PCR reaction at the W homologs. Moreover, large insertions in the amplified target fragments could result in failure of PCR amplification for the extended W sequences. Alternatively, CTNNB1 and WAC may have been lost independently from the W chromosomes in the species.
This result agreed with the chromosome and genomic map of snakes which showed the absence of these two genes in some snake lineages (Matsubara et al., 2006;Vicoso et al., 2013). The third explanation is that there is none or little difference in sizes of PCR products between the Z and W homologs, as a consequence of a very small region of the nonrecombining portion of the W chromosome which did not have time to diverge significantly from the Z chromosome in these genes.
Snake W sex chromosomes are known to degenerate at varying rates and undergo substantial reorganization over short periods of evolutionary time (Oguiura et al., 2009). The observation of DNA bands of F I G U R E 1 Agarose gel electrophoresis of PCR products in males and females of seven snake species using Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R (a), Eq-WAC-int9-F and Eq-WAC-int9-R (b), and CTNNB1W-F and Eq-CTNNB1-13-R (c  Figure 1). Such individual variation in the Z form of the gametologous gene within the same species was also reported by nucleotide analysis for the CHD1Z gene of birds (Casey, Jones, Sandercock, & Wisely, 2009;Friesen, Congdon, Walsh, & Birt, 1997;Trimbos et al., 2013) and for the CTNNB1 gene of snakes (Laopichienpong et al., 2017).
To avoid the misidentification of sexes through size variation, we designed the additional primer at the exon of the CTNNB1W gene which yielded a female-specific 250-bp PCR products in 16 caenophidian snakes and indicated dominant DNA pattern type.
Internal PCR control with the BDNF gene was also used to avoid amplification failure which might result in a misidentification of a female as a male. The present results of the three sexing markers could also be applied to identify the sex of individuals in new snake species as a simple procedure for males and females. We highlight that all the three sexing markers, with an additional BDNF control marker, could be simultaneously used to identify sex with higher diagnostic accuracy.
In the two henophidian snakes P. bivittatus and B. constrictor, CTNNB1 and WAC genes were located on both Z and W chromosomes, and no sex-specific sequence was detected for the two gametologous F I G U R E 2 Phylogenetic relationships among sampled snake species illustrating the sex-specific amplification of CTNNB1 and WAC genes using the primers: Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R; Eq-WAC-int9-F and Eq-WAC-int9-R; and CTNNB1W-F and Eq-CTNNB1-13-R. Phylogeny was partially derived from Vidal, Rage, Couloux, and Hedges (2009). Agarose gel electrophoresis of PCR products in males and females of twenty-two snake species using three sexing markers: Eq-CTNNB1-11-F1 and Eq-CTNNB1-13-R (column A), Eq-WAC-int9-F and Eq-WAC-int9-R (column B), and CTNNB1W-F and Eq-CTNNB1-13-R (column C) are indicated at the right edge of the tree. M, male; F, female genes (Matsubara et al., 2006(Matsubara et al., , 2016Vicoso et al., 2013). In this study, five henophidian snakes (C. ruffus, E. maurus, X. unicolor, P. bivittatus, and P. regius) showed the absence of female-specific PCR products in molecular sexing with our three primer sets. This implies that the CTNNB1 and WAC genes on the Z and W chromosomes have not been differentiated in henophidian snakes (Laopichienpong et al., 2017;Matsubara et al., 2006Matsubara et al., , 2016Vicoso et al., 2013 This study reveals the potential and usefulness of gametologous sequences which develop new candidate sexing markers in caenophidian snakes. Our work opens a new tool kit for further application with the genomic DNA of small and degraded tissues such as slough as noninvasive sampling. To facilitate application of this tool, we provide a "know-how" guide to apply Loop-mediated lsothermal Amplification (LAMP) or aptamer to conservation program where species identification or sex determination is needed in real time and in situ. Future work in this direction will enormously facilitate the work of wildlife managers, researchers, and exotic snake breeders with important conservationist, economic, and commercial benefits.

ACKNOWLEDGMENTS
This study was financially supported by grants from Professor