An assessment of genetic variation in vulnerable Borneo Ironwood Eusideroxylon zwageri Teijsm. & Binn. in Sarawak using SSR markers

: Borneo Ironwood Eusideroxylon zwageri Teijsm. & Binn. has high market value for its valuable and durable timber, which has put it at risk due to illegal logging. This study analysed E . zwageri genetic variation using four microsatellite markers in populations at Nirwana Rehabilitation Forest (NRF), and Tatau, Sarawak. We found that 20.1% of total genetic variation corresponded to differences between populations, while 79.9% was attributed to differences among individuals from the same population. The Tatau population had lower genetic diversity compared to NRF, and both populations showed depressed heterozygosity indicative of inbreeding. Allelic data were also used to confirm variety level differences proposed by earlier workers, and three informal varieties: zwageri, grandis , and exilis were recognized in the study area. It is expected that the results from this study could serve as baseline data for conservation of this vulnerable species.


INTRODUCTION
The Borneo Ironwood Eusideroxylon zwageri Teijsm.& Binn. is one of the most treasured and crucial commercial timber trees endemic to the Asian forest (Malaysia, Indonesia, Brunei, and The Philippines).It is known locally as 'Belian' in Malaysia, 'Ulin' in Indonesia, and 'Tambulan' in The Philippines.The species belongs to the Lauraceae family that includes the avocado, bay laurel and cinnamon tree.It has been listed as 'Vulnerable' on the IUCN Red List of Threatened Species due to over-exploitation and habitat destruction (Asian Regional Workshop 1998).This species will remain endangered unless circumstances threatening its survival and reproduction improve.Slow growth rates (mean radial growth rate is 0.058cm per year; Kurokawa et al. 2003) and slow regeneration in logging areas also contribute to lower numbers of forest trees.
Many taxonomists have reported variety level morphological differences in E. zwageri (Teijsmann 1858;Teijsmann & Binnendijk 1863;Koopman & Verhoef 1938;Kostermans et al. 1994;), however, no valid taxonomic treatment for these varieties has been proposed so far.Some of the vernacular names given to the varieties were 'Belian telor', 'Belian kapur', 'Belian sirap', 'Belian tanduk', 'Belian tembaga', and 'Belian lilin'.Nevertheless, Irawan et al. (2016) tried to confirm the presence of these varieties using amplified fragment length polymorphism (AFLPs) studies.The study showed promising results, with 98% out of a total of 50 samples clustered according to the varieties recognized by the local people in Indonesia.The four varieties were informally recognized as var.zwageri, var.exilis, var.grandis, and var.ovoidus (Irawan 2005a,b).Local people in Sarawak also recognized these varieties based on the differences in fruit's form, bark or wood structure (Marzuki pers. comm. 18.ix.2017),however, no valid taxonomic treatment has been published.
This tree can reach a height of up to 50m and may live over 1,000 years (Global Trees Campaign 2020).Mature trees produce large fruits that, although poisonous to humans, are important food sources for foraging animals.The species is also valued for cultural reasons.The wood is dense (0.85-1.1 g/cm 3 ) (Irawan 2016), strong and resistant to decay, making it preferred by indigenous people of Borneo for building houses.The Dayak people of Borneo believe the tree protects them from dangerous animals, while 'Murut' (Borneo headhunter) use it to make blowpipes and Dusun ancestors used it to create coffins.The black pepper industry in Borneo has traditionally used Belian wood as a support to the creeping herbs.In the famous Murut Cultural Center of Malaysia, Belian wood pillars are used.
The revised status of E. zwageri in Sarawak under Criteria B of IUCN Red List (Md-Isa et al. 2021) indicates the need to formulate conservation plans to protect this species from extinction.In addition to economic or social information, occurrences and distribution patterns of species, genetic information of the species is another important aspect that needs to be considered in conservation action plans.Habitat fragmentation can contribute to the reduction of genetic diversity of this species.Although transplantation of E. zwageri from other locations is a practical conservation strategy, an accurate understanding of the genetic structure of natural population of E. zwageri is necessary for conserving them.This is because relocation of the species or reduction in size for other reasons will cause the loss of genetic diversity in the new population through genetic drift (Lowe et al. 2005;Finlay et al. 2017).
Little is known about the genetics of E. zwageri, especially in Sarawak.Two studies in Indonesia using randomized amplified polymorphic DNA (RAPD) markers revealed 96% genetic diversity of E. zwageri within populations, and the remaining attributed to population differences (Harkingto et al. 2006;Rimbawanto et al. 2006).This was probably due to the samples that originated from the same population.Nurtjahjaningsih et al. (2017a,b) showed high genetic diversity for E. zwageri in Indonesia, and suggested transplantation among different populations should be conducted with careful consideration.A few studies on genotyping of E. zwageri using direct amplified minisatellite DNA (DAMD) marker (Yoon 2006), M13 universal marker (Siew 2005) and RAPD marker (Hong 2005) were mainly aimed at identifying and genotyping the two genera (Eusideroxylon and Potoxylon) known by similar common names Belian in Sarawak.
In this study, two different habitats were chosen to study the genetic variation of E. zwageri in Sarawak.One was a restoration forest, Nirwana Rehabilitation Forest (NRF) located in Universiti Putra Malaysia (UPM) Bintulu Campus.The other was a fragmented forest in Tatau, Bintulu, Sarawak.This study was conducted to assess the genetic variation of E. zwageri by using four highly polymorphic microsatellite markers recently developed for the species (Kurokochi et al. 2014).We compared genetic variation between the two populations and within each population, and determined the level and pattern of genetic variation in both areas.Allelic data were also analyzed for the presence of variety level J TT differences in collected samples.

