Phylogenetics and Biogeography of Cobra (Squamata: Naja) in Java, Sumatra, and Other Asian Region

The separation of Sumatra and Java occurred at the end of the Miocene (10 mya) to the end of Pliocene (1.8 mya). The existence of ecological variations and geographic barriers inhibits gene flow through the isolation of adaptation, geography, reproduction, inbreeding, and leading to population segregation. Cobra (genus Naja) distribution became greatly influenced by the geologic condition and sea level. This study was conducted by phylogenetic analysis towards the 16S rRNA gene. Survey was done with Visual Encounter Survey (VES) method on 6 locality. There were 3 location in Sumatra Island and the others in Java Island. Sample from other Asian region was obtained from Genbank, which are 11 individuals from China, Thailand, and Nepal. DNA extraction was done according to the QIAmp® DNA Mini Kit standard protocol. The forward and reverse 16S sequences are combined with the SequencherTM version 4.1.4 program, then in BLAST (Blast Local Alignment Search Tool) at www.ncbi.nlm.nih.gov. Phylogenetic analyzes of clade A (MP = 60, ML = 54, BI = 88) indicate the presence of division into two monophyletic subclade (AI and AII). Subclade AI consists of groups of Cobra from Sunda (Thailand, Sumatra and Java). Subclade AII is a group of species N. kaouthia originating from Chumpon Province, Southern Thailand with (MP = 96, ML = 95, BI = 100). Clade B is divided into two subclasses (BI and BII). The result is supported by bootstrap value MP = 93, ML = 99, BI = 100. N. atra of Fujian Province is a sister lineage of the same species from Jiangxi Province (MP = 86, ML = 86, BI = 100).


INTRODUCTION 
The evolution of South Asia's diversity was occurred in Paleocene epoch (60 mya-million years ago), when the Malay Peninsula, Java, Sumatra, and Borneo still merged as a big continent, called as sundaland. Then, the formation of the volcanic rings during Eocene (40 mya) and sea level fluctuations at the beginning of Miocene (23.03 mya) initiated the separation between the Malay Peninsula and Borneo. On the late Miocene, sea level fluctuations occured again and half of Java Island were flooded by water. The separation of Sumatra and Java as well as Borneo rotational changes occurred at the end of the Miocene (10 mya) to the end of Pliocene (1.8 mya) [1]. The movement of the earth's plate during the land evolution was separate several regions of Asia. In addition, Sumatra Island is estimated to occured in last 25 million years ago (late Miocene) and divided into west and east regions by Bukit Barisan mountains. The existence of ecological variations and geographic barriers inhibits gene flow through the isolation of adaptation, geography, reproduction, inbreeding, and leading to population segregation [2]. Cobra (genus Naja) distribution became greatly influenced by the geologic condition and sea level [3]. This study was conducted by phylogenetic analysis towards the 16S rRNA gene. The purpose of this study is to estimate the relationship of Naja genus in the island of Java, Sumatra, and other Asian region based on the 16S rRNA mithocondrial gene. Zoogeography analysis is also examined for more advance understanding between the correlation and evolution of Sundaland.

