Genetic diversity and population structure of muscovy duck (Cairina moschata) from Nigeria

Ethical considerations

All experimental procedures in the present study were performed in accordance with the Research Guidelines for the Institutional Review Board of Kunming Institute of Zoology, Chinese Academy of Sciences (SMKY-20160105-11) and the University Ethical Review Committee, University of Ilorin, Nigeria (UERC/ASN/2021/2161).

Sampling and data collection

We sampled a total of 287 domestic Muscovy ducks (155 males and 132 females) from farmer’s flocks in six different States in Nigeria as follows (Fig. 1A; Table S1); Sokoto State (n = 7 males; n = 9 females), Kano State (n = 36 males; n = 26 females), Taraba State (n = 4 males; n = 5 females), Oyo State (n = 44 males; n = 29 females), Ekiti State (n = 26 males; n = 34 females), and Bayelsa State (n = 38 males; n = 29 females). Our sample collection considered genetically unrelated domestic Muscovy ducks with geographical and ecological perspectives from the six states in Nigeria. This included random sampling of at least two matured birds per flock from domestic Muscovy ducks raised on free-range scavenging system with a distance of at least 500 m apart. Information on the pedigree of the birds were obtained from their owners before blood collection. Blood samples were collected from each bird through their wing venipuncture with the assistance of trained veterinary personnel. Blood samples were preserved in 95% ethanol at room temperature thereafter transported to the laboratory and stored at 4 °C for instant use, or −80 °C for extended storage.

DNA extraction, PCR and sequencing

Genomic DNA was extracted from whole blood using the standard phenol/chloroform method (Sambrook & Russell, 2001). We examined genomic DNA for both quality and quantity using 2% Agarose Gel Electrophoresis against a 2 kb DNA ladder marker and Thermo Scientific™ NanoDrop 2000 spectrophotometer respectively.

Two overlapping primers were used to amplify 940 bp of cytochrome b: forward for the first set DUKCYB1L (5′-ATCTTTCGCCCTATCCATCC-3′) and reverse DUKCYB2R (5′-TTTGGTTTACAAGACCAATGTTTT-3′), and the second, P105 (5′-GCCTCCTGCTAGCCATACAC-3′) and P106 (5′-TACGGCGGGAAAGAGAAATA-3′) (Kameshpandian, Thomas & Nagarajan, 2018). The cytochrome b amplification reaction was performed in a 50 µl PCR mixture containing 5 µl 10× reaction buffer, 1.5 mM MgCl₂, 0.2mM dNTPs, 0.2 µM each primer, 1.5 U Takara Taq DNA polymerase, and approximately 30 ng genomic DNA. The amplification conditions: 95 °C for 5 min, 30 cycles at 94 °C for 30 s, 54 °C for 45 s, 1 min for 72 °C and a final extension of 72 °C for 5 min (Kameshpandian, Thomas & Nagarajan, 2018). After, the quality of PCR products was confirmed on 2% agarose gel.

In addition, 747 bp fragment of nuclear DNA CYP2U1 was amplified using forward primer P103 (5′-GTTATTTGGTTATGCATATCGTG-3′) and reverse primer P104 (5′-GAGACGGTTGGCGTATATGG-3′) (Kameshpandian, Thomas & Nagarajan, 2018). Similar PCR conditions as above was used except for annealing temperature of 45 °C. The same quality check on 2% agarose gel was done for nuclear DNA CYP2U1 PCR products. The amplified cytochrome b and CYP2U1 fragments were purified with Exo-SAP-IT cleanup kit as per manufacturer’s instructions (Affymetrix). Sequencing reactions were performed in both directions using the BigDye™ Terminator Cycle Sequence Kit 3.1 Ready Reaction Cycle Sequencing Kit (ABI Applied Biosystems, Waltham, MA, USA) and the products afterwards were purified through alcohol precipitation. Finally, ABI PRISM 3,730 automated DNA sequencer (ABI Applied Biosystems, Waltham, MA, USA) was used to analyze the purified products. We assembled all sequences in SEQMAN PRO of LASERGENE 7.1.0 (DNAStar Inc., Madison, WI, USA).

Data analysis

Genetic diversity

All the sequences were compared, and haplotypes identified using DnaSP 5.10.1(Librado & Rozas, 2009). Further, we calculated genetic diversity across all studied populations in terms of Number of haplotypes (nHT), Haplotype diversity (HTdiv), Nucleotide diversity (ndiv) together with their respective standard deviations using Arlequin version 3.5.2.2 (Excoffier & Lischer, 2010).