Sample collection
Two sampling sites were selected as model habitats for this study; 1) Nirwana Rehabilitation Forest (NRF) in UPM Bintulu Campus, and 2) fragmented forest area in Tatau, Bintulu, Sarawak (Figure 1).Samplings were conducted in April 2016, August 2016, and September 2017, over the periods of two weeks each.Leaves were collected from 52 trees, of which 39 from NRF and 13 from Tatau forest.A single leaf was collected per tree.The leaf materials were kept in silica gel prior to DNA extraction.

DNA extraction
Total genomic DNA was isolated using conventional Cetyl-Trimethyl Ammonium Bromide (CTAB) method (Doyle & Doyle 1987) with some modifications.The leaf materials were ground with CTAB buffer until fine paste and transferred into sterilized 1.5mL microcentrifuge tube.In the fume hood, 500μL of 2X preheated (at 65 o C) CTAB extraction buffer [2% (w/v) hexadecyltrimethylammonium bromide (CTAB); 1.4 M sodium chloride (NaCl); 100 mM Tris-HCl, pH 8.0; 20 mM ethylenediamine tetra-acetic acid (EDTA), pH 8.0; 1-2% (w/v) polyvinyl-pyrrolidone (PVP-40T)] and 2μL of 1% β-mercaptoethanol were added to each sample in the 1.5mL micro-centrifuge tube and vortex gently until mixed well.The homogenized mixture was then incubated for 20-30 minutes at 50 o C in a water bath and inverted every five minutes (Md-Isa 2020).The tubes were then transferred to another water bath and incubated at 65 o C for another 15 minutes.The samples were allowed to cool slightly before adding 500μL of chloroform: isoamyl alcohol (24 : 1) and inverted to mix.At room temperature, the samples were gently shaken for 15 minutes and centrifuged for 10 minutes at 12,500rpm.About 400μL of supernatant was transferred into a new sterile 2.0mL screw cap micro-centrifuge tube (Md-Isa 2020).
The DNA was precipitated by adding 800-1,000 μL of 95% cold ethanol and inverted gently, and allowed to precipitate up to three hours or longer in -20 o C freezer.The tubes were then centrifuged for 10 minutes at 12,500rpm to pellet the DNA.The supernatant was J TT discarded and the pellet was washed with 500μL of 80% cold ethanol and gently mixed for 10 minutes.Then, the pellet was centrifuged for five minutes at 12,500rpm and the supernatant was discarded again.The DNA pellet was allowed to air dry at room temperature before resuspending with 100μL of TE buffer or distilled water.The DNA samples were then kept at -20 o C for further usage.
Afterward, the DNA quantity was estimated through electrophoresis on 1.0% (w/v) agarose gel.The gel was run in tris-boric acid-EDTA (TBE) buffer with EtBr "Out" Staining Solution (YEASTERN BIOTECH Co. Ltd) at 90V for 30-45 minutes and quantified in comparison to the 1kb DNA ladder (Promega, Madison, WI, USA) with known concentrations.The band was visualized using UV transilluminator.The band intensity of the DNA product was quantified to the intensity in the ladder (Promega 2020).The images were captured with DOC PRINT system (Vilber Lourmat, USA).