MATERIALS AND METHODS Sampling and Preservation
Survey was done with Visual Encounter Survey (VES) method on 6 locality. There were 3 location in Sumatra Island and the others in Java Island (Table 1). Sample from other Asian region was obtained from Genbank, which are 11 individuals from China, Thailand, and Nepal (

DNA Extraction
DNA extraction was done according to the QIAmp ® DNA Mini Kit standard protocol. A total of 25 mg of tissue was cut in tiny pieces and placed in 1.5 ml microtube. The tissue was given 180 µL of ATL buffer and 20 µL of K-proteinase. Samples were incubated in waterbath-shaker, and was done in 56 o C for 1-3 hours, then added 200 µL of AL buffer. Sample incubation at a temperature of 56 o C was done for 10 minutes, then added 200 µL of 90-100% ethanol. Each procedure is terminated with vortex. Samples were pipetted and inserted into 2 mL collection tube and placed in Dneasy Mini Spin. Centrifugation was done at 8000 rpm (one minute), supernatant is then discarded. Spin Column was placed in new collection tube and was given 500 µL of AW1 buffer. Sample then centrifuged on 8000 rpm (one minute), supernatant is then discarded. Spin column then placed in new collection tube (2 mL) and was added 500 µL of AW2 buffer. Centrifugation was done at 14.000 rpm (3 minute), supernatant is then discarded. Spin column then placed in 1.5 or 2 mL microtube. DNA was thawed with 200 µL AE buffer in the middle of spin column. Incubation was done at one minute (15-25 o C), then centrifugated at 8000 rpm (1 minute).

Qualitative Analysis
Qualitative analysis was done by agarose gel (1%) electrophoresis [4]. Gel was made using 0.3 g agarose and 30 mL Tris-Boric-EDTA (TBE) with pH 8. Gel was made inside Erlenmeyer flask until boiled a little. Warm agarose was given 1 µL EtBr and then homogenized. Homogenate was poured slowly into gel cast, and allowed to harden. The hardened gel was inserted into electrophoresis chamber and soaked in TBE (pH 8). DNA sample is mixed wit loading dye, then inserted into wells as well as the DNA marker ladder 1 kb (1:1). Electrophoresis then conducted at half voltage for 45-60 minute. The results were documented with Gel Doc.

Phylogenetic Analysis
The forward and reverse 16S sequences are combined with the Sequencher TM version 4.1.4 program, then in BLAST (Blast Local Alignment Search Tool) at www.ncbi.nlm.nih.gov. As outgroup, Ophiophagus hannah (Family Elapidae) is used with DNA sequence data that is accessed from the site www.ncbi.nlm.nih.gov. Both sample sequences and outgroup species were aligned with ClustalW analysis on the MEGA6.0 program.
The reconstruction of phylogenetic tree was done by Maximum Parsimony (MP) analysis, with bootstrap 1000 in PAUP 4.0b10a program. The Jmodeltest 3.06 program is used for the best model of likelihood analysis. Analysis of MP, Maximum Likelihood, and Bayesian Inference were performed on Mr.Bayes program ver 3.1.2 MCMC 10 million generation. Haplotype Network Analysis is done with DNAsp program and Bandelt Network.

RESULT AND DISCUSSION Qualitative Analysis
The results of electrophoresis of all N. sputatrix DNA samples and N. sumatrana showed the presence of DNA with an ampicton size of 500-600 bp (Fig. 1). An amplycone length range of 550-590 bp [6] is the result obtained from the amplification of 16S rRNA gene with 16 SBR primary and 16 SAR as reverse and forward respectively.Visualization of DNA bands describes the purity of DNA resulting from isolation process, hence the visualized DNA can be used as a further molecular analysis material [6].

Genetic Distance of Naja
The determination of new species on N. haje is indicated by the percentage of p-Distance value of the 16S rRNA gene ≥5% [7]. The genetic characteristics of a species are influenced by differences in geographical conditions, while geographical barriers can alter the genetic distance of a population [8]. The p-Distance results showed a low value (1-2%) between the species of N. sputatrix (Java) and N. sumatrana (Sumatra) . The obtained value is supported with a similarity rate of 96-97% in the BLAST analysis. The p-Distance value between N. sumatrana from Sumatra and Thailand is low (p-Distance 1%), which is supported 97-98% similarity level on BLAST result. The p-Distance value between N. atra (China), N. sputatrix (Java), and N. sumatrana (Sumatra) is high (p-Distance 4-16%). The p-Distance value between N. kaothia and N. sputatrix (Java) and N. sumatrana (Sumatra) is low (p-Distance 2-3%), while p-Distance value between N. kaouthia with N.atra and N. naja classified as high (p-Distance 9-12%). The p-Distance analysis concluded that N. sputatrix, N. sumatrana, and N. kaouthia were spread in the Sundaland and had close phylogenetic relationships. Based on the p-Distance value of ≥5% for the determination of new species, the three Naja species may be possible still within the same species (Table 3).

Phylogenetic Analysis
The phylogenetic analysis of Naja genus based on 16S rRNA (MP=98, ML=80, BI=100) shows the formation of two clades of Naja (A and B), and one outgroup, which is Ophiophagus hannah (C) (Fig. 2). The results of this study led to the idea to review the results of research by Philip [9] which states that Ophiophagus hannah is a sister lineage of the genus Naja.    According to BLAST results, the AIa group (N. sputatrix from Java) is similar to Naja sumatrana and Naja kouthia with 96-97% similarity level. Wuster and Thrope [10] show that N. kouthia and N. sumatrana in Asia are one group with N. sputatrix and represent a lineage with appropriately similar morphology or distribution and ecological characteristics. The results showed that N. sputatrix from Malang and Bondowoso (Fig. 3) had closer relationship than N. sputatrix from Cilacap, in accordance with similar morphological characters (brown scales on N. sputatrix from Malang and Bondowoso, while black on N. sputatrix from Cilacap).
Das [3] revealed that N. sputatrix has 1189 mm (SVL) and155 mm (TL). Meristic characters are two nasal scales, no loreal scales, one supraocular scales, one preocular scoop, three postocular scales, supralabial scales, nine infralabial scales, 19 dorsal central scales, smooth surfaces, 180 ventral scales, and anal scales are divided. Generally, there is no pattern behind the hood, but if there is, it usually resembles the shape "V". The presence of small scale accessories [11] on the head scales is not directly related to the process of adaptation, but in contrast, color is possible to be one of the adaptation factor. Populations in Central Java to West Java have dominant blackish color and occupy a thick humid rainforest, while the population from East Java has a silver to brownish scales, occupying habitat with dry soil conditions. The AIb group (N. sumatrana) originating from Sumatra has a similar likelihood of 97-98% with the same species from Thailand, while it has a similarity of 95-96% with N. kaouthia (AII) of affect the distribution process of this species in Southern Thailand resulting in competition between N. sumatrana and N. kaouthia. N. sumatrana has a body length of 1.5 m with a strong body, large head differs from the neck, round muzzle, and round hood (adult) or oval (juvenile). This species has a single preocular scales, postocular scales of two to three, seven supralabial, four infralabial, 179-206 ventral scales, 40-57 subcaudal scales. The back color is related to geographical origin and size [3].
Subclade AII is a group of species N. kaouthia originating from Chumpon Province, Southern Thailand with (MP = 96, ML = 95, BI = 100). The results of the phylogenetic reconstruction trees in accordance with the results of research Wuster [12] which showed that three evolutionary lines separated from the lineage of the genus Naja. The main sub-lineage representing the Naja genus of Asia represented by N. naja where the reconstructed phylogenetic tree by Bayesian analysis shows that N. kaouthia is a sister lineage of N. atra supported by a valid matching value (bootstrap 95%).

Clade B
Clade B is divided into two subclasses (BI and BII). Sequences with valid matching values from the BI subclass N. atra are from China (Fujian, Guangdong, Jiangxi, Yunnan province) are taken from GenBank. The result is supported by bootstrap value MP = 93, ML = 99, BI = 100. N. atra of Fujian Province is a sister lineage of the same species from Jiangxi Province (MP = 86, ML = 86, BI = 100). N. atra of Guangdong Province is a sister lineage of the same species of Jiangxi Province with a strong matching value (MP = 100, ML = 100, BI = 100).
Subclade BII is N. naja from Nepal (India), and has a valid topology (MP = 100, ML = 100, BI = 100). N. naja according to the results of the reconstruction of phylogenetic trees is the species with the closest kinship with outgroup species Ophiophagus hannah. The results of the reconstruction are in accordance with research from Wuster [12] and Wallach [13], where N. naja are separated from clade of N. kaouthia, N. sputatrix, and N. Sumatrana.

Naja Biogeography
The divergence time of Elapidae occurs at the beginning of the Neogenic Period (approximately 16 mya) to the end of the Cretaceous subperiod [14]. The distribution of this family is related to the evolution of Sundaland. There are two main factors affecting the distribution of cobra in the Malay Peninsula and Indonesia, i.e. ecological conditions and sea surface changes [3].
The separation of Elapidae occurred in about 16 mya [16]. N. sputatrix is recorded only in Indonesia (Java, Bali, Lombok, Sumbawa, Flores, and possibly the surrounding islands). The abundance of N. sputatrix in Java is varies depend on locality [17]. Ecological conditions is found to be correlated with patterns of diversity among N. sputatrix, and occur in various climatic zones in Java and the Sunda Islands. The phylogenetic trees reconstruction and the value of p-Distance show that N. sputatrix in Malang (East Java) with Bondowoso (East Java) is closely related.
Separation of Java begins at late Miocene, which accompanied by sea water fluctuations that cause half of Java to be submerged by wate (Fig. 5). Separation of Java and Sumatra occurs during the end of Miocene (10 mya) to the end of Pliocene (1.8 mya) [1]. The process of land separation causes a geographical isolation on N. sputatrix. Geographical isolation ultimately limits the spread of this species due to its inability to adapt to the Sumatra, Kalimantan, and other Asian environments. This insights explain why N. sputatrix is only found in Java. Furthermore, at the beginning of the Pleistocene (250,000 years ago) there was a decrease of sea water surface, allowing the distribution of N. sumatrana in southern Thailand. This species is able to adapt to climatic and ecological conditions in Thailand and is able to spread from Sumatra and Kalimantan [18].
N. kaouthia that used in this study was came from Southern Thailand. N. kaouthia experienced separation with the divergence time of 1-0.8 mya [16]. The species is distributed almost in various regions: Myanmar, Thai Peninsula, South Laos, Cambodia, South Vietnam, northern Peninsular Malaysia, India, Nepal, Bangladesh, South China, Sichuan Province, and Yunnan [3]. N. kaouthia throughout the mentioned region is the result of the evolution of the Sundaland. The Distribution was began in the early Eocene to the late Oligocene (25 million years ago), when Burma, Malay Peninsula, Java and Sumatra were still in the mainland [2]. N. naja is separated from Elapidae family with divergence time of 1-0.8 mya [16]. The relation of this species with another Naja in Sundalandis quite distant. This evidence can be explained by the separation proccess of Indian and African land that occurred at the end of the Mesozoic era (120 mya) to the beginning of the Cenozoic era (65 mya). The expansion of the ocean floor that occurred during the Tertiary period led to the movement of the Indian mainland towards the Asian Continent (Fig. 5). The movement of mainland India resulted in the occurrence of collisions and continues on the formation of the Himalayas and the Tibet plateau [19]. The N. naja species subsequently experienced an allopathic specia-tion due to the separation of population by the ecological barrier (Himalayan Mountains).
N. atra has a divergence time of about 3 mya. Nanling and Luoxiao mountains do not become a barrier in the distribution of N. atra in mainland China, where in this study is represented by the provinces used as sampling data (Fujian, Yunan, Jiangxi and Guangdong). The existence of N. atra in Taiwan, an island separated from mainland China was seen as any divergence split progress. However, N. atra in Taiwan differs from the species of mainland China. N. atra in Taiwan has a distinct genetic character and shows the separation of the phylogenetic tree lineage between populations in mainland China and Taiwan. Formation of Taiwan Island occurred during the Miocene, thenit was separated from mainland from the southeast coast of mainland China by Taiwan Strait. Climate change during the Pleistocene caused most of the Taiwan Strait to become a landmass, thus providing an opportunity for N. atra to migrate from mainland China [20].