Genetic Variation of Eight Indonesian Swamp-Buffalo Populations Based on Cytochrome b Gene Marker

Genetic variation is a major concern in animal genetic resources conservation program. This study aimed to analyze genetic variation and phylogeography of eight Indonesian swamp-buffalo populations based on cytochrome b gene marker. A total of 78 DNA fragment samples originating from eight Indonesian swamp-buffalo populations were used in this study, namely Bombana Island, Bombana mainland, Kolaka, Konawe, North Toraja, West Nusa Tenggara, Banten, and Aceh with 11, 10, 13, 14, 10, 10, 5, and 5 samples, respectively. The cytochrome b gene sequence and genetic variation parameters were analyzed in MEGA software (ver 6), and DnaSP software (ver 5.10.01). The results of this study showed that all DNA-fragment samples were successfully amplified by PCR technique with the size target (906 bp). Based on the distribution of all samples, it was found 9 polymorphic sites, and 10 haplotypes with the haplotype diversities were 0.6590. The average of genetic distances between populations ranged from 0.0000-0.002. They were grouped into two main clusters. The first cluster consisted of Aceh, North Toraja, West Nusa Tenggara, Banten, Kolaka, and Konawe populations, meanwhile, the second cluster consisted of Bombana Island, Bombana mainland, Kolaka, and Konawe populations. The results of the study were concluded that eight Indonesian local swamp-buffalo populations were grouped into two main clusters where Bombana Island and Bombana mainland populations were specific breeds because they were only found in the second cluster and also had specific nucleotides sites on the 57 nucleotides which C base changed to T. The results of this study were useful in formulating the program of conservation and utilization of Indonesian buffalo genetic resources, especially in the buffalo population with specific breeds.


INTRODUCTION
The buffalo population has decreased drastically in almost all countries in Southeast Asia. By the years of 1990-2016 period, for example, the decreasing buffalo populations were successively for the highest in Thailand (82.59%), followed by Indonesia (58.43%), and Malaysia (42.19%), meanwhile the Philippines has slightly increased (3.58%) (FAOSTAT 2015;FAOSTAT 2018). The total population and productivity of buffaloes in Indonesia were much lower than that of cattle. In 2018, the local buffalo population in Indonesia was 1,356,390 heads; meanwhile the cattle population was 17,050,006 heads. The potency of buffaloes as national meat resources was only 0.88%; meanwhile the potency of cattle was 13.81% (Ditjen of PKH, 2018). The trend of population declining is a serious threat to the sustainability of buffalo genetic resources. are more than eight months per year. Simelue buffaloes adapt to coastal and hilly areas, meanwhile Toraya buffaloes adapt to agricultural areas (Talib et al., 2014).
Mitochondrial DNA marker has been used as one approach in molecular genetic studies to identify genetic variation and phylogeny in the animal. The main reason for using mitochondrial DNA marker because the mitochondrial genome has a relatively small size (± 16500 bp), and a rapid evolutionary rate, which causes high intra-species diversity (Avise, 1994). Mammalian mitochondrial DNAs show several special features such as an absence of intron, maternal inheritance, all offsprings from the same female parent having the same mitochondrial DNA sequence, the existence of single copy of orthologous genes, the lack of recombination events and a high mutation rate (Linacre & Tobe, 2011;Kim et al., 2013). Mitochondrial DNAs of the animal contain thirteen protein-coding genes and are considered as strong markers for genetic variation analysis such as at the family, genus, and species levels (Arif & Khan, 2009;Doosti et al., 2014). The cytochrome b (Cyt b) gene is one of the protein-coding genes in mitochondrial containing abundant phylogenetic information among intra-species and inter-species and it is considered to be a good marker to study the genetic differentiation and phylogenetic relationships among species within the same genus or the same family. Cyt b gene markers have been proved as an efficient and powerful tool for breed characterization and species identification (Browers et al., 1994;Zardoya & Meyer, 1996;Saif et al., 2012).
Previous studies on genetic variation and phylogeography based on mitochondrial Cyt b gene marker in livestock have been reported, such as in goats (Pakpahan et al., 2016), Romanian cattle (Xuan et al., 2010), Korean cattle (Kim et al., 2013), Pakistani cattle (Hussain et al., 2018), Indonesian local cattle (Hartatik et al., 2015;Hartatik et al., 2018), Nili-Ravi and Kundi buffaloes in Pakistan (Saif et al., 2012), and Southeast Asian buffaloes (Lau et al., 1998;Yindee et al., 2010;Lei et al., 2011;Zhang et al., 2011;Zhang et al., 2016). Study related to genetic variation and genetic relationships of Indonesian local swamp buffaloes was reported by Sukri et al. (2014) that from eight Indonesian local swamp-buffalo populations (Madiu, Aceh, Blitar, Riau, Kalimantan, Tator, Bali, Lombok, and Bima), it was found 19.87% variable sites consisted of 16 haplotypes and 7 different haplogroups unique based on their geographical regions. Meanwhile, Amin et al. (2015) reported that four local swamp-buffalo populations in Central Indonesia (Bali, Toraja, Lombok, and Bima) had close genetic relationships, which they only consisted of one main cluster and two sub-clusters. Also, there are still many local swamp buffaloes from the other regions in Indonesia that have unique characteristic with specific geographical spread, such as swamp buffaloes from Banten, Riau, Sumenep, and Southeast Sulawesi (Sumantri et al., 2017), but their genetic variations and genetic relationships have not been identified yet.
Information about genetic variation within and between population or breed has been needed in conservation and breeding programs of buffalo. The data of genetic differences and genetic relationships between and among populations or breeds are important to determine the best approach of AnGR conservation . Therefore, this study was carried out to analyze the genetic variation and phylogeography of eight Indonesian swamp-buffalo populations based on mitochondrial Cyt b gene marker.