Polymerase chain reaction (PCR)
Polymerase chain reactions (PCR) were performed in a volume of 25µL using 2µL of DNA template (5-10 ng), 10µM of each primer, 5X GoTaq buffer, 10mM dNTPs, and 25mM MgCl using an Eppendort AG 22331 Mastercycler.PCR cycling conditions were as follows: a single cycle of pre-denaturation for two min a 95 o C, followed by 40 cycles, each consisting of 30 sec denaturing at 95 o C, 30 sec at annealing temperature 50 o C, and 30 sec elongation at 72 o C, the last cycle ending with a single cycle of final extension at 72 o C for five min.
The PCR amplification of the DNA sample was carried out using four highly polymorphic microsatellites developed for E. zwageri (Kurokochi et al. 2014).Forward primer of each marker was labeled with fluorescent dye (Table 1).One representative of PCR product from each primer was subjected to one direction sequencing done by First Base Laboratories Sdn.Bhd.(Seri Kembangan, Selangor, Malaysia) to confirm the amplification of the microsatellite repeat region.The products were then subjected to bi-directional sequencing for fragment analysis (FA).FA were carried out to detect changes in the length of a specific DNA sequence to indicate the presence or absence of the microsatellite marker through detection of fluorescent label in the PCR product.The size standard (500-ROX) was combined with the sample of interest and co-injected on the capillary electrophoresis system (ABI3730XL Applied Biosystems Genetic Analyzer).

Data analysis
For data analysis, a total of 52 leaves samples were distributed into two populations with Tatau (13 samples) as one population and NRF (39 samples) as one population.The scoring of allele sizes were performed using GeneMapper® version 4.0 analysis software using the service provided by First Base Laboratories Sdn.Bhd.The alleles nearest to the expected PCR product size were recorded; while the non-specific products, which were out of range, were ignored.Allele frequencies per locus and per population were analyzed by FSTAT version 2.9.3 (Goudet 1995).Number of alleles (Na), expected heterozygosity (He), observed heterozygosity (Ho) and polymorphism information content (PIC) were also estimated using Cervus version 3.0.F-statistics, including inbreeding coefficient of an individual relative to the subpopulations (F IS ), inbreeding coefficient of an individual relative to the total population (F IT ), and genetic differentiation index between population (F ST ) were calculated using GenAlex version 6.501 (Peakall & Smouse 2012).The software was also employed to determine genetic diversity within each population (Na, Ho, He, and F IS ).
Further, the allelic data of the 52 samples were subjected to estimation of genetic distance among genotypes using simple matching coefficients by bootstrapping 1,000 times and then were clustered using unweighted neighbor-joining method by using

J TT
Dissimilarity Analysis and Representation for Windows (DARwin) version 6.0.21 (Perrier & Jacquemoud-Collet 2006).This analysis was done to see if there are any variety level differences in the collected samples.The names of varieties used in the present study were adopted based on earlier studies (Irawan 2005a,b;Irawan et al. 2016).

RESULTS
The 52 samples from E. zwageri trees were scored for all four microsatellite DNA loci.The amplification of the microsatellite repeats was confirmed by sequencing, where specific repeat motifs were successfully identified.The PCR reactions which failed to produce sufficient product for genotyping were recorded as missing data for all analyses, however, locus Ez-04 yielded less than 20% amplification and was discarded for the subsequent analysis.
Genetic variations among three markers tested for all 52 individuals are summarized in Table 2.In total, 25 alleles were detected at these three loci in 52 individuals, with the number of alleles per locus ranging from 3 (Ez-09) to 12 (Ez-05), with an average of 8.33 alleles per locus.Ez-09 showed lowest number of alleles being amplified but it was detected in 50 individuals with 38 individuals from NRF and 12 individuals from Tatau compared to Ez-05 detected in 42 individuals with 35 individuals from NRF and seven individuals from Tatau (Appendix 1).
Meanwhile, observed and expected heterozygosity values of all three loci ranged 0.511-0.720and 0.664-0.867,respectively.Whereas, the average observed heterozygosity for both populations (Ho= 0.593) was lower than the average expected heterozygosity (He= 0.791), which may indicate moderate levels of genetic variation in the Belian populations studied.While, the polymorphic information content (PIC) value for all three loci were higher than 0.5 ranging from 0.583 to 0.841, which suggested all three loci used in this study were highly polymorphic.
Besides, locus Ez-05 showed the highest value of allelic richness (A R = 7.117) among the three markers with total of 12 alleles being amplified from both populations; however, among the 12 alleles amplified, only six alleles were identified in Tatau compared to 10 alleles in NRF.The distributions and allele frequencies of the three loci in both populations are shown in Figure 2A-2C and listed in Appendix 1.
Furthermore, F-statistics were estimated in a fixation index as shown in Table 2.The average inbreeding coefficient of the individuals to the total population (F IT ) was 0.295 while the average inbreeding coefficient of the individuals to the subpopulation (F IS ) was 0.048.Whereas, the average genetic differentiation (F ST ) of the subpopulation compared to the total populations was 0.201.The average value of F ST indicated that about 20.1% of total genetic variation corresponded to differences between populations, while 79.9% was explained by differences between individuals of the same populations.And the lower F ST value, which is less than 0.25, may suggest that there is gene flow between the two populations.
The genetic indicators within each population are summarized in Table 3.The data showed higher number of alleles (Na= 7.333) being amplified in 36 individuals in NRF populations.The observed heterozygosity value (Ho= 0.659), however, was lower than the expected heterozygosity (He= 0.739), which indicated moderate level of genetic variation in NRF populations.Comparatively, in Tatau, lower number of alleles (Na= 4.000) was found in nine individuals.Moderate level of genetic variation was also observed in Tatau population by the lower number of observed heterozygosity (Ho= 3.99) than the expected heterozygosity (He= 0.563).
Generally, the lower number of observed heterozygosity than the expected heterozygosity was J TT also evidence that both populations deviated from Hardy-Weinberg equilibrium.Both populations also showed a deficiency of heterozygosity, indicated by positive F IS values (NRF= 0.054; Tatau= 0.165).Additionally, the allelic data for 52 samples were analyzed based on the estimation of genetic distance among genotypes to segregate the individuals according to their varieties.The unweighted neighbor-joining dendogram grouped the 52 samples of the two populations into three varieties (Figure 3).Of the 52 samples, 15 samples were grouped together as variety exilis, 16 samples as variety grandis and 21 samples as variety zwageri.