Phylogenetic tree analyses and haplotype network

We investigated the association of Nigerian Muscovy duck with other Muscovy duck samples retrieved from public database (Table S1). MEGA version 6.06 (Tamura et al., 2013) was used to construct a rooted neighbor-joining (NJ) tree (Saitou & Nei, 1987) with the Maximum Composite Likelihood evolutionary distance approach (Tamura, Nei & Kumar, 2004). The bootstrap test was employed at 1,000 replications to assess the confidence of each node (Felsenstein, 1985). Further, we visualized the relationships between haplotypes by constructing median-joining networks of 397 Muscovy duck cytochrome b (rooted with chicken sequence; GenBank: KM886937) and 396 Muscovy duck nuclear DNA CYP2U1 sequences. Both networks were constructed using the median-joining (MJ) algorithm (Bandelt, Forster & Röhl, 1999) implemented in Network 5.0.1.1 (www.fluxus-engineering.com). The Median networks with all possible short trees were simplified by running the maximum parsimony (MP) calculation option to eliminate superfluous nodes and links.

Demographic dynamic profiles and population genetic structure

The demographic patterns and population dynamics, demographic statistical parameters for Tajima’s D (Tajima, 1989), Fu’s Fs (Fu, 1997) tests, and population F_ST of Nigerian Muscovy duck population were calculated using Arlequin version 3.5.2.2 (Excoffier & Lischer, 2010). The relationship between geographic distance and genetic distance was analyzed using linear regression model (LM) carried out under R statistical environment (R Core Team, 2021). Further, mismatch distribution patterns were estimated (Rogers & Harpending, 1992) with their geographical regions. We inferred genetic variation within populations, among populations, and groups of populations, analysis of molecular variance (AMOVA) with 50,000 permutations in Arlequin version 3.5.2.2 software. This analysis was done for Nigerian Muscovy duck populations at different hierarchical levels using the F_ST parameter at a significant P level of 0.05.

Cytochrome b genetic diversity

The 397 Muscovy duck cytochrome b sequences from Nigeria and India classify into 42 haplotypes (Table S2). The 319 Nigerian individuals (248 de novo (GenBank: MZ383528–MZ383776) and 71 from the GenBank) constitute 31 distinct haplotypes (Table S3), resulting from 57 polymorphic sites. The first most frequent haplotype, H4 (ADE15) occurs only in Nigerian samples while the second, H8 (BAA106) appears both in Nigerian and Indian populations (Table S2). The estimated haplotype diversity (HTdiv) across all Nigerian individuals is 0.266 ± 0.033, which is lower than the Indian estimate in almost three folds (Table 1). Within Nigeria, fourteen haplotypes occur in the Kwara population with the diversity value 0.458 ± 0.083, while three in the Oyo population with a diversity value of 0.060 ± 0.040 (Table 1). There were differences among populations with respect to their nucleotide compositions.

Cytochrome b gene Phylogenetic tree and haplotype network

A Phylogenetic tree was constructed using 71 Muscovy duck sequences selected from haplotypes in Table S2. The 71 sequences are Nigerian Muscovy (n = 42) newly generated in this study together with other samples from Nigeria (n = 16), India (n = 11), the USA (n = 1), and China (n = 1) retrieved from the GenBank (Table S1). Chicken sequence (GenBank: KM886937) was used to root the phylogenetic tree. The phylogenetic tree showed significant clustering of Nigerian Muscovy duck together with Indian, USA, and Chinese Muscovy. A few minor clustering was also noticed within the major clade (Fig. S1).

The median-joining network of 397 Muscovy duck together with chicken cytochrome b revealed 42 different haplotypes, most of which were singleton (Fig. 1B (i)). Of the 42 haplotypes observed, H4 (ADE15) showed high frequency, which was present in all eight populations. The next prominent haplotype, H8 (BAA106) was present in five Nigerian Muscovy duck populations studied and shared values with Indian populations (Table S2). H4 (ADE15) exhibits a star-like pattern indicating rapid population expansion (Fig. 1B (i)).

Cytochrome b population genetic structure and historical demographics

The Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) are −2.409 (P < 0.05) and −25.209 (P < 0.05), respectively (Table 1). These values also signal demographic expansion. The extent of haplotype-sharing in the network indicates the absence of a definite population structure in Nigeria Muscovy duck.

