Development and parentage analysis of SNP markers for Chestnut-vented Nuthatch (Sitta nagaensis) based on ddRAD-seq data

1.   Introduction

Monogamy is one of the most common mating systems in birds, with 92% of species exhibiting monogamous behavior (Griffith et al., 2002; Reitsma et al., 2018). With the application of parentage testing techniques in ornithology, researchers have found that extra-pair copulations (EPCs) are prevalent in monogamous birds, which poses a new challenge to the traditional concept of sexual selection and mating systems (Griffith et al., 2002; Cockburn, 2006; Jetz and Rubenstein, 2011). EPCs refer to the phenomenon of females mating with males outside of their social pair bond. This behavior can lead to extra-pair fertilities (EPFs) and ultimately result in the production of extra-pair offspring (EPO), which is referred to as extra-pair paternity (EPP) (Birkhead, 1987; Wan et al., 2013). Birds have played a crucial role in the study of EPP across various animal mating systems (Kaiser et al., 2017).

Among socially monogamous passerine species, nearly 76% of species exhibit EPP, accounting for 11% of offspring (Griffith et al., 2002; Brouwer and Griffith, 2019; Cousseau et al., 2020). In previous studies, parentage analysis was largely examined by microsatellites. Microsatellites have been the most widely used molecular markers for studying parentage relationships in wild populations of fish, reptiles, birds, and mammals (Vignal et al., 2002; Uller and Olsson, 2008; Coleman and Jones, 2011; Dawson et al., 2013). The high mutation rates of microsatellite markers result in a high level of polymorphism, allowing the identification of individuals and estimation of kinship with only a few loci (Kaiser et al., 2017). However, there are still some problems with using microsatellites for parentage analysis. For example, microsatellites have high genotyping errors and are difficult to score accurately (Pompanon et al., 2005; Kalinowski et al., 2007). It is extremely successful in species with abundant and highly polymorphic genomes (DeWoody and Avise, 2000); however, the lack of polymorphism in some species makes robust parentage analysis difficult (Flanagan and Jones, 2019).

Thus, the availability of genomes and the emergence of methods to simultaneously discover multiple single nucleotide polymorphisms (SNPs) have led researchers to shift from using microsatellites to SNPs for parentage analysis. Compared with microsatellites, SNPs have the advantages of broader genome coverage (Brumfield et al., 2003), high-throughput, known mutational processes, and low-cost genotyping (Brumfield et al., 2003; Ellegren, 2004; Anderson and Garza, 2006). Although most SNPs are biallelic with relatively lower polymorphism and heterozygosity, which weakens the potential of each SNP loci for parentage analysis, this problem can be alleviated by increasing the number of SNPs involved in the analysis (Glaubitz et al., 2003). Selecting the loci with higher minor allele frequency (MAF) can enhance the parentage exclusion ability of markers (Dussault and Boulding, 2018). Many studies have demonstrated that SNP markers also provide a more accurate method for parentage analysis (Flanagan and Jones, 2019). A significant number of studies comparing the efficacy of SNPs and microsatellites for parentage analysis indicate that a relatively small number of variable SNP markers generally provide an equivalent or superior level of accuracy in parentage assignment compared to the available microsatellite markers (Kaiser et al., 2017; Thrasher et al., 2018; Weng et al., 2021). Furthermore, SNPs are considered more efficient in parentage analysis due to their high reproducibility and low genotyping error rates (Kaiser et al., 2017; Flanagan and Jones, 2019).

Double-digest restriction site-associated DNA sequencing (ddRAD-seq) is a RAD-seq protocol that could provide a more general approach for SNP discovery by selecting a smaller fraction of the genome (Peterson et al., 2012; Kess et al., 2016). Thrasher et al. (2018) described a modified ddRAD-seq method that can simultaneously discover and screen a large number of SNPs with high power. By applying this method to analyze parentage and relatedness based on SNPs in the socially complex and highly promiscuous Variegated Fairy-wren (Malurus lamberti), it was demonstrated that the ddRAD-seq method successfully recovered a substantial number of SNP loci, enabling confident determination of relationships within a species in a complex social system. The number of reports on parentage analysis using ddRAD-seq is increasing, and some studies have demonstrated that SNPs identified through this method are superior to highly polymorphic and species-specific microsatellite loci in terms of accurately assigning paternity and estimating relatedness (Thrasher et al., 2018; Miller et al., 2022).

The Chestnut-vented Nuthatch (Sitta nagaensis) (Passeriformes, Sittidae) is a typical secondary cavity-nesting species that inhabits coniferous and mixed-coniferous forests at altitudes of 1500–3000 m (John et al., 2000; Mo et al., 2023). It is a socially monogamous bird and both sexes provide biparental care during breeding. Paternity analysis of Sittidae species has been studied only in Eurasian Nuthatch (S. europaea) and Brown-headed Nuthatch (S. pusilla).