Sample Collection and Genomic DNA Extraction
A total of 78 blood samples of local swamp-buffalo from eight different populations were used in this study (Table 1). Blood samples (3 mL) were collected from jugular veins and put in vacutainer tubes containing EDTA, then absolute ethanol was added and the mixture was stored at 4 o C. Genom DNAs were extracted from the blood. The Genomic DNA extraction protocol was performed according to Sambrook & Russel (2001) with a minor modification.

Primer Design
The forward and reverse primers were designed using the Primer3 software (ver 4.0) (http://primer3. ut.ee) based on the Cyt b gene sequence of swamp buffalo from GenBank with acc. No. D88637. The primer forward was for: 5'-CAT TCA TTG ACC TCC CTG CT -3', and the reverse primer was for 5'-GCC GGA ACA TCA TAC TTC GT -3'.

DNA Amplification
The DNA fragment was amplified using the Polymerase Chain Reaction (PCR) technique using a thermocycler machine (GeneAmp® PCR System 9700, Applied Bio SystemTM, Foster City, CA, USA). DNA amplification was carried out at a total volume of 50 μL, containing 50 ng/μL DNA templates, 25 pmol primers (IDT, Singapore), 1 unit of Go Taq® Green Master Mix (Promega, Madison, WI, USA), and water. The PCR process was run with 30-35 cycles consisting of predenaturation at 95ºC for 5 min, denaturation at 95ºC for 10 s, annealing at 52ºC for 20 s, extension at 72ºC for 30 s, and the final extension at 72ºC for 5 min. The PCR products were visualized through 1.5% agarose

Sequence Analysis
PCR products were sequenced using the Sanger method through analysis services from the 1st Base, Malaysia. All sequence results (ABI trace files) were edited in FinchTV software (ver 1.4.0) (www.geospiza. com). The Basic Local Alignment Search Tool (BLAST) was used to identify the similarity (homology) with Cyt b gene sequences in GenBank (https://blast.ncbi. nlm.nih.gov). Sequence analysis was performed using the Multiple Sequence Alignment Analysis Technique by using the Clustal program in MEGA software (ver 6) (Tamura et al., 2013). Cyt b gene sequence from GenBank with acc. No. D88637 was used as a reference. Cyt b gene sequence that was analyzed did not include sequences of forward and reverse primers with base length was 866 bp.

Data Analysis
Variables observed included conservation sites, variable sites, transition/transversion substitutions, and genetic distances were then analyzed by MEGA software (ver 6) (Tamura et al., 2013). The phylogenetic tree was reconstructed based on the Neighbor-Joining method with 1000 bootstraps. The haplotype diversity was analyzed by DnaSP software (ver 5.10.01) (Rozas et al., 2010). The Network software (ver 5.0) (Available at http://www.Fluxux-engineering.com) was used to reconstruct a median-joining network of haplotypes.

DNA Amplification and Sequence Analysis
A total of 78 DNA fragment samples in this study were successfully amplified using PCR (Polymerase Chain Reaction) technique. The size of PCR product target was 906 bp. Bands pattern visualization of Cyt b fragment in 1.5% agarose gel electrophoresis is presented in Figure 1, and the primers annealing position is presented in Figure 2. The BLAST results showed that the similarity percentage of Cyt b gene sequences was at the range of 99-100% with query cover percentage was 100%.

Nucleotide Variation
Nucleotide variations of the Cyt b genes of eight Indonesian swamp-buffalo populations are presented in Table 2. The results of multiple alignment analysis of 78 Cyt b genes sequences of Indonesian swamp buffalo were found 857 (98.96%) conserved sites and 9 (1.04%) variable (polymorphic) sites. The highest variable sites were found 5 sites (0.8%) in Konawe and Bombana Island populations; meanwhile the lowest one was found 2 sites (0.23%) in North Toraja, NTB, and Banten populations. Nucleotides position at variable sites of swamp buffalo Cyt b genes is presented in Figure 3.
Variable sites consisted of 3 (0.35%) parsimony sites-informative and 6 (0.69%) singleton sites. A site is parsimony-informative if it contains at least two types of nucleotides, and at least two of them occur with a minimum frequency of two; meanwhile a singleton site contains at least two types of nucleotides with, at most, one occurring multiple times (Tamura et al., 2013). The highest parsimony site was found 3 sites (0.35%) in   2 3 4 4 6 6  1 5 6 8 1 7 3 7  6 2 7 2 1 7 2 1

Genetic Distance
Genetic distance was calculated based on the Kimura-2 parameter method (Tamura et al., 2013). The average genetic distances among eight Indonesian swamp-buffalo populations were at the range from 0.0000-0.0022 (Table 3). The closest genetic distance was identified among North Toraja, NTB, and Banten populations (0.0000), meanwhile the farthest was between Konawe and NAD populations (0.0022).

Nucleotide Variations
A total of 78 Cyt b genes sequences of local swamp buffaloes from 8 different populations were observed in this study. Previous studies reported that from 17 Cyt b gene sequences of Indonesian swamp buffaloes, 80.13% were found to be conserved sites and 19.87% were found to be variable sites (Sukri et al., 2014). In the present study, it was found 856 (98.96%) conservation sites and 9 (1.04%) variable sites. Cyt b genes sequences of Indonesian swamp buffaloes in this study were more conservative than those reported by Sukri et al. (2014) that found lower variable sites. The use of different Cyt b genes references from the GenBank influenced the results. Cyt b gene sequence reference used in this study was from swamp buffalo (Bubalus bubalis carabanesis); meanwhile Sukri et al. (2014) used Cyt b genes sequences reference from African buffalo (Syncerus caffer).
The variable sites found in this study consists of 3 (0.35%) parsimony-informative sites and 3 (0.35%) singleton sites. The parsimony-informative site was only found in Kolaka (0.12%) and Konawe populations (0.35%), meanwhile singleton site was found in all populations. The nucleotide variations of Cyt b gene of local swamp buffaloes in this study were lower than that  of cattle. In Romanian cattle populations, 95 polymorphic sites were detected consisted of 73 parsimony-informative sites and 21 singleton sites (Xuan et al., 2010), meanwhile in Ethiopian cattle populations, 18 polymorphic sites were detected consisting of seven parsimonyinformative sites and 11 singleton sites (Tarekegn et al., 2018). The positions of the variable site were found in the 6,12,57,262,831,417,472,631, and 678 nucleotides ( Figure 4). In the 57 nucleotide, C base changed to T. This mutation (g.57C>T) was only found in the Cyt b gene sequences of swamp buffaloes from Southeast Sulawesi, especially Bombana Island and Bombana mainland populations. This nucleotide mutation was found in all swamp buffaloes from the Bombana Island and Bombana mainland populations. The site of g.57C>T could be claimed as a specific nucleotide site (geographic origin) to the Island Bombana and main-land Bombana populations. The specific nucleotide site with the same nucleotide position (g.57C>T) was also found in several other swamp buffaloes from Southeast Sulawesi, consisted of Kolaka (6 individuals) and Konawe (10 individuals) populations. The emergence of the specific nucleotides sites was an influence of forming a new population (founder effect) as a result of genetic drift which was a situation where a small number of individuals move to a new place and to form a new population. The founder effect was an example of the population bottleneck which was a situation where the number of parents in the population becomes very small for one generation or more (Nicholas, 2010).
The nucleotide variations could be due to the mutations of nucleotide bases. The results of previous studies reported that there were 4 transition substitutions (C>T) in Cyt b gene sequence of swamp buffaloes in Asian mainland (Thailand and Peninsular Malaysia), but transversion substitution was not found (Lau et al., 1998). Also, Kim et al. (2013) found fifteen transition substitutions in the Korean cattle from fifteen polymorphic sites identified. In contrast, in the Pakistani river buffaloes (Nili-Ravi, Kundi) there were 79 transition substitutions (g.33A>G and g.46C>T) and 43 transver-  Bbc  BK  BD  KL  KN  TR  NTB  BTN  NAD  1  1  --------1  1.27  2  -9  10  6  2  ----27  34.18  3  -1  ---- . T sion substitutions with a ratio of 1.83: 1 (Saif et al., 2012). In this study, 7 (0.81%) transitional substitution and 2 (0.23%) transversion substitution were found in Cyt b gene sequence of Indonesian swamp buffaloes with a ratio of 6.67: 1. The nucleotide substitutions in this study were silent mutations because it only occurred in the third codon translating the same amino acid even though it has different DNA sequences (Saif et al., 2012).

Phylogeography
In this study, there were identified two major clusters from eight Indonesian swamp-buffalo populations based on their topology phylogenetic trees ( Figure  4). This indicates that the Indonesian swamp-buffalo populations have two different maternal lineages. The formation of these lineages was possibly related to the entry path of swamp buffaloes from their domestication center. Molecular genetic studies confirmed that swamp buffaloes were surmised of becoming domesticated in South China (Lei et al., 2007) or in the border regions between South China and North Vietnam (Wang et al., 2017) or Indo-Chinese and North Thailand Southeast Asia (Lau et al., 1998;Zhang et al., 2011). Swamp buffaloes from the domestication centers migrated to Indonesia through two routes. The first route was through Vietnam, Thailand, the Malay Peninsula, Sumatera, Java, Sulawesi, and Nusa Tenggara; meanwhile the second route was through China, Taiwan, Philippines, and Kalimantan Wang et al., 2017).
The phylogenetic tree (Figure 4) was shown that the first cluster (lineage I) consisted of NAD, Banten, NTB, and North Toraja populations and also a portion of Kolaka and Konawe populations. This indicated that they came from the same maternal lineage. Meanwhile, all of individual swamp buffaloes in the second cluster (lineage II) originated from Southeast Sulawesi consisted of Bombana Island, Bombana mainland, Kolaka, and Konawe populations. This indicated that swamp buffaloes spreading in various regions of Southeast Sulawesi originated from the same maternal lineage and they were different from the other swamp-buffaloe populations in Indonesia. The results of previous studies also found two different maternal lineages in Chinese native swamp buffalo based on DNA Cyt b gene and D-loop of mitochondrial (Lei et al., 2011;Yue et al., 2013), and swamp buffalo in Southeast Asia based on mitochondrial D-loop (Lei et al., 2007).
The highest genetic similarity (100%) was found between Banten, NTB, and North Toraja populations with a genetic similarity were 100%; meanwhile the furthest one was between Konawe and NAD populations, which they had the farthest geographical distance from all populations. The genetic distances among Indonesian swamp-buffaloe populations in this study were relatively closer compared to the swamp-buffalo populations in Southeast Asia based on D-loop sequences of mitochon-drial DNA at the ranges of -0.0041 to 0.0486 (Lau et al., 1998). Closer genetic distance (0.0000-0.0005) based on Cyt b gene sequence of mitochondrial DNA was detected in wild Java Banteng populations in the four different national parks (Qiptiyah et al., 2019). Meanwhile, further genetic distances based on Cyt b gene sequence of mitochondrial DNA was detected in Ethiopian native cattle populations, at the ranges of 0.00010 to 0.12446 (Tarekegn et al., 2018).
The genetic relationships among Indonesian swamp-buffaloe populations in this study were not correlated with the geographical distance of their habitats. For instance, geographical distance between North Toraja and Banten populations was further than that of geographical distance between North Toraja and Kolaka populations, but genetic distance between North Toraja and Banten populations was closer than that of between North Toraja and Kolaka populations, although swamp-buffalo populations of North Toraja, Kolaka, and Konawe originated from the same Island (Sulawesi Island). That a number of local swamp buffaloes may have a fairly close genetic relationships, and others have far genetic relationships. Genetic relationships level was often found being uncorrelated with the geographical distance of buffalo habitat. This result might be related to the fact that there were some populations that were more open because transportation routes were easily accessible so that it made easier for buffalo migration (Anggraeni et al., 2011).
Haplotypes and nucleotides diversities are alternative measures of mitochondrial DNA variations (Yindee et al., 2010). In this study, 10 haplotypes were found with haplotypes diversity of 0.6590 from 79 Cyt b genes sequences of swamp buffalo (including buffalo from GenBank). The number of haplotypes of Cyt b gene sequences of eight Indonesian swamp-buffaloe populations in this study was higher than that reported by Saputra et al. (2013) that they only found three haplotypes of seven Indonesian swamp-buffaloe populations based on the 21 Cytochrome Oxidase subunit I (COI) gene sequences. The higher haplotypes number of Cyt b gene sequences of Indonesian swamp buffalo was reported by Sukri et al. (2014) as many as 16 haplotypes from 17 Cyt b gene sequences. Meanwhile, in Thailand swamp buffalo there were 14 haplotypes from 62 Cyt b gene sequences (Yindee et al., 2010).
Based on the results of the haplotype analysis and the median-joining network of haplotypes ( Figure 5), two major haplotype groups were identified in the Indonesian swamp buffalo, namely haplotype 5 (H_5) and haplotype 2 (H_2). The major haplotypes act as radiating nodes for many other haplotypes (Kathiravan et al., 2011). H_5 is the haplotype as a representative of the Indonesian swamp buffalo which is widespread in all Indonesia regions, both in NAD, Java, and NTB or Sulawesi, especially North Toraja and partly Kolaka and Konawe. Meanwhile, H_2 is a haplotype as a representative of Southeast Sulawesi swamp buffalo. The haplotypes in the Bombana Island region and Bombana mainland region were specific breeds because they were only found in the second cluster and H_2 (see Figures 4  and 5). This condition was different from Konawe and Kolaka haplotypes which could be found in the first clusters and H_5 and also in the second clusters and H_2. Therefore, the haplotypes in the population of Bombana Island region and Bombana mainland region were feasible to be grouped as different breeds/strains from the other populations. The populations of swamp buffaloes in Bombana Island and Bombana mainland have the possibility of originating from an expansion of H_5 especially from the populations either from Kolaka or Konawe buffaloes migrating to Bombana Regency and experienced a mutation.
The other seven haplotypes (H-3, H_4, and H_6 to H_10) were expansion haplotypes from two major haplotypes. The formation of separate haplotypes from the major haplotype indicates the existence of groups or strains of Indonesian swamp buffalo that were genetically different. According to Coroian et al. (2015), the formation of a haplotype indicates that the population isolated from the viewpoint of geography and reproduction. Based on the polymorphic sites number ( Table 2) and percentage of haplotypes (Tabel 4), it can be claimed that the population of Southeast Sulawesi swamp buffaloes was relatively higher than the populations of buffaloes from Aceh, Banten, NTB, and North Toraja, where the highest polymorphic site haplotypes number were found in the populations of Konawe and Bombana Island.

CONCLUSION
The results of the study concluded that eight Indonesian local swamp-buffaloe populations were grouped into two main clusters. Bombana Island and Bombana mainland populations were specific breeds because they were only found in the second cluster and also have specific nucleotides sites on the 57 nucleotides which C base changed to T. The results of this study were useful in formulating the program of conservation and utilization of Indonesian buffalo genetic resources, especially in the buffalo population with specific breeds.

CONFLICT OF INTEREST
We certify that there is no conflict of interest with any financial, personal, or other relationships with other people or organizations related to the material discussed in the manuscript.