A total of eight populations of LRI were investigated in this study. Six of these populations were Tsaiya ducks. There are two major breeds of Tsaiya duck: White Tsaiya and Brown Tsaiya duck. After a period of selection and breeding, two populations of White Tsaiya and four populations of Brown Tsaiya were raised in LRI, including the GBT (Figure 1E), BI (Figure 1B), BII (Figure 1C), BIII (Figure 1D), Ilan White Tsaiya TLRI NO. 1, L101 (GWT) (Figure 1F) and Ilan White Tsaiya TLRI NO. 1, L102 (WI). Among the populations, the hybrid (H) was the progeny of GWT crossed with BI. The P (Figure 1G) was a breed which was raised by LRI. The population was also investigated in this study as a reference. The experimental protocols used in the present study were approved by the Experimental Animal Care and Use Committee of Ilan Branch, Livestock Research Institute, Council of Agriculture, Executive Yuan, Taiwan.

Blood samples were collected from 566 individuals from these 8 populations, including 43 GBT, 69 BI, 95 BII, 72 BIII, 80 GWT, 72 WI, 96 H, and 39 P individuals. For each bird, 3 mL of blood were drawn from the superficial plantar metatarsal vein or wing vein. Whole genomic DNA was then extracted with Genomic DNA Isolation Reagent (GenePure Technology CO., LTD, Taiwan) using the standard phenol-chloroform method.