S. nagaensis is a closely related species of S. europaea and S. pusilla. Given that EPP has been found in both species, we hypothesize that S. nagaensis also exhibits phenomenon of EPP. For this study, we developed SNP markers from the population of S. nagaensis using ddRAD sequencing. Individuals were genotyped with each set of markers and parentage tests were performed. The objectives of this study were to (1) identify viable SNPs for parentage test in S. nagaensis and assess the effectiveness of the resulting marker panel in accurately determining parentage relationships; (2) verify whether phenomenon of EPP exists in this species. Our study will contribute to understanding the mating system of S. nagaensis. More generally, it provides a case study for improving animal mating system theory.

2.   Materials and methods

2.1   Study area and data collection

The study area is located in Zixi Mountain Provincial Nature Reserve of Yunnan Province, China (24°58′58″–25°04′01″ N, 101°22′49″–101°26′07″ E), with a total area of 160 km² (Fig. 1) (Mo et al., 2023). This area is located in the monsoon belt of the northern subtropical plateau and has an average annual temperature ranging from 12.1 to 14.9 ℃ and receives an annual average precipitation of ~900 mm (Feng et al., 2019; Mo et al., 2023). The main vegetation types in the reserve are secondary coniferous forest and coniferous and broad-leaved mixed forest. Castanopsis kawakamii, Pinus yunnanensis, Pinus armandii, and Alnus nepalensis forests are the main forest types (Mo et al., 2023).

Figure  1.  Overview of the study area: (A) map of Yunnan Province; (B) Chuxiong Yi Nationality Autonomous Prefecture in Yunnan Province; (C) Overview of the study site and the layout of artificial nest boxes.

DownLoad: Full-Size Img PowerPoint

Three hundred artificial nest boxes were established from 2021 to 2023 to attract S. nagaensis. The nest boxes were made from pine wood and were 25 cm × 14 cm × 12 cm (height × length × width) in size, with a 5 cm diameter hole. They were installed on trees approximately 2–4 m above the ground, with random orientation and on random tree species. The linear distance between each nest box was more than 50 m (Zhang et al., 2021; Lecce et al., 2023), and the locations of nests were recorded using GPS. We checked the nest boxes every 1–3 days from March to June each year. When S. nagaensis bred in the nest boxes, the clutch size and nestling size were monitored. Adults were captured using mist nets when the nestlings at the age of 6–8 days old, and the nestlings were captured at 10–12 days old (Wang et al., 2021; Arrieta et al., 2022). We marked all individuals with colored plastic rings and assigned corresponding numbers for identification. We also collected a blood sample (10–15 μL) from all individuals for parentage analysis. Blood samples were collected from the brachial vein and stored in ethanol until DNA extraction. After sampling, all birds were released at the capture sites. A total of 14 males, 11 females, and 40 offspring from 13 nests were captured in three years. However, we did not collect information from nests of three female adults, resulting in a lack of information on relevant social mothers.

2.2   DNA extraction and ddRAD sequencing

Total genomic DNA was extracted from the blood samples using the Trelief ® Hi-Pure Animal Genomic DNA Kit (TSP202, TSINGKE, Beijing, China) according to the manufacturer's protocol. To assess the integrity of DNA, 1% agarose gel electrophoresis was performed, and the purity and concentration of DNA were determined using a NanoDrop 2000 (NanoDrop Technologies, Wilmington, DE, USA). After the completion of DNA extraction, we used the sex identification primer (P2/P8) to determine the gender of all adult birds (Valenzuela-Guerra et al., 2013).

There was no reference sequence available for S. nagaensis, so the ddRAD sequencing method was taken since this is suitable for the development of SNP loci in non-reference genome species (Baird et al., 2008; Gutierrez et al., 2017). This method could simultaneously discover and screen a large number of SNP loci, providing high power for addressing questions related to parentage testing (Thrasher et al., 2018). To evaluate the dependability of our SNP panel for parentage analysis, we selected one individual from any individual nest in all the years of our study. The DNA samples from the first two years were difficult to sequence due to prolonged storage time so, in rare cases, we selected two individuals from one nest. A total of 10 qualified DNA samples of S. nagaensis were selected for ddRAD sequencing by Shanghai Personal Biotechnology Co., Ltd. (Shanghai, China), to develop SNP markers. The sequencing samples were obtained from 1 nest (1 individual) in 2021, 3 nests (4 individuals) in 2022, and 4 nests (5 individuals) in 2023 (Table 1).

Table  1.  A list of 10 individuals of S. nagaensis based on ddRAD-seq.