Mismatch distribution and pairwise F_ST were calculated to infer the population dynamics and maternal genetic structure of Nigerian Muscovy duck. Mismatch distribution patterns were unimodal for Nigerian Muscovy duck populations (see Fig. S3), which is in accordance with the signal of demographic expansion. The F_ST distance between most Nigerian Muscovy duck populations studied was low (Table 2). Among the eight populations, the highest genetic distance (F_ST) was observed between Oyo and Niger populations (Fig. 1A and Table 2), whereas the lowest (negative F_ST) were observed between Taraba and the six populations except Ekiti, Niger and Sokoto, Bayelsa and Kano, Oyo and Ekiti. Negative F_ST values is effectively seen as zero values. This indicates absence of genetic subdivision between the populations. The negative F_ST may possibly be product of sampling size and low haplotype diversity. The regression analysis results showed a negative correlation between geographic distance and genetic distance (R² = 0.2268, P < 0.05), suggesting the possibility of “absence” or “insignificant change” in genetic distance associated with increase in geographic distance (Fig. 2A).

For the analysis of molecular variance (AMOVA) among populations, F_ST was significant (F_ST = 0.02550, P = 0.01075; Table S6). Other analyses to assess differentiation within Nigeria were statistically non-significant (P > 0.05; Table S6).

Nuclear DNA CYP2U1 gene genetic diversity

The nuclear DNA CYP2U1 gene sequences of 396 Nigerian and Indian Muscovy ducks classify into 95 haplotypes, and five haplotypes were common in both countries (Table S4). The 316 Nigerian individuals (237 de novo [GenBank: MZ383202–MZ383439] and 79 from the GenBank) constitute 91 distinct haplotypes (Table S5), resulting from 84 variable sites. The most frequent haplotype, H1 (ABO14) occurs only in Nigerian samples while the second, H82 (MUSCO02), occurs in Nigerian and Indian populations (Table S4). The estimated haplotype diversity (HTdiv) across all Nigerian individuals is 0.800 ± 0.023, which is higher than the Indian estimate in almost one and half folds (Table 1). Within Nigeria, 32 haplotypes occur in the Oyo population with a diversity value of 0.869 ± 0.037, while two in the Taraba population with a diversity value of 0.667 ± 0.314 (Table 1).

Nuclear DNA CYP2U1 gene phylogenetic tree and haplotype network

The phylogenetic tree was constructed using 168 Muscovy duck nuclear DNA CYP2U1 sequences selected from haplotypes in Table S4. This comprises of Nigerian Muscovy duck sequences newly generated in this study together with other samples from India retrieved from the GenBank. The phylogenetic tree showed a major clustering of Nigerian Muscovy duck together with Indian Muscovy. A few minor clustering was also noticed within the major clade (Fig. S2).

The haplotype network of 396 Muscovy duck nuclear DNA CYP2U1 gene revealed 95 different haplotypes (Fig. 1B (ii)). Of the 95 haplotypes, H1 (ABO14) showed high frequency and was present in all eight populations (Table S4). The next prominent haplotype, H82 (MUSCO02), is dominated by Indian Muscovy duck populations and is also shared with Nigerian Muscovy duck from Oyo (Table S4). Some Nigerian populations showed the star-like pattern in Fig. 1B (ii), indicating rapid population expansion. The Tajima’s D (Tajima, 1989) and Fu’s Fs (Fu, 1997) are −0.707 (P = 0.249) and −5.228 (P < 0.05), respectively (Table 1). These values also signal demographic expansion. The extent of haplotype-sharing in the network indicates the absence of population structure in Nigeria.

Nuclear DNA CYP2U1 gene population genetic structure and historical demographics

Mismatch distribution and pairwise F_ST were calculated to infer the population dynamics and maternal genetic structure of Nigerian Muscovy duck. Mismatch distribution patterns were unimodal for Nigerian Muscovy duck populations (Fig. S3), which is in accordance with the signal of demographic expansion. The F_ST distance between most Nigerian Muscovy duck sampling localities was low (Table 2). Among the eight populations, the highest genetic distance (F_ST) was observed between Ekiti and Sokoto populations (Table 2), whereas the lowest (negative F_ST) were observed between Taraba and the six populations except Ekiti. The negative F_ST may possibly be product of sampling size and low haplotype diversity. The regression analysis results (R² = 0.2666, P < 0.005) also showed consistency with cytochrome b gene (Fig. 2B).

For the analysis of molecular variance (AMOVA) among populations, F_ST was significant (F_ST = 0.02318, P = 0.00978; Table S7). Other analyses to assess differentiation within Nigeria were statistically non-significant (P > 0.05; Table S7). Our findings with nuclear DNA CYP2U1 support the general absence of genetic structure in Nigeria as earlier revealed by the cytochrome b gene.