DISCUSSION
The informativeness of observed loci across the two populations was measured based on polymorphic information content (PIC).Theoretically, PIC values range 0-1 (Hilderbrand et al. 1994).At a PIC of 0, the marker has only one allele, and at a PIC of 1 the marker would have an infinite number of alleles.Thus values of PIC greater than 0.5 are considered to be highly informative.Data from the current study resulted in an average value of PIC= 0.746.Therefore, all the three loci used in this study can be classified as highly informative loci (PIC >0.5) and appropriate for assessing genetic variation.
Based on the average value of expected heterozygosity (He= 0.791), a moderate level of genetic variation among the 52 individuals studied was obtained.Nevertheless, when compared the genetic variation between NRF and Tatau, NRF showed higher genetic variation within population compared to Tatau, where theoretically, the wild population should have higher genetic variation (Pandey et al. 2004;Gauli et al. 2009).This may be due to the source of NRF trees, which originated from several places.It may not only limit to Bintulu area but from several places, which contribute to the high level of genetic variation in the populations, However, inbreeding depression may still occur because of the small size population.As for wild population of Tatau, the genetic variation was lower than expected, probably due to the small population size in an island forest fragment within palm oil plantation, which may contribute to the low genetic variation in the population.
In addition, heterozygote deficiency was detected, which was depicted by the lower average value of observed compared to expected heterozygosity (Ho < He).It suggested that both populations might be inbred.This was also evidenced by the positive average value of F IS (Table 2), which observed a stronger inbreeding in the Tatau population than the NRF population.Heterozygote deficit can be explained by various factors such as    J TT non-random mating, unamplified alleles (null alleles) and inappropriate sampling (population admixture Wahlund's effects) (Borsa et al. 1991;Castric et al. 2002;Dharmarajan et al. 2012;Waples 2015).
Inbreeding often results from a population bottleneck (genetic drift) due to anthropogenic or environmental events.In this study, higher inbreeding depression was observed in Tatau.The small population size and limited number of standing trees in the area might be the main cause of the inbreeding depression.Other factors could be founder events, where a population has reduced genetic variation compared to the original population and thus produces an apparent high level of inbreeding.This phenomenon may happen in the NRF population where seeds from several unknown populations in Bintulu and outside Bintulu added the small size of the NRF population.The different source of seeds, however, may contribute to the high genetic variation in NRF compared to Tatau.
Furthermore, based on Wright (1978), F ST = 0-0.05indicates little population differentiation, F ST = 0.05-0.15indicates moderate differentiation, F ST = 0.15-0.25 indicates high differentiation, and F ST >0.25 indicates highest differentiation.Current study revealed high genetic differentiation between the two populations (F ST = 0.201).The genetic differences might be due to the J TT extensive geographic range (population isolation) of the species and the small population size.
This fixation index (F ST ) may also provide approach for estimating inter-population gene flow (Nm) (Avise 2004).Average gene flow (Nm= 1.826, Table 2) in the current study may indicate one incoming migrant per generation in each population when F ST = 0.201 (Wright 1931).This indicates some migration (gene flow) at low rate between the two populations and low levels of interbreeding.
Additionally, both populations showed a possible deviation from Hardy-Weinberg equilibrium (HWE), which can be observed through the lower number of observed heterozygosity (Ho) than expected heterozygosity (He).Another possible parameter to look for the deviation from HWE in this study is the F-statistics.Populations are in Hardy-Weinberg equilibrium if F IS = 0 and F IT = F ST (Guries & Ledig 1981) The deviation from HWE can be caused by several factors such as mutation, migration, random mating, selection and small size population (Keats & Sherman 2013;Johnston et al. 2019).In this study, the disequilibrium was probably caused by the small population size of both study sites.Consequently, sampling error is unavoidable in this study.Small sample size contributes to the violations of the HWE principles.Other factors such as migration or seed dispersal by human or animal were observed in NRF, which also contribute to the deviation from HWE.While in Tatau, dispersal of seed was observed mainly through the river where the population is situated.Inbreeding and population isolation as discussed above may also contribute to the deviation from HWE.
Overall, results from this study should be used with caution, as the number of marker and samples tested in this particular study are not sufficient to make any conclusive statement.The PIC of all the three markers, however, is very high, in this case the average is 0.746, which is more than 0.5 and suitable to use for the analysis.Nonetheless, future study where more samples and more markers can be included should be carried out to have more comprehensive understanding on the genetic variation of E. zwageri especially in Sarawak.
The allelic data obtained in this study were also used to segregate the samples according to their varieties based on estimation of genetic distance among genotypes.The result showed that the 52 samples were grouped into three varieties.Morphological study conducted in the same sampling area documented similar variations (Md-Isa 2020).The present result demonstrates the existence of variety level differences in E. zwageri.Nevertheless, further study on using more markers, samples and DNA barcoding will verify and give more concrete answer for this finding.Thorough work on taxonomy classification is ongoing to validate the taxonomic status of the varieties in E. zwageri.

CONCLUSION
Genetic analysis of the two populations of E. zwageri in Bintulu, Sarawak shows that 20.1% of total genetic variation corresponded to differences between populations.The Tatau population was observed to have relatively lower genetic diversity compared to NRF area.Therefore, it was suggested that establishment of the restored population (NRF) from a limited number of individuals originated from unknown places has higher level of genetic variation compared to the wild population (Tatau).This occurrence, however, resulted in the inbreeding due to genetic drift (population bottleneck and founder effect).Furthermore, inbreeding can lead to fixation of deleterious alleles and may lead to a decrease in the genetic variation of the population.
Consequently, it is important to obtain more detailed information on the genetic variation of this species in order to form an effective conservation strategy.We suggest further conservation efforts focused on ensuring suitable habitat for the continued recovery of this species.Effort can be made to identify locations of E. zwageri as protected forest areas.This will facilitate natural regeneration without disturbance.Alternatively, sprouting and cutting techniques can be proposed, as the sprouts tend to grow faster and may reach mature stage faster than regenerating them from seedling (Mostacedo et al. 2009).
In addition, the findings on the segregation of the species into three varieties based on the allelic data are promising.Validation of the names of the varieties can now be proposed with effective taxonomic publication.It would also be interesting to establish DNA barcodes for this species to confirm our observations.More samples from different locations and different genetic markers are suggested to be included to support future on this topic.

Figure 1 .
Figure 1.The locations of the sampling sites in Sarawak, Malaysia.

Figure 3 .
Figure 3. Unweighted Neighbor-joining trees using simple matching similarity coefficient based on three microsatellite markers for the 52 individuals of Eusideroxylon zwageri from Nirwana Rehabilitation Forest and Tatau, Bintulu.The tree shows the clustering pattern of three varieties of E. zwageri namely; variety exilis, variety grandis and variety zwageri.The scale bar (0-0.2) represents the level of dissimilarity.
. The value of F IS ranges between -1 and +1.Negative F IS values indicate heterozygote excess (outbreeding) and positive values indicate heterozygote deficiency (inbreeding) compared with HWE expectations.Result from this study showed a positive F IS value, (NRF= 0.054; Tatau= 0.165) which indicates deficiency of heterozygosity, hence, a possible deviation from HWE.
unweighted (UW) allele frequencies per locus and per population.N is allele size (base pair).

Table 1 . Characteristics of four polymorphic microsatellite loci in Eusideroxylon zwageri (Kurokochi et al. 2014) used in the current study.
Na-number of allele | Ho-observed heterozygosity | He-expected heterozygosity | Tm-melting temperature.