Year	Nest box No.	Individual No.
2021	E09	2021E0922
2022	G29	2022G2921
	G29	2022G2922
	G51	2022G5111
	G66	2022G6612
2023	E09	2023E0921
	E86	2023E8621
	H91	2023H9124
	G75	2023G7511
	G75	2023G7512

 | Show Table

DownLoad: CSV

The extracted Genomic DNA from 10 individuals was digested with a restriction enzyme, 5 μL of each sample was evaluated for agarose gel electrophoresis, and then the oversized and undersized fragments were removed through magnetic bead capture (VAHTSTM DNA clean beads). Adapters P1 and P2 were ligated to the fragments, where adapter P1 contained the amplification primer site, sample barcode, and restriction enzyme site, and adapter P2 contained the amplification primer site and restriction enzyme site (Wang et al., 2017). PCR amplification was performed on DNA fragments connected with P1- and P2-specific primers to enrich sequencing library templates. Then, the PCR product was purified using magnetic bead capture. Before sequencing, the library was subjected to quality inspection on the Agilent Bioanalyzer using the Agilent High Sensitivity DNA Kit (Personalbio, Shanghai, China). After the qualified sequencing libraries were gradient-diluted, the 2 × 150 bp paired-end sequencing workflow was performed on the Illumina NovaSeq sequencer (Personalbio, Shanghai, China).

2.3   Development and genotyping of SNP markers

