Diversity of Subterranean Termites in South India Based on COI Gene

Termites (Isoptera) represent up to 95% of soil insect biomass [1,2] show an elaborated morphology and complex behaviour. The family Termitidae is considered to be the most evolved group with about 85% of all known genera and nearly 70% known species Ohkuma et al. Termites are considered as the most abundant invertebrates that are classified in about 280 genera, and over 2800 species within 14 subfamilies [1,3-5]. In India about 300 species within seven families have been reported [6] The subterranean termites which are of economic importance Wang et al. [7] to agriculture are mostly found in temperate climates Groves et al. [8] Termites are often referred as “ecosystem engineers” [9,10] as they play a vital role in recycling of plant materials and wood, modifying and improving the soil condition and composition, and providing food for other animals [3,7,11] and are also considered as potent catalyst due to their role involved in converting lignocellulose into biofuels Manjula et al. [12] which is of industrial value.


Introduction
Termites (Isoptera) represent up to 95% of soil insect biomass [1,2] show an elaborated morphology and complex behaviour. The family Termitidae is considered to be the most evolved group with about 85% of all known genera and nearly 70% known species Ohkuma et al. Termites are considered as the most abundant invertebrates that are classified in about 280 genera, and over 2800 species within 14 subfamilies [1,[3][4][5]. In India about 300 species within seven families have been reported [6] The subterranean termites which are of economic importance Wang et al. [7] to agriculture are mostly found in temperate climates Groves et al. [8] Termites are often referred as "ecosystem engineers" [9,10] as they play a vital role in recycling of plant materials and wood, modifying and improving the soil condition and composition, and providing food for other animals [3,7,11] and are also considered as potent catalyst due to their role involved in converting lignocellulose into biofuels Manjula et al. [12] which is of industrial value.
Termites are the serious pests of agricultural and horticultural crops that mainly destroys the roots and above ground parts and feed on paper, wood and timber Murthy et al. [10]. They are found to considerably damage the artificial structures and commodities due to their large colony size, varying feeding preferences and nesting behaviour Wang et al. [7]. It has been estimated that worldwide the overall cost annually for the control of damage caused by termites is more than $20 billion NY [13].
The understanding of termites in various biological processes is extremely limited regardless of their importance in agriculture and hence it is very important to understand their biology and ecology, which greatly relies on accurate species identification Singla et al. [14]. Termite on systematics was exclusively based on taxonomical identification based morphological characters of individuals belonging to various castes (e.g., soldiers or workers) Kambhampati and Eggleton [3]. Nevertheless, the caste differentiation, eusocial behaviour, varying physiological functions and crypto-biotic structure had contributed to ambiguities in their morphological identification [7,10] and species diagnosis has become a challenging task Kirton et al. [15].
The use of molecular methods which are fast and reliable, complementary to the morphological identification [2,16], are helpful in estimating evolutionary relatedness between the species Singla et al. [14]. Studies on mitochondrial genome sequences such as the ATrich region, 16S rDNA and cytochrome oxidase genes have shown an efficient alternative for species identification and phylogenetic studies [2,17,16,[19][20][21][22][23]35]. The mitochondrial DNA is more abundant as the mitochondrial genes evolve more rapidly, than the nuclear genome Wang et al. [7], therefore at species level mitochondrial DNA is more suitable Masters et al. [24] and various other regions also can be sequenced [8,25,26].
Cytochrome c oxidase subunit I gene (CO1) is one of the three mitochondrial DNA (mtDNA) encoded subunits of respiratory complex IV. It is a key enzyme in aerobic metabolism. CO1 gene is the most conservative protein-coding gene in the mitochondrial genome. Since, the mutation rate in CO 1 gene is fast enough it can differentiate precisely the closely related termite species and assess their phylogeny, in understanding the evolutionary relationships, the gene has been extensively used to understand the diversity and genetic relatedness.
In the present study, we characterized twelve species of termites of family termitidae obtained from various locations in South India based on the mitochondrial CO 1 gene and studied their phylogenetic relation.

Collection and identification of termite samples
Termite specimens collected from different locations in South India were preserved in absolute alcohol and stored at -80°C at the Division of Molecular Entomology, NBAIR-ICAR Bangalore, India. Taxonomical identification of these specimens was done at the Division of Entomology Indian Agricultural Research Institute, New Delhi, Institute of Wood Science Technology Bangalore, and Centre for Insect taxonomy, University of Agricultural Sciences, Bangalore.
water followed by dissection of the head region. The dissected head was dried and collected in 1.5ml eppendorf tube to which 180 µl of ATL buffer was added and homogenised using micro pestle. 20 µl of proteinase K was added to it and mixed methodically by vortexing and then incubated in water bath at 56°C overnight. 100 µl of AL buffer was added and kept further for incubation at 56°C for 10 minutes, after incubation 100 µl of 100% ethanol was added to the eppendorf tube and vortexed. The solution was then transferred into the mini spin columns with silica member that binds the genomic DNA. The columns were then centrifuged at 8000 rpm for 5 minutes and the flowthrough collected in a tube was discarded. The column was then further given two washes with two different buffers present in kit i.e., AW1 and AW2 and centrifuged at 8000 rpm for 5 minutes simultaneously. These columns were then transferred to a new fresh 1.5 µl eppendorf tubes in order to elute out the DNA bound to the Silica member in the column. The elution was carried out by pipetting 100 µl of sterile water to the columns, centrifuged at 8000 rpm for 5 minutes. The DNA was then checked at 1% agarose gel and stored in 4°C until PCR was done.

Quantification of DNA by nanodrop spectrophotometer
The extracted DNA was quantified using nanodrop spectrophotometer and the amount of DNA present in the samples was identified. Since, all the organic compounds show a characteristic absorption, the nitrogenous bases in DNA show a strong absorption at a wavelength of 260 nm.

PCR amplification of CO1 gene
The DNA obtained was then amplified for a portion of mitochondrial CO1 gene fragment, using the universal primers CO-1 F 5' GGTCAACAAATCATAAAGATATTGG 3' and CO-1 R 5' TAACTTCAGGCTGACCAAAAAATCA 3' , which were obtained from M/S Bioserve biotechnologies (India) Pvt Ltd. Each PCR reaction mixture of 25 µl consisted 2.5 µl of 10x PCR buffer with 15 mM MgCl 2 , 2.0 µl of dNTP's mix, 1 µl of each forward and reverse primer, 1 µl of Taq Polymerase (1U/µl), 2.5 µl of template DNA and 15 µl of sterile water. The PCR was carried out in thermal-cycler (BioRad, USA) and the conditions set were initial denaturation at 950 C for 4 minutes, denaturation at 950 C for 30 seconds, annealing and extension at 500ºC and 720ºC for 1 minute, respectively and final extension for 7 minutes at 720ºC. The PCR was carried out for a total of 34 cycles. Amplified DNA was then checked on 1% agarose gel using a DNA ladder of 250 bp and the gel was visualized in gel dock ( Figure 1).

Sequencing
Dideoxy method or chain termination method was used for sequencing of DNA, wherein, the modified bases dideoxy bases i.e., ddNTP's were used. The sequencing of amplified CO1 product was carried out at M/S. Eurofins Pvt Ltd, Bangalore. The sequence data was retrieved in the form of chromatograms which was then submitted to genbank for obtaining the accession numbers.

Sequence analysis and data interpretation
Chromatograms were edited in order to remove the ambiguous bases. These edited sequences were then aligned using Basic Local Alignment Search Tool (BLAST), with the sequences of same or related genera retrieved from the nucleotide database (PUBMED) of National Centre for Biological Information (NCBI). The CO1 nucleotide sequences of the termite species included in our present study were aligned and compared with the species obtained from PUBMED, using CLUSTAL W alignment Thompson et al. [27].

Phylogenetic analysis
Phylogenetic tree was constructed using neighbour-joining method and the evolutionary distances were computed using pdistance method with a boot strap consensus of 1000 replicates Tamura and Nie [28]. Constructed phylogenetic tree was visualised using tree viewer program.

Results and Discussion
The genomic DNA was collected from 12 populations of termites and the CO1 gene was characterized, using the universal primers CO-1 F 5' GGTCAACAAATCATAAAGATATTGG 3' and CO-1 R 5' TAACTTCAGGCTGACCAAAAAATCA 3' . The amplified CO1 product was sequenced at M/S. Eurofins Pvt Ltd, Bangalore. The sequence data was submitted to NCBI genbank and the accession numbers were obtained for the populations of termites (Table 1).

Nucleotide analysis
The complete gene analysis of nucleotide sequence for each of the collected termite species showed a considerably high percentage of A+T base composition content, with an average composition of A+T=54.88% and G+C=45.11% of various species ( Table 2).
The variation in A+T (%) among the different populations was 4.63% and 5.68% with respect to G+C (%). The populations from Attur, Sivaganaga, Thrissur and Marthalli had the least variation in the base composition of A+T (%) (54.17 to 54.89) and (45.28 to 45.83), respectively. The estimated Transition/Transversion bias (R) is 1.65. Substitution pattern and rates were estimated under the Tamura-Nei [28] model. The nucleotide frequencies were A=32.20%, T/U=24.96%, C=27.27%, and G=15.57%. A tree topology was automatically computed for estimating the ML values W et al. [29]. The maximum Log likelihood for this computation was -1053.298, which involved analysis of 12 nucleotide sequences.
Codon positions included were 1st+2nd+3rd+Noncoding. All positions containing gaps and missing data were eliminated. There were a total of 206 positions in the final dataset. Evolutionary analyses were conducted in MEGA6 Szalanski et al. [30] software.

Divergence and percent identity
The divergence and percent identity was calculated using MEGA 6 software based on the sequence alignment. The overall average was found to be 0.06 for the total 48 nucleotide sequences. The number of base differences per site from between sequences is shown. Standard error estimate(s) are shown above the diagonal Tamura et al. [31]. There were a total of 190 positions in the final dataset.
Odontotermes longignathus population from Attur, Dasarahalli and Marthahalli had negligible divergence with the similarity co-efficient ranging between 0.01-0.03%, indicating the proximity in geographical location and topography. However, the population from Thrissur had low divergence of 0.03% Tamura et al. [31]. Therefore the occurrence of termite populations of species similarity across the geographical barriers did not vary, notwithstanding the habitat. Population from Mysore (Hypotermes xenotermitis) and Marthahalli indicated similar trend (Figure 2). Our observations are in broad conformity with the reports by earlier workers [6,10]. The greater occurrence of Odontotermes longignathus among the collections could be due to widest niche breadth [6,31,33].

Phylogenetic analysis
A total of 48 strains were used for the phylogenetic analysis. The phylogenetic tree built using Hypotermes xenotermitis, Hypotermes makhamensis, Odontotermes longignathus, Nasutitermes octopilis, Odontotermes escherichi, Hypoteres sp were divided into two major clusters and three smaller clusters. The phylogenetic tree displayed in (Figure 3), shows that all the species belonging to Hypotermes xenotermitis were clubbed together with bootstrap score of 97, 53, 1 and 68. All the Hypotermes xenotermitis collected from Ratnagiri, India were grouped together with the bootstrap score of 97. The first major cluster was further divided into 6 sub-clusters. In all the 6 sub-clusters only the strains collected from Ratnagiri, India were clubbed together.
The Hypotermes xenotermitis strains collected from Pune, India formed two clusters. In one cluster, most of the Hypotermes xenotermitis species collected from Pune, India were grouped together with bootstrap score of 53. In another cluster, only four sequences with the ID of KT879848, KT879846, KT879847 and KT879845 were grouped together with the bootstrap score of 63. These four gene sequences of Hypotermes xenotermitis species collected from Pune showed divergent from the other gene sequences collected from the same place. Apart from Hypotermes xenotermitis strains, four Hypotermes makhamensis strains were also collected from Pune and Sivaganga. All the four sequences along with the Hypotermes xenotermitis strain collected from Marat were grouped together and formed a cluster. The two strains (KT879848 and KT879846) collected from Pune were closely related and grouped together with bootstrap score of 87, whereas the two strains collected from Sivaganga were also closely related and grouped together to form a sub-cluster with bootstrap score of 98. The Hypotermes xenotermitis strain collected from Marat region was formed a single clade in the cluster. The Odontotermes longignathus strain collected from Thrissur and Attur were grouped together with a good bootstrap score of 70. This cluster was sub-divided into two small sub-clusters, where the Odontotermes longignathus strain collected from Attur were formed a one cluster with a bootstrap score of 71 and Odontotermes longignathus collected from Thrissur formed a single clade with a bootstrap score of 98. Apart from these, the Nasutitermes octopilis, Hypoteres sp, Odontotermes escherichi and only one strain of Odontotermes

Percent identity
Divergence longignathus collected from DAST were formed a single clade which was appeared to be very divergent from other strains used to construct the phylogenetic tree. Overall, the phylogenetic tree analysis suggested that the Hypotermes xenotermitis strains collected from Ratnagiri were appeared to be closely related with the Hypotermes xenotermitis strains collected from Pune, India.
The occurrence of different species of termites could be attributed