Genomic DNA was isolated from blood taken from one male and one female Tsaiya duck using a previously described method; the following steps were slightly modified from the procedure previously described by Glenn and Schable [17]. The isolated genomic DNA was partially digested with RsaI and XmnI (NEB, USA) until most of the DNA fragments were between 300–1,000 bp in length. Re-naturing of single strand SuperSNX24 forward (5′GTTTAAGGCCTAGCTAGCAG AATC3′) and SuperSNX24+4P reverse (5′pGATTCTGCTA GCTAGGCCTTAAACAAAA3′) formed SuperSNX24 linkers that were ligated to the digested DNA fragments. Linked DNA fragments were then amplified with SuperSNX24. Next, microsatellite markers including fragments with biotinated probes were isolated and enriched. Three probe compositions were used: i) (TG)₁₂, (ACT)₁₂, (ACTG)₆, and (ACAG)₆; ii) (AG)₁₂, (ACAT)₈, (AACT)₈, and (AAGT)₈; and iii) (AAG)₈, (AAAC)₆, (AATC)₆, and (AGAT)₆. The biotinated probes then were annealed to fragments of gDNA containing complementary regions. Finally, the microsatellite-containing fragments were enriched using streptavidin-labeled metal beads (Dynabeads M-280 Streptavidin, catalog #11205D, Invitrogen, Carisbad, CA, USA). The enriched segments were then cloned into pGEM-T Easy vector (pGEM-T Easy Vector system, Promega, Madison, WI, USA) and sequenced with an ABI PRISM 3730XL DNA Analyzer (Applied Biosystems, Foster City, CA, USA). The fragments with higher repeats were then selected for polymorphism testing.

Seven published microsatellite markers were added to the marker set in this study, including APT004, APT005, APT008, APT015, APT020, APT025, and APT031[16].

The selected highly repeated fragments were subjected to polymerase chain reaction (PCR) and polymorphism testing which would verify if the microsatellite loci could be amplified and used to show diversity in the investigated populations. Primers for loci amplification were designed using Primer3plus [18]. CAG-tag (5′-CAGTCGGG CGTCATCA-3′) or M13Reverse (5′-GGAAACAGCTATGACCAT-3′) was added to the 5′ end of one of each primer pair [19]. Following the protocol described by Schuelke [20], a fluorescent dye-labeled tag, as a third primer, was used with the primer pair to amplify the target fragments which were detectable on the capillary electrophoresis. Thirty unrelated Tsaiya ducks from the LRI were tested. PCR was performed on a 20 μL volume using a thermalcycler (GeneAmp PCR system 9700, Applied Biosystems, USA) containing 0.5 U Taq DNA polymerase (TAKARA, Kyoto, Japan), 1×PCR buffer (1.5 mM MgCl₂), 0.2 mM dNTP, 0.2 μM unlabeled primer, 0.04 μM tag-labeled primer, 0.16 μM dye-labeled tag, and 50 ng gDNA. The PCR cycling program was as follows: 95°C for 5 min, 35 cycles of 95°C for 30 s, 50°C to 65°C for 40 s, 72°C for 40 s, and a final elongation at 72°C for 7 min. The amplified microsatellite PCR products were analyzed with a DNA analyzer (ABI PRISM 3730 DNA analyzer, Applied Biosystem, USA). Allelic sizes of all loci were estimated relative to in-line GeneScan500 LIZ Size Standard marker (ABI PRISM, Applied Biosystem, USA). The fragment size was calibrated and analyzed with Peak Scanner Software version 1.0 (ABI PRISM, Applied Biosystem, USA). The loci which had an allele number (N_a) greater than two and similar annealing temperatures were selected for whole population analysis.

For each locus and population and across populations, commonly derived statistics from the microsatellite genotypic data including allele frequencies, the observed number of alleles (No), the observed heterozygosity (H_O), the expected heterozygosity (H_E), and the polymorphic information content (PIC) were calculated with the Microsatellite Toolkit. The Hardy-Weinberg equilibrium (HWE) test was performed using the GENEPOP computer program, which also was used to estimate F-statistics (F_IT, F_IS, and F_ST) for each locus, the pairwise F_ST between populations, and the average inbreeding coefficient (F_IS). Nei’s genetic distance (D_A) between populations was measured with Microsatellite Analyzer. The phylogenetic tree was generated via the PHYLIP program using the neighbor-joining (NJ) method with bootstrap test of 1,000 resamplings of loci with replacement. The genetic distances of the proportion of shared alleles (POSA) was used to estimate and draw a POSA individual phylogenetic tree [21].

A hierarchical analysis of variance was carried out to allow the partitioning of total genetic variance into components owing to region, population, and individuals. Computations were carried out using a hierarchical analysis of molecular variance (AMOVA) procedure, as implemented in the ARLEQUIN 3.5 package. Principal component analysis (PCA) was performed with GENALEX v.6.501 software in order to spatially plot clusters and individuals based on the distance matrix with data standardization [21].

The model-based approach proposed for population structure analysis of the eight populations was carried out with the software STRUCTURE 2.3.1 [22] which assessed the genomic clustering (K) of the sample. To obtain a representative value of K for data modeling, ten independent runs were performed for each value from one to seven. The run length was set to 100,000 burn-ins following by 100,000 iterations. The ΔK estimated the most likely number of K that represented the population structure [23]. Subsequently, the CLUMPP software 1.1.2 [24] was used to demonstrate the optimal alignment of the 20 replicates for the same K value. The mean membership matrix across replicates was calculated with DISTRUCT v.1.1 [25].

Over 50 DNA fragments containing STRs were cloned. A polymorphism test was performed using the LRI populations, and 17 loci which had allele numbers more than two were picked. The 17 novel markers included three complex repeats (T1P4017, T3P1027, and T3P2090). In simple repeat markers, there was one dinucleotide unit (T1P1075), and four trinucleotide units (TYD037, TYD038, TYD042, and T1P2013), while others contained tetranucleotide units. These novel markers were then combined with the seven previously known markers [16], and together formed a set of markers for this study (Table 1).

Polymorphism was clearly observed at all the microsatellite loci in all eight of the LRI duck populations. The genetic characteristics of the 24 microsatellite loci are listed in Table 2. The average number of alleles per locus (N_a) was 11.292. The actual number of alleles ranged from 2 (T3P1027) to 36 (TYD029). The average number of effective alleles per locus (N_e) ranged from 1.816 (T3P2090) to 16.588 (TYD029), with an average across loci of 5.370. The PIC value ranged from 0.375 (T3P1027) to 0.937 (TYD029), with an overall average of 0.708. All of the selected microsatellite loci in this study were sufficiently polymorphic, indicating that these loci were suitable for the genetic analysis of ducks.

The HE among the 24 microsatellite loci had a range of 0.450 (T3P2090) to 0.941 (TYD029), with the average value of H_E being 0.746. The H_O among the 24 microsatellite loci had a range of 0.347 (T3P2080) to 0.808 (TYD024), with the average value of H_O being 0.591 (Table 2). However, there were four loci, namely, TYD021, T1P2013, T3P2090, and APT004, that significantly departed from the HWE (p<0.01).

The Wright’s F-statistic values (F_IS, F_IT, and F_ST) for each locus are shown in Table 2. The average F_IS for all the loci was 0.043, and the F_IS per locus varied from −0.093 (APT020) to 0.304 (TYD025). The average F_IT for all the loci was 0.233, and the F_IT per locus varied from 0.095 (T1P1075) to 0.397 (TYD006). The mean F_ST for all the loci was 0.197. This value implied that around 19.7% of the total genetic variation was caused by population differences and that 80.3% of total genetic variation was due to genetic differentiation among individuals within each population.

The genetic statistics relating to polymorphism, including H_E, H_O, PIC, the mean observed number of alleles, and the mean effective number of alleles (N_e), were calculated to estimate the allelic diversity at each locus of population. These genetic parameters across the 24 loci for the eight duck populations are listed in Table 3. H_E varied from 0.535 (WI) to 0.686 (BIII), whereas H_O varied from 0.482 (WI) to 0.707 (H) and PIC ranged from 0.461 (WI) to 0.633 (BIII). The WI population had the lowest values of H_O, H_E, and PIC. The BIII population exhibited the highest values of H_E, and PIC. The H population exhibited the highest values of H_O, and the smallest value of F_IS (Table 2).

Among the eight populations, the P population had the highest observed mean number of alleles (MNA) (6.542), followed by the hybrid (6.500) and BIII (5.833) populations, while the WI population had the smallest observed MNA (3.708). Negative F_IS values were observed only in one population (H), indicating an insufficient degree of inbreeding. The deviation from the Hardy-Weinberg proportions within populations (F_IS) varied from −0.047 to 0.103. The highest inbreeding effects were found in the BI population (0.103). The H population contained 16 loci that were significantly deviated from the HWE (p<0.01) (Table 3), and were much higher than other populations.

The values of F_ST and D for each test population pair are summarized in Table 4. The F_ST for each population pair was highly significant (p<0.05) The F_ST values of the population pairs were varied from 0.083 (for the BI and H population pair) to 0.322 (for the WI and BII population pair). There were several pairs of population have sub-lowest F_ST of which around 0.09 (BI and GBT pair, BI and BIII pair, H and BIII pair and H and GWT pair).

The genetic distances between the duck population pairs varied from 0.157 (for the BI and H population pair) to 0.582 (for the BII and WI population pair). There were relatively high genetic distance among the white duck populations (GWT, WI, and P) and between the other populations (p>0.3) except for GWT and H pair while GWT was paternal of H.

In this study, the results of the differentiation of the eight duck populations using the NJ tree are shown in Figure 2. In the NJ tree, the brown duck populations (GBT, BI, BII, and BIII) and the white duck populations (GWT, WI, and P) were independently grouped into two clads. The closest relationship in brown duck group was GBT and BIII. The following population was BI, but their bootstrap confidence value was only 56%. The furthermost population among brown ducks were BII. In the white duck clad, P located apart from GWT and WI, of which bootstrap value was 68%. The H populations were depicted as independent taxa.

A PCA of pair-wise genetic distances among the eight examined goose populations was used to represent the relative positions of the populations. The first (PC1), second (PC2), and third (PC3) principal components accounted for 25.95%, 23.62%, and 20.60% of the total variation, respectively (Figure 3). The results of the PCA were similar to the phylogenetic tree drawn up via the NJ method except that GWT and WI was separated in PC3.

The results of the degree of variance in the duck populations via AMOVA are summarized in Table 5. In the analysis of all the populations, the largest variation was found within individuals (76%), followed by variation among populations (19%). The variation among individuals within populations accounted for 5% of the total variation. In the Brown Tsaiya duck populations, the variation among individuals within the populations accounted for 13% of the total variation, while that among populations accounted for 7% of the total variation. In the White Tsaiya duck populations, the variation within individuals decreased to 70% while that among populations reached to 24%, while the variation among populations in Brown Tsaiya duck populations was only 13%.

The results of the STRUCTRE program analysis of all the populations are shown in Figure 4A. White Tsaiya and Brown Tsaiya were definitively separated into different clusters just at K = 2, while H and P became mixed groups. At K = 3, the GWT population was split into a new cluster. At K = 4, the BII population formed another new cluster. All populations were assigned into their cluster except that GBT bound to BII and H just has half BI and half GWT while K was set at six. Seven population excluding H were almost have their clusters until K = 8. No new cluster appeared, however, and the figures were not much different for values of K larger than 8. In order to survey the structures on intra-clad basis accurately, analyses of the white ducks and brown ducks, respectively, were carried out (Figure 4B, 4C). Among the brown ducks and P, the BII population became a subgroups first at K = 2, while the P population subsequently became a new subpopulation at K = 3. Four separately clusters were form at K = 4, that only GBT was a mixed population with the components of BI and BIII. At K = 5, the brown ducks and P were split into 5 clusters exactly. BII population tended to have two subgroups at K = 6. The white ducks were relatively simple. At K = 2, GWT formed a completed cluster and WI and P were the similar mixed clusters. The white ducks were split into 3 clusters exactly at K = 3. When K was over 4, the components of P had more complicated than GWT and WI.

In Table 3, the results present the intra-population genetic variation in eight of the Tsaiya duck populations. Deficiencies of heterozygosity (that is, H_O lower than H_E) were exhibited by almost all the populations. Only the H population had excess heterozygosity, while the GBT and GWT populations each had roughly equal H_O and H_E values. The sources of heterozygous deficit probably included inbreeding, null allele, aneuploidy, misscoring and Wahlund effect [32]. There were all small and closed populations in LRI, however, all populations but GWT and GBT were established by selection of specific traits. GWT and GBT were derived from the original populations without any major selection since before 1980 [7]. The major causes of deficiencies of heterozygosity perhaps were selection of specific traits that accompanied specific alleles. The BII population tended to split into two Brown Tsaiya subpopulations when K = 6 (Figure 4B). The deficiencies of heterozygosity on BII maybe combined a Wahlund effect. The first generation of BI crossed GWT have some characters predictably: the excess heterozygosity (H_E = 0.676, H_O = 0.707), far from inbreeding (F_IS = −0.047) and departure HWE (16 of 24 markers HWE departure). In general, except for H, every Tsaiya duck population in LRI was a relatively low F_IS (almost lower than 0.1), and the numbers of loci departing from the HWE were less than 7. It can thus be concluded that the Tsaiya ducks in LRI probably have good breeding programs. The results of the intra-population diversity analyses for the Pekin populations from the LRI were similar, including relatively high variation (H_O = 0.561, PIC = 0.589 and the highest observed MNA = 6.542) and F_IS values of 0.08. WI was selected for crossing Muscovy drakes to improve the carcass appearances of mule ducks. WI had related low diversity (PIC = 0.461; H_O = 0.482) and little more inbreeding (F_IS = 0.098). This was probably caused by the multiple and long-term selection.

Brown Tsaiya and White Tsaiya ducks are the two major domestic breeds in this study, in addition of P and a hybrid population. According to the AMOVA results (Table 5), the proportion of genetic variation attributed to population differences in all populations was about 19%, while the proportion of genetic variation attributed to individuals within populations was 5% of the total genetic variation. Compared to other researches, differentiation was analyzed with 15 microsatellite loci in 26 breeds of Chinese indigenous ducks, and 12.43% of total genetic variation was attributed to differences among populations [33]; four breeds of Indian indigenous ducks were analyzed with 24 microsatellite loci, and 19.76% of the total genetic variation was attributed to differences among populations [34]. However, the Tsaiya populations of LRI were derived from unique species: Brown Tsaiya ducks. The formation of these breeds was just a series of selections and different periods of breeding. The differentiation among populations of Tsaiya ducks in LRI was larger than many different breeds in China and India. The only inference was that higher variance of markers were used in this study [35]. In addition. The variation attributed to individuals within populations was 5%. The high inter-population and low intra-population variations was very benefit to discriminate the different populations that even had similar appearances. The individual phylogenetic tree (Figure 5) revealed that not only separated entirely eight populations, but the individuals were assigned exactly to their own population. Although all Brown Tsaiya populations were very close in PCA (Figure 3). Such high resolution, the 24 microsatellite loci became an excellent tool to monitor the genetic structure of each Tsaiya ducks population, and to consult a breeding program for breed conservation. When AMOVA analysis focused on Tsaiya ducks (six populations), the variance among populations (22%) were comparatively more heterogeneous than all populations (19%), and phenomenon which might have been caused by the differences between the Brown and White Tsaiya ducks. This hypothesis was confirmed with the drawing of the NJ phylogenetic tree, which showed two major clads of brown and white (including P) ducks (Figure 2), and through the STRUCTURE analysis, which revealed that the brown and white populations were divided into two parts when K = 2, with the H population split into two contributions (Figure 4A). The PCA plot showed that the brown and white ducks were significantly separated in terms of PC1, while the hybrid ducks were located between them (Figure 3). The individual phylogenetic tree (Figure 5) also divided the brown and white ducks into two clusters, with the hybrid ducks and P inserted between them. In our distinct analyses, the H population was always situated between two major groups or bisected individuals at K = 2 to K = 6, especially two parts just composed GWT and BI while almost populations were assigned to specific clusters. The 24 microsatellite markers used in this study have the ability to detect any crosses between two breeds.

The Brown Tsaiya populations in LRI were originated from GBT population and selected and derived to BI population, then BII and BIII populations were derived from BI. Moreover, LRI have been bred these four populations and have been kept their characters for becoming a breed. During such period, Monitoring their genetic characteristics and structures was very important. In AMOVA, the proportion of genetic variation of brown duck populations that came from population differences was only 13%, less than all Tsaiya duck populations, 22%. The four populations were close and compact as a group that far from other populations in 3D PCA. It’s perhaps all BI, BII, and BIII were derived from GBT and all were selected by only one trait, moreover their separation was not for long time, the selective programs began at 1984, 1992, and 1996 respectively [6,8,9]. Such low diversity between populations and each population would breed for a breed, the monitoring and analyzing system should be high definition. The unrooted individual phylogenetic tree (Figure 5) showed that four Brown Tsaiya duck populations were close but parted. The STRUCTURE analysis of the Brown Tsaiya ducks (Figure 4B) also revealed that every individual that belonged to the same line was simply assigned to the same population, even when the portions in an individual were almost close to 100%. Our microsatellite markers had high definition and were suitable for analyzing the genetically close populations. Unexpectedly, the Brown Tsaiya and P populations were not isolated until K = 3. BII was isolated first at K = 2. However, P located on different clad of Brown Tsaiya ducks in NJ phylogenetic tree and separated from Brown Tsaiya duck cluster at PC2 in PCA. In our analysis, P still was a very different breed from Tsaiya duck. As the expectation, BII was split to a specific group first at K = 2, matching the results of neighbor joining phylogenetic tree and PCA that BII was a peripheral group of Brown Tsaiya ducks. It was maybe that BII raised and bred on Tainan, south of Taiwan, nevertheless, other Brown Tsaiya duck populations were on Ilan branch, north of Taiwan. The geographic separation and continuous selection [9] probably were the causes of this phenomenon. As BII, BIII was selected from BI, but its variation or distance to GBT was even less than BI’s. It could be caused by that genetic structure of BIII, picked with blue shell eggs, was closer to GBT. Otherwise, BIII and GBT have the highest H_O (BIII = 0.631, GBT = 0.661) and PIC (BIII = 0.633, GBT = 0.623) (Table 3). High polymorphism perhaps make both of them similar genetic content. The BII population was even split into two clusters at K = 6 (Figure 4B). Selection in a closed population may have been the cause of this.

In contrast, belonging a clad in phylogenetic tree, the white duck populations (including P) were dispersal in PCA than brown ducks, that WI was split away from GWT at PC3 (Figure 3). The separating far away their co-original population maybe was caused by the progeny test for the white hybrid mule that needed long term selection. The diversity of WI was lost during the breeding procedure that was expressed in the lowest H_O (0.482) and observed MNA (3.708) (Table 3). Chang et al [10] reported the same trend that there were just three clearly distinguishable populations at K = 3, and that even at K = 4, two individual White Tsaiya ducks were still assigned fully to their own population (Figure 4C). It showed that two white ducks populations were totally differentiation under analysis of our microsatellite markers. On the other hand, P had split into two sub-populations while K = 4.