Raw sequencing data was filtered using fastp v0.20.0 (https://github.com/OpenGene/fastp) (Thrasher et al., 2018) based on the sliding windows algorithm and the strict criteria listed below were used to eliminate low-quality reads from the analysis: (1) removed reads with adapter contamination from 3' end; (2) removed reads with a Q value < 20; (3) removed reads if the length of any read was ≤50 bp in the paired-end; (4) removed reads if the number of unidentified nucleotides (N) was ≥5 in the paired-end (Thrasher et al., 2018; Weng et al., 2021). Reads in a single sample were clustered according to sequence similarity using the "ustacks" command in the Stacks version 2.55 package (Catchen et al., 2013). We set the minimum stack depth of 3 (m parameter) and the maximum stack difference of 3 (M parameter). When building the catalog, a maximum of 4 mismatches (n parameter) were allowed between loci of different individuals. The set of parameters mentioned has demonstrated applicability for other polymorphic passerine species (Campagna et al., 2015), but the optimal combination of parameters may differ for various datasets (Thrasher et al., 2018). SNPs obtained in Stacks are filtered through the populations' module and outputted. Our approach involved grouping all samples within a population, and we outputted a locus only if it was present in at least 50% of the individuals within that population (r parameter) and possessed a minor allele frequency of at least 0.05 (–min_maf). To avoid objective errors affecting the accuracy of the experiment during sequencing, SNPs were only searched for in differential contigs containing at least 30 reads, and SNPs with more than a 30% difference between samples were selected.

We randomly selected 100 SNP markers based on filtered data and designed SNP primers using Premier 5.0 (Lalitha, 2000). Then we performed three rounds of PCR amplification screening on SNP primers. For the first screening, we randomly selected 12 DNA samples of S. nagaensis and set six annealing temperature gradients (55 ℃, 57 ℃, 59 ℃, 61 ℃, 63 ℃, 65 ℃), and then PCR amplification was performed on the six gradients of each sample using 100 pairs of primers. The PCR was performed in 25 μL of reaction volume containing 12.5 μL of 2 × T5 Super PCR Mix (Tsingke, Beijing, China), 1 μL of each primer (10 μmol/L, synthesized by Tsingke, Beijing, China), 1–2 μL of template DNA, and added ddH₂O up to 25 μL. The amplification reaction was performed by initial denaturation at 98 ℃ for 3 min, followed by 35 cycles of denaturing at 98 ℃ for 10 s, annealing temperature for 10 s, and extending at 72 ℃ for 15 s and a final extension at 72 ℃ for 1 min. The PCR products were detected by 2% agarose gel electrophoresis. To further screen the primers with higher versatility, 8 DNA samples were randomly selected and PCR amplified by the SNP primers obtained from the preliminary screening. The PCR reaction system and PCR conditions (annealing temperature set at 59 ℃) were the same as the first screening. The PCR products were detected by 2% agarose gel electrophoresis. For the third screening of SNP primers, the primers obtained from the second screening were used for PCR amplification in all individuals. The total volume of PCR reaction and reaction conditions were the same as those of the second screening. The PCR products were detected by 1.5% agarose gel electrophoresis, and 54 loci were successfully amplified (Table 2).

Table  2.  Information of SNP loci used in this study.

Locus	PCR amplification forward primer	PCR amplification reverse primer	Variation	Tm (℃)
LTSNP1	GGGACTCTCCTTCCTGACCT	AGGGACTGGGTGAGAGGATT	G/T	59
LTSNP4	GCCCATTTTGTGTGAACCTC	GCACATTAACATGGACAAAAGC	G/A	59
LTSNP5	AGGAGGATGCTCTGAAGCTG	CAGGGATGGTCACTCCAACT	T/C	59
LTSNP6	ATGAGAAGCCCTGGCCC	ACAGTGTCACAGCCCTGCTT	C/T	59
LTSNP14	CAGCTGGAGGAGACCTTGTG	GCCTGAACACCTGTTGTGTG	C/T	59
LTSNP15	AGGCTGTGTTTCTCAGTCCC	GCCTTTTGAAGGAGAGTGTTCA	C/G	59
LTSNP17	GGCCCCAGATTAATTTCAGA	CTGGTGGCAGAAGGTCTCTC	C/A	59
LTSNP20	GAGTGGTGCCACAGGTGC	CTGCAGATGCTGTCCCCT	A/C	59
LTSNP24	CTGCAGAGGGATGTTTCCA	GGCAGATGCTGTGGGAAT	A/T	59
LTSNP25	GCTCTTCCACAAGGTCCAAC	CCTGCAGGCAAAAAGAAAAG	A/C	59
LTSNP27	CAGTGCAAGTGCAGCAGATT	TTTCCTGACAAGTCTCTCTCCTG	T/C	59
LTSNP31	GGGACACAGGGGTGAGC	GGCTAGCTGGCACAGAGTTC	G/T	59
LTSNP34	AGCTTTGAACACTCTCCTCCC	TTCCTGTGAACTCCAAGCCT	A/C	59
LTSNP35	GCAAGGTGGATTCTCCATCA	TAAATACACAGCCACGTCCC	C/G	59
LTSNP36	GCTCCAGGCACACACATTTA	AGCTCTTCAGCTTCATGGGA	A/G	59
LTSNP37	CTGACAATGCCTCCCCAAG	TAAGAATCCCCCTCTGGCTT	G/A	59
LTSNP39	AGCCACCTGTTACCTCCGT	AAAGTCCTGGAAATCCTCCC	G/C	59
LTSNP40	GGGGTCTTCTCTGAGGCAAC	TCTCCATCTCCCTGAGCAAC	G/A	59
LTSNP41	TGCTGTGCTAACTAGACTGACCA	AACCAAGCTGGGAAGAGGTT	C/T	59
LTSNP43	ATGATGCCACTCCATGAACA	GACCACAAGGACTGTGGCTT	T/G	59
LTSNP47	TTTCTAAAATTTGCGGGCAC	ACCTGCACCAAGCAGGC	A/T	59
LTSNP48	GAGGATGCTGCTGGGAACT	CCACACTGGTGACAGCTTTG	T/C	59
LTSNP49	TGTGTCCCCTGGCATTTTA	GGCAATGTCAAGGCAGTTCT	T/C	59
LTSNP50	ACTGCACCAACGGGGAA	TACTTGTGTTTGGCAAGGCA	C/T	59
LTSNP51	GGAAGCGTGAGTATAGGCCA	GGGAACTGCACTACTGAGGG	G/A	59
LTSNP52	GGCCTGTGGCTTACTTTCTC	TCTTTGTCACCTCATGCAGC	A/G	59
LTSNP53	AGAGCCAGCACTCCCAACT	CACAAAATCCCTAGGAGGGC	C/A	59
LTSNP54	CTTCTTGGGAAGAAAGCTCG	CACGTGAGCCAACTCTCATT	G/A	59
LTSNP55	GGAAACTTCCTGTGCCTGAC	CCTCCCTGACTCACCTGGAG	G/C	59
LTSNP56	ACCTCGTTTTTCCCACTTCC	AGCAGGCCTGAGAATAGCAA	A/G	59
LTSNP58	TGCCTCTCACACTAGAGGGG	CCACAGCACTTGCCCTTTAT	C/T	59
LTSNP60	TGCCATAGCAGATTGCCG	TTCACGCACTGACCTTCTTG	T/C	59
LTSNP62	TCAGCTTGATGGTCAGCAGT	CTTGCCATTGCCATCTAACC	C/T	59
LTSNP63	GGGCTTAGAGCACTCCCAAT	GCAGAAGTGAATGACACCAACT	A/T	59
LTSNP65	CCCTGGCTTTGTAGGGGAG	ATTTGTAGTCGTTGACCCGC	C/T	59
LTSNP67	CTCTGCCGCTCCATTTCTAC	GGAGTTATGCAGAGCACAACAA	C/A	59
LTSNP68	GCAGAGCAGAGCAGTTTGTG	TGAGGCTTCAGCTGAGTTTG	T/G	59
LTSNP69	AGTTTAACCCCTCTCCCCAC	ATCAGCTTGAGCTGTTGGGT	C/T	59
LTSNP70	GGCTCTGTGCTATTCCCCTC	ATTCTCAGGGCTTGGATGG	T/C	59
LTSNP71	AGGGTACATTGTTCCTGGGG	CCCAGTAGAGAGGGGTGGTG	A/G	59
LTSNP72	TACGTGACCCCAATCTCCTG	ACGTGCCCAGGATTAACACA	G/C	59
LTSNP73	CCAAAGGCAAACCACAAAAT	CCTCAGAAGTCCACTCATCCA	A/G	59
LTSNP74	TGTGCCTTTCCCTCCTAAAA	TCTAAGTACGCAATGCAGCG	G/C	59
LTSNP76	GGGAGCATTTATTGCAGCAT	CCCCTGAAGTGACAACCACT	C/T	59
LTSNP77	GGGGACAGTGAGTTCTGCTG	TGGTGAGTGTCCATGAGAGC	G/A	59
LTSNP78	TGGTAGGAAACCACTGGAGG	AGGGCTAGTTGGATGGCTTT	T/A	59
LTSNP81	TGAGTTGACAAGTGACAGCCT	GGACAGCAAAGCAGACTGAG	C/T	59
LTSNP82	GTGTGGGTGAATTTTGCTCA	CCCCTTGCCACTCTGTAAAT	T/A	59
LTSNP88	AAATGCCTAAAGAGCTGGGG	ACCACCGTTTGAGTAGCAGC	T/C	59
LTSNP90	ATTTTGTCAATCCGCTCCTG	AATGAGGGGAATTTTCTGGC	C/A	59
LTSNP92	CTACCCAAGCATTTGCCGT	AAATTTTATTCGGCTGGGCT	C/G	59
LTSNP94	CATTCACCTCCCTCCCTG	ATCCATGGTCAGTCCTCTGC	T/A	59
LTSNP95	ACACAGCCCCCGTCCTC	GAGACGTGTCCTGCCAGC	T/C	59
LTSNP99	TGATGGGTGGAGCTGGAG	CCCAGTGTGGGAAGAGATTC	C/A	59

 | Show Table

DownLoad: CSV

The PCR products of all individuals were sent to Tsingke Biotechnology Co., Ltd. (Beijing, China) for direct sequencing. We used SnapGene software (Version 5.1.5, Applied Biosystem) to determine the genotype according to the corresponding base of the SNP loci. To eliminate genotyping errors, loci that were unable to be used for matching offspring with their parents were subjected to re-PCR amplification and re-genotyped.

2.4   Parentage analysis

After three screenings, the SNP marker panel set for parentage analysis of S. nagaensis contained 54 markers, developed from 100 SNPs originally developed and supplied. We first examined Hardy-Weinberg equilibrium (HWE) and Linkage disequilibrium (LD) at 54 loci using PLINK v1.90 software (Purcell et al., 2007). SNPs deviated from HWE with a p-value < 0.05 (Weng et al., 2021). The LD was obtained by calculating r² values for every pairwise combination of the 54 SNPs in the final dataset. We removed one locus from each pair of loci with r² > 0.2 (Weng et al., 2021). Then, we calculated minor allele frequency and removed loci with MAF < 0.05. In addition, we calculated the polymorphism parameters including allele number (N), observed heterozygosity (H_o), expected heterozygosity (H_e), polymorphic information content (PIC), and cumulative exclusion probability (Excl) using CERVUS version 3.0.7 (Kalinowski et al., 2007).

We conducted population parentage analysis on the combined samples collected over three years, using CERVUS version 3.0.7 for parentage assignments with SNP genotype data for all nestlings. This method assigned maternity or paternity using the natural logarithm of the likelihood ratio (LOD score), providing the likelihood of maternity or paternity of each candidate parent relative to a random parent in the population (Marshall et al., 1998; Kalinowski et al., 2007). To determine the threshold log-likelihood scores for assigning parental pairs with known sex, we employed the simulation module. This was done by simulating 10, 000 offspring with 11 candidate mothers and 14 candidate fathers. The proportion of gene typing is determined by the frequency analysis of alleles. We used a default genotype error rate of 1% and default levels of 80% (relaxed) and 95% (strict) confidence (Liu et al., 2015; Weng et al., 2021). The CERVUS assignment was considered acceptable if the highest-ranked parent demonstrated a pair LOD score that was positive or slightly negative, and the number of genotype mismatches between the assigned parent and the offspring was ≤8 (Arrieta et al., 2022). To match both parents, each offspring was tested for all possible maternity and paternity combinations. For the offspring of three nests lacking female adults, all female adults were selected as candidate mothers, and possible maternity testing was also conducted. The remaining unresolved offspring were matched with single parents. In addition, we compared the parentage assignment results with those recorded in the field.

3.   Results

3.1   SNP marker characteristics and identification power

After ddRAD sequencing, a total of 78, 476, 448 reads were obtained from 10 samples, producing about 10.81 Gb of high-quality data. The selection process for parentage verification primarily relied on assessing the quality of genotyping and the informativeness of SNPs. Two markers were homozygous (LTSNP88, LTSNP99), possibly due to the small sample sizes of S. nagaensis—the remaining 52 markers were heterozygous. After removing the loci significantly deviated from HWE and the loci in LD (LTSNP71, LTSNP99) and MAF < 0.05 (LTSNP67, LTSNP88, LTSNP99), we retained 50 SNPs for further analyses (Table 3).

Table  3.  Information of genetic indices of the 54 SNP markers developed in S. nagaensis.

Locus	Position	Ho	He	PIC	P(HWE)	F(Null)	MAF
LTSNP1	18, 703	0.477	0.410	0.324	0.625	−0.079	0.05
LTSNP4	62, 612	0.338	0.504	0.375	0.042	0.193	0.32
LTSNP5	142, 062	0.200	0.349	0.287	0.286	0.268	0.12
LTSNP6	206, 612	0.108	0.456	0.350	0.056	0.616	0.29
LTSNP14	789, 062	0.508	0.481	0.363	0.677	−0.031	0.14
LTSNP15	811, 009	0.231	0.382	0.307	0.056	0.243	0.14
LTSNP17	952, 451	0.523	0.451	0.348	0.663	−0.078	0.08
LTSNP20	1, 362, 742	0.338	0.484	0.365	0.421	0.173	0.23
LTSNP24	2, 221, 852	0.538	0.500	0.373	0.691	−0.041	0.18
LTSNP25	2, 633, 639	0.338	0.440	0.348	0.640	0.132	0.14
LTSNP27	4, 112, 627	0.138	0.273	0.234	0.083	0.323	0.09
LTSNP31	8, 878, 820	0.338	0.389	0.312	0.286	0.066	0.09
LTSNP34	10, 250, 983	0.508	0.500	0.373	1.000	−0.012	0.20
LTSNP35	11, 523, 630	0.169	0.313	0.262	0.287	0.295	0.11
LTSNP36	12, 753, 698	0.492	0.451	0.348	1.000	−0.047	0.09
LTSNP37	13, 457, 669	0.169	0.273	0.234	0.056	0.231	0.08
LTSNP39	15, 684, 561	0.523	0.503	0.374	1.000	−0.024	0.22
LTSNP40	16, 187, 116	0.277	0.303	0.256	0.640	0.042	0.05
LTSNP41	16, 850, 242	0.400	0.429	0.335	1.000	0.032	0.11
LTSNP43	17, 500, 689	0.292	0.456	0.350	0.041	0.215	0.20
LTSNP47	19, 874, 297	0.246	0.389	0.312	0.101	0.222	0.14
LTSNP48	19, 915, 571	0.354	0.349	0.287	0.056	−0.010	0.05
LTSNP49	20, 569, 442	0.492	0.503	0.374	0.115	0.007	0.23
LTSNP50	23, 276, 657	0.108	0.456	0.350	0.083	0.616	0.29
LTSNP51	24, 372, 580	0.262	0.465	0.355	0.187	0.277	0.23
LTSNP52	25, 434, 035	0.262	0.456	0.350	0.187	0.268	0.22
LTSNP53	26, 695, 409	0.523	0.451	0.348	0.663	−0.078	0.08
LTSNP54	26, 771, 405	0.323	0.410	0.324	0.330	0.115	0.12
LTSNP55	27, 166, 256	0.369	0.374	0.302	1.000	0.003	0.06
LTSNP56	28, 836, 102	0.338	0.341	0.281	0.056	−0.001	0.05
LTSNP58	31, 402, 100	0.508	0.481	0.363	0.677	−0.031	0.14
LTSNP60	32, 773, 918	0.185	0.262	0.226	0.198	0.170	0.06
LTSNP62	34, 068, 575	0.108	0.456	0.350	0.101	0.616	0.29
LTSNP63	35, 245, 807	0.538	0.500	0.373	0.691	−0.041	0.18
LTSNP65	38, 186, 382	0.523	0.489	0.368	1.000	−0.037	0.15
LTSNP67	40, 697, 049	0.154	0.194	0.174	0.121	0.111	0.03
LTSNP68	43, 162, 206	0.262	0.435	0.339	0.083	0.246	0.18
LTSNP69	43, 175, 829	0.308	0.404	0.320	0.663	0.131	0.12
LTSNP70	43, 583, 831	0.323	0.349	0.287	0.586	0.035	0.06
LTSNP71	43, 809, 004	0.185	0.240	0.210	0.287	0.128	0.05
LTSNP72	43, 894, 360	0.308	0.341	0.281	0.586	0.047	0.06
LTSNP73	44, 217, 752	0.246	0.283	0.242	0.141	0.066	0.05
LTSNP74	44, 757, 469	0.446	0.500	0.373	0.688	0.053	0.23
LTSNP76	45, 049, 724	0.400	0.429	0.335	1.000	0.032	0.11
LTSNP77	45, 581, 157	0.262	0.465	0.355	0.187	0.277	0.23
LTSNP78	47, 025, 931	0.462	0.501	0.374	0.422	0.037	0.23
LTSNP81	52, 738, 243	0.308	0.404	0.320	0.663	0.131	0.12
LTSNP82	53, 716, 212	0.138	0.332	0.275	0.056	0.408	0.14
LTSNP88	57, 074, 644	0.000	0.000	0.000	1.000	1.000	0.00
LTSNP90	57, 541, 813	0.523	0.477	0.361	1.000	−0.050	0.12
LTSNP92	61, 666, 209	0.169	0.313	0.262	0.287	0.295	0.11
LTSNP94	62, 729, 015	0.462	0.501	0.374	0.422	0.037	0.23
LTSNP95	62, 767, 856	0.492	0.503	0.374	0.115	0.007	0.23
LTSNP99	78, 038, 803	0.000	0.031	0.030	0.020	0.427	0.02

 | Show Table

DownLoad: CSV

For the 50 variable loci, the average values of genetic variability parameters H_o, H_e, and PIC of the 50 markers were 0.343 (ranging from 0.108 to 0.538), 0.419 (0.262–0.504), and 0.327 (0.226–0.375), respectively. Four of the 50 SNP markers were low informative (PIC < 0.25) and the rest were polymorphic (0.25 < PIC < 0.5, moderate polymorphism) (Table 3). The minor allele frequency ranged from 0.05 to 0.32. When one parent was known, the cumulative exclusion probability of the combination was between 0.991 and 0.9999, and when both parents were known, the cumulative exclusion probability of the combination was as high as 0.999999.

3.2   Parentage assignment

The parents of 13 broods and their 40 nestlings were genotyped at 50 SNP loci (Table 4). For the maternity analysis using SNPs, 31 of 40 (77.5%) offspring were assigned with 95% confidence, 3 of 40 (7.5%) offspring were assigned with 80% confidence, and the remaining 6 of 40 (15%) offspring were not assigned a mother. In the paternity analysis, 36 of 40 (90%) offspring were assigned with 95% confidence, and 2 of 40 (5%) offspring with 80% confidence. However, 2 of 40 (5%) offspring were excluded from paternity assignment due to missing genotyping data or typing errors. In addition, for the offspring-identified maternal bond, all offspring were assigned mothers consistent with social mothers recorded in the field. The identification of paternity had different results. SNP panels assigned 3 out of 38 nestlings to males who were not their social fathers. That is, in 3 out of 13 broods (23.1%), we found one EPO per brood (Table 4), with the social father of 3 offspring (7.9%, a total of 38 offspring in paternity assignment) being excluded from paternity. The remaining offspring were assigned fathers consistent with the social fathers recorded in the field. Overall, no nestlings had pair-wise mismatches with their social mothers, but three nestlings had pair-wise mismatches with their social fathers.

Table  4.  The number of nests, adults, and nestlings in populations of S. nagaensis during the breeding seasons (2021–2023). Also, the number of EPO and the number of nests with EPO are given.

	Year			Total	EPO
	2021	2022	2023
Nests	4	3	6	13	3
Adults	7	6	12	25	–
Nestlings	11	10	19	40	3

 | Show Table

DownLoad: CSV

4.   Discussion

We have generated a new collection of SNP loci specific to S. nagaensis by performing ddRAD sequencing on 10 individuals, as no reference genome is currently available for this species. To our knowledge, this is the first study to develop SNPs based on the ddRAD-seq method for parentage analysis of Sittidae species. Previous studies have shown that SNPs obtained through the ddRAD-seq method possess considerably greater effectiveness in assigning paternity and estimating relatedness among individuals compared to a moderate number of microsatellite loci that are specific to the species (Thrasher et al., 2018).

We assessed the polymorphism of the used loci and the potential for parentage analysis. Most of the 50 SNPs (mean H_o = 0.343, mean PIC = 0.327) used in this study were moderately, not highly, polymorphic. This is understandable because SNPs are biallelic, meaning that they have a lower level of polymorphism per locus compared to multi-allelic microsatellites, and more SNPs help increase the discerning ability of individuals (Weng et al., 2021). The number of SNPs needed for parentage analysis is strongly dependent on the average expected heterozygosity of the panel and the observed heterozygosity (Morin et al., 2004; Weinman et al., 2015; Kaiser et al., 2017) and minor allele frequency of individual SNPs (Strucken et al., 2016). Through a series of simulations, Baruch and Weller (2008) found that as the number of SNPs increased while maintaining the same MAF cut-off, the exclusion probability also increased. Premachandra et al. (2019) also found the number of SNP loci utilized had a significant influence on the accuracy of parentage assignment, and the accuracy usually increased with increasing MAF. When the average MAF is low, a higher number of SNPs is required. By contrast, the MAF in our study was relatively low, possibly due to the small sample size, which may reduce the accuracy of using SNP loci for paternity analysis. Second, a major factor that reduces the accuracy of paternity assignments is the occurrence of missing data, leading to incomplete genotypes of offspring or candidate fathers. The presence of missing data presents a challenge in any application involving high-throughput SNP genotyping (Nielsen et al., 2011; Toews et al., 2015; Kaiser et al., 2017).

A higher degree of heterozygosity in the markers necessitates a smaller number of loci (Morin et al., 2004). The SNPs we selected for parentage assignments were relatively heterozygous (mean H_o = 0.343), and the paternity exclusion probability reached 99.9%, which could estimate the accuracy of parentage assignments. For the socially monogamous Setophaga caerulescens, Kaiser et al. (2017) demonstrated that through the selection of the most heterozygous SNPs, a panel of only 40 SNPs (with a mean Ho = 0.37) was enough to assign paternity to the same proportion of offspring as the complete panel of 97 SNPs. For EPP of Thryothorus pleurostictus, Cramer et al. (2011) showed that the power of exclusion of 41 SNPs (95%) was better than that of 7 microsatellites (88%), and the paternity exclusion probability of the combination of the two loci was 99.9%. If there were more SNPs available to fill the panels and the average minor allele frequency of the SNPs was closer to 0.50, it is expected that the difference in performance would have been even greater (Gudex et al., 2014).

We conducted parentage analyses for 40 nestlings from 13 nests and examined the results of parentage assignment based on field records. Our research results indicate that three offspring did not match their social father and were considered EPO, accounting for 7.9%. The EPP occurrence probability is 23.1%. In previous studies regarding the genus Sitta, Segelbacher et al. (2005) analyzed 32 nests of S. europaea using 5 microsatellites, finding 18 EPO, which accounted for 10% of all offspring. Notably, genetic analysis revealed that 12 of these EPO were sired by males from neighboring territories. In another study, Han et al. (2015) examined 59 nests of S. pusilla using 9 microsatellite loci, finding EPO in all 24 nests. Interestingly, two nests suggested that female helpers may have produced EPO. Compared to the two Sitta species, the proportion of EPP occurrence in S. nagaensis was lower than that in S. europaea (the probability of EPP was 37.5%; probability of EPO was 10%) and S. pusilla (the probability of EPP was 41%), and the proportion of EPO is also lower.

There are some limitations to our study. First, the sample size we collected over the past three years was relatively small, which is directly related to the low utilization rate of artificial nest boxes by S. nagaensis. The unknown risks associated with artificial nest boxes and human interference may be among the factors contributing to the low utilization rate of the nest boxes (Potti et al., 2021; Mo et al., 2023). Second, compared to other researchers studying avian parentage analysis, the MAF of the SNP markers we used for parentage testing was relatively low. Additionally, a few nests had incomplete data due to a lack of sample data from adult individuals, leading to incomplete nest records. The issues with data missing or typing errors in genotyping during parentage identification rendered accurate allocation of these individuals to their respective mothers or fathers impossible. In future research, we will address this issue by selecting more heterozygous SNPs. We expect to continue increasing the field samples to further clarify the proportion and mechanisms of EPP.

Ethics statement

The study was approved by the Ethics Committee of IACUC of Southwest Forestry University (IACUC-SWFU L20201220 and December 20, 2020).

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper..

CRediT authorship contribution statement

Qingmiao Yuan: Writing – original draft preparation, Methodology, Investigation, Conceptualization. Xi Lu: Investigation. Ruixin Mo: Investigation. Xianyin Xu: Investigation. Xu Luo: Writing – review & editing, Conceptualization. Yubao Duan: Writing – review & editing, Methodology, Conceptualization.

Acknowledgements

We are grateful to all those who have helped us with our field and experimental work. We also thank the reviewers for their valuable comments and suggestions to improve the manuscript.