Characterization of novel microsatellite markers of the Emei Shan Liocichla using restriction site-associated DNA sequencing

Background

Genetic diversity is crucial for the sustainability of species, especially for threatened species (Meffe and Carroll 1994). Habitat fragmentation, a common cause of endangering threatened species, can lead to reductions in the population size and gene flow among patches, which leads to a loss of genetic diversity and an increase in inbreeding (Frankham et al. 2002). The loss of genetic diversity associated with inbreeding has two important consequences: it reduces the ability of populations to adapt to environmental/climate change, and it can result in a direct and immediate loss in fitness through inbreeding depression (Keller and Waller 2002; Armstrong and Seddon 2008; Jamieson 2011). It is, therefore, necessary to assess the genetic diversity and inbreeding level, and to determine whether a population is genetically fragmented before undertaking the genetic management actions for threatened species, which are important in biodiversity conservation (Frankham et al. 2002).

Microsatellites have a higher mutation rate compared to mitochondrial DNA and nuclear introns (Ryman and Leimar 2008), a feature that can provide a powerful tool for detecting genetic differences, especially in small and fragmented populations (Li et al. 2002; Duncan et al. 2016). Microsatellite markers have been widely applied in molecular ecology and conservation genetic studies since the 1990s (Bruford and Wayne 1993; Jarne and Lagoda 1996; Chapuis et al. 2015; Faria et al. 2016). Although cross-species amplification works for a few closely related taxa (Dawson et al. 2010; Gu et al. 2012), the major drawback of microsatellites is that they need to be independently isolated from most species (Sudheer et al. 2010). Most microsatellites must be identified by screening genomic libraries with appropriate probes or enrichment protocols, which are inefficient and time-consuming (Zane et al. 2002). Whereas the Illumina next generation sequencing (NGS) platform, which can produce moderately long paired-end reads (up to 150 bp with the HiSeq 2500 platform), is an effective method for obtaining reliable microsatellites, even for birds that have few microsatellite loci (Castoe et al. 2012). The restriction-site-associated DNA (RAD) method (Miller et al. 2007) with Illumina sequencing (RAD-Seq) can further reduce the cost and improve the average depth per locus (Baird et al. 2008). Unlike other methods for generating genome-wide data, RAD-Seq does not require any prior genomic information of the taxa being studied (Andrews et al. 2016). To date, several studies have shown that paired-end RAD-Seq is a powerful tool for obtaining microsatellites in both model and non-model species (Miller et al. 2007; Pfender et al. 2011).

The vulnerable (VU) Emei Shan Liocichla (Liocichla omeiensis) is an endemic passerine bird found on Daxiangling Mountain and Xiaoliangshan Mountain in southwestern China. It lives on the edges or gaps of natural broadleaf forests between an elevation of 1450 and 2150 m (Fu et al. 2011). This bird is resident and moves vertically between seasons (Fu and Zhang 2011). Additionally, the Emei Shan Liocichla is socially monogamous and protects its territory continually throughout the year, with a territory range of approximately 10, 000 m² (Fu, unpublished data). However, the population of the Emei Shan Liocichla is suspected to be declining at a moderate rate, which is in line with the accelerated habitat loss and degradation within its range (BirdLife International 2017). The natural broadleaf forests in its habitat have been cleared or replaced by non-native coniferous forests and tea plantations, which might lead to habitat fragmentation as a previous study has shown the Emei Shan Liocichla did not use these types of vegetation (Fu et al. 2011). Habitat loss and fragmentation are regarded as the major factors contributing to the population decline in many bird species (Ribon et al. 2003; Gill 2007). In addition, during our sampling period, we witnessed a number of captures of these birds for sale, which will also damage the population. Although conservation actions have been taken since 1999, little is known about the genetic status of this vulnerable species and its populations. Therefore, to carry out genetic diagnoses before undertaking genetic management actions (Frankham et al. 2002), a set of suitable genetic markers, such as microsatellites, are urgently needed.

In this study, we utilized Illumina paired-end RAD-Seq to isolate highly polymorphic microsatellites of the Emei Shan Liocichla that can be useful for further studies, and tried to understand some preliminary genetic status to help the management and conservation of this vulnerable species using our incomplete sampling by answering the following questions: (1) how much genetic diversity is there in our sampling populations? (2) are the populations suffering from inbreeding? (3) are the populations from two mountains already genetically fragmented? and (4) what we can do for conservation of the Emei Shan Liocichla?

Methods

Sample collection and DNA purification

Fifty-one samples containing 38 blood and 13 tissue samples, were collected from six sites distributed in Daxiangling Mountain and Xiaoliangshan Mountain located in southwestern Sichuan, China (Fig. 1). Blood samples were collected from the brachial vein and stored in absolute ethanol. All the captured individuals were released immediately after the blood samples were obtained. Tissue samples were collected opportunistically from individuals that died naturally and dead chicks in nests. Permission for the sampling was granted by the Sichuan Forestry Department and the Laojunshan National Nature Reserve of Sichuan Province, China. The samples were frozen in the field and preserved at -80 ℃ in the lab. Total genomic DNA was extracted using TIANamp Genomic DNA Kit DP304-2 (Tiangen, China) following the manufacturer's recommendations.

Figure  1.  Sample locations of the Emei Shan Liocichla (Liocichla omeiensis)

DownLoad: Full-Size Img PowerPoint

RAD-Seq and microsatellite characterization

In total, 3 µg of genomic DNA from a single male was sent to Novogene Bioinformatics Technology Co., Ltd. (Beijing, China) for the DNA library preparation and sequencing to obtain sequencing data for microsatellite screening. We used P1 and P2 adaptors for the restriction enzyme sites to construct a DNA sequencing library with a size range of 300-700 bp. Paired-end 125 bp reads were obtained from separate lanes of the Hiseq 2500 Genome Analyzer. Poor-quality reads were cleaned and adapter reads were trimmed using Trimmomatic 0.36 (Bolger et al. 2014) with default settings. Reads containing RAD tags were clustered by different RAD tags using cd-hit-est (Li and Durbin 2009), and RAD tags containing 10‒400 reads were considered repetitive and removed. The remaining paired-end reads from each RAD site were then sent to Velvet (Zerbino and Birney 2008) to assemble contigs with the following parameters settings: VelvetOptimiser.pl -s 23 -e 31 -x 4. An SR search was developed by Novogene Bioinformatics Technology (Beijing, China) to detect the microsatellites from contigs using the following criteria: 2-6-bp motifs, longer than 12 bp, and the distance between two microsatellites was greater than 12 bp.

Primer design and microsatellite evaluation

We identified 11, 564 putative microsatellite sequences after excluding compound and interrupted microsatellites. We randomly designed 600 pairs of primers ranging in length from 12 to 36 bp using Primer Premier 5 (Lalitha 2000). The predicted PCR products included trimers, quadmers, pentamers and hexamers. PCR was performed using these primers on four separate Emei Shan Liocichla DNA samples. PCR was conducted in 10 μL reaction volumes containing ~20 ng of genomic DNA, 1 μL of 20 μM of each primer, and 3‒5 μL 2× Taq polymerase (Tiangen, China). Following a denaturation step of 5 min at 95 ℃, the PCR mixture was subjected to 30 cycles of 94 ℃ for 30 s, Tm for 30 s, and 72 ℃ for 45 s, followed by 72 ℃ for 10 min and storage at 4 ℃. A total of 99 primer pairs were chosen with clear straps, correct placement, and no non-specific amplification by agarose gel electrophoresis. Each locus was examined in eight Emei Shan Liocichla individuals by sequencing the PCR products (Tsingke, Beijing), and MEGA 7.0 (Kumar et al. 2016) was used to assess polymorphisms at each locus. In total, 47 microsatellite loci with one or two alleles were excluded, and each of the remaining 52 microsatellite loci was labelled with one fluorescent dye (FAM, HEX, or TAMRA) to the forward primers. We then amplified all these loci in 51 samples and performed genotyping using Genemapper 4.0 (Applied Biosystems).

Statistical and genetic data analysis

The number of alleles per locus (N_a), observed heterozygosity (H_o), expected heterozygosity (H_e), and polymorphic information content (PIC) were computed via CERVUS 3.0 (Marshall et al. 1998). Tests for deviations from the Hardy-Weinberg equilibrium (HWE), linkage-disequilibrium (LD), and the inbreeding coefficient (F_IS) were conducted by Arlequin 3.5.2.2 (Excoffier and Lischer 2010). The resulting p values were adjusted by applying a sequential Bonferroni correction using SGoF+ 3.8 (Carvajal-Rodriguez and de Uña-Alvarez 2011). A maximum-likelihood estimate of the frequency of null alleles and allelic dropout was calculated for each locus using Micro-Checker 2.2.3 (Van Oosterhout et al. 2004). Genetic differential index (F_ST) for the populations between two mountains was estimated in GENETIX 4.05.2 using 1000 bootstraps to calculate significance (Belkhir et al. 2004). We used BayesAss 3.03 to estimate rates of recent immigration between populations from two mountains (Wilson and Rannala 2003). Model parameters were set to default values, and the MCMC was run for 4, 000, 000 iterations, with a burn-in of 1, 000, 000 iterations and a sampling frequency of 2000. The chromosome location of each microsatellite was obtained by comparison with the Zebra Finch (Taeniopygia guttata) genome (http://www.ncbi.nlm.nih.gov/genome/seq/BlastGen.cgi?taxid=59729).

Results

Identification of microsatellite loci

A total of 61, 704, 680 NGS reads (7.18 G bp) were generated in this study. After removing the adapter-related, low quality, and duplicated reads, 56, 901, 588 reads remained (92.22%). We obtained 24, 408, 786 clean reads after RAD tag clustering. These reads were assembled into 973, 893 contigs. The average length and the N50 of contigs were approximately 361 and 510 bp, respectively. The average GC content of retained contigs was 40.41% and the average sequencing depth was 17.03.

In total, 11, 564 putative microsatellites contained 3152 di-nucleotides, 6030 tri-nucleotides, 1886 tetra-nucleotides, 395 penta-nucleotides, and 101 hexa-nucleotides. Within the 600 random microsatellite-containing sequences, there were 119 tri-nucleotides, 250 tetra-nucleotides, 136 penta-nucleotides, and 95 hexa-nucleotides.

Based on the criteria that a microsatellite locus must contain more than two alleles, high heterozygote peak height ratios, minimal stuttering and split peaks (Guichoux et al. 2011), 24 microsatellite primers (Table 1) were selected to genotype 51 Emei Shan Liocichla samples.

Table  1.  Characterization of 24 microsatellite markers from the Emei Shan Liocichla (Liocichla omeiensis)

Locus	Size (bp)	Temperature (℃)	Repeat motif	Primer sequence (5ʹ-3ʹ)	Na	Ho	He	PIC	FIS	r	Chromosome location
LO125	143	61	(ATCT)4	F: AAAAGGTATCATGCCTCTAAGTGG R: GCAAAAATGGGCTGGGTTAGA	9	0.8824	0.8288	0.7970	-0.0653	0.0000	3
LO131	136	56	(CAGA)4	F: CCTCAGCACAGAACTGGACCTT R: AAATCGCAGCCTGAATCTCC	4	0.5294	0.5962	0.5320*	0.1130	0.0517	3
LO141	166	58	(AGAT)4	F: CCTCTTAGAAGGCTGCTGTCC R: CCCTGGTCTTTTCATCTATGTAGTT	6	0.8235	0.7690	0.7230	-0.0717	0.0000	3
LO180	204	55	(AAT)3	F: GGAGAGGAAAATCAGAAAAGGC R: CTGGGCATCAGTCAGTGGA	10	0.9020	0.8393	0.8110	-0.0755	0.0000	3
LO183	182	57	(GAT)3	F: CATCAACATCATCCAACCCAT R: AGGCTGGGTCTCAACAATAAAA	5	0.4510	0.5871	0.4930	0.2336	0.1082	3
LO220	168	61	(TTGG)4	F: GATTTGTTTTTCCTTTTGCCTGC R: GTCGATTGTGGAACTGAGTGAAAA	5	0.5686	0.6327	0.5530	0.1022	0.0395	Z
LO244	145	57	(CAG)3	F: AAAAGAGGGAAACTGGAAGACC R: ACAGCAGGCAACTGGTATCACT	6	0.7451	0.6544	0.6100	-0.1401	0.0000	Z
LO271	191	58	(AAGA)4	F: TATTGATATTTCCTCCCTCACTGG R: CAGAATCATTTTCCCCTGGCTA	5	0.7451	0.7177	0.6560	-0.0385	0.0062	Z
LO299	154	58	(ATAG)4	F: TGTAGACGACCATCTTCTGCTGT R: AGATTGCTGGCTCTTTGGTTATG	6	0.4902	0.7191	0.6700**	0.3205	0.1581	3
LO313	134	57	(TAGG)4	F: CTACGTTACCTTTTTTTCTTCTGC R: CCCACAAATAAGCAACAAAACA	8	0.5294	0.8464	0.8180**	0.3769	0.1783	2
LO321	160	58	(TTCTA)5	F: GTCAGCCCTGAAGTCGACAAGT R: GGATTCTTCAGTTGAGGCTGTATTC	6	0.8431	0.7569	0.7090	-0.1152	0.0000	1
LO348	168	61	(GTTG)4	F: GGCACTGGAAAACAGAGGCA R: GCCTGCCCCTGACATCTGTAT	5	0.7647	0.6564	0.6080	-0.1670	0.0000	3
LO403	162	56	(TGA)3	F: GTGCTTTTTAGGCAGATCTCTCA R: AAGAACAGGACCCCAGAAAAA	6	0.5686	0.5844	0.5470	0.0272	0.0053	12
LO437	174	57	(AAATA)5	F: CTCTTCACCCACTTCACATCTCTT R: GAAAAAAATGATGGACAGAAAAGC	10	0.7843	0.7717	0.7440	-0.0165	0.0016	11
LO438	177	57	(TGC)3	F: AGGCTGCCAGGATGTTGCT R: CAGAGAAGAAGAAAGGGAGGTC	5	0.7451	0.6851	0.6230	-0.0885	0.0000	1
LO445	180	59	(TTAG)4	F: AGATGCAGAAGATTTAGGAGTGTGG R: TGCCCCCAGCACCTGACTA	6	0.6667	0.7593	0.7110	0.1230	0.0573	8
LO464	217	55	(GCA)3	F: AAAGGGGTGGAGGAAGATAACT R: TATGCTCATTTTGCTGTTTGC	4	0.7059	0.7173	0.6620	0.0161	0.0023	20
LO514	195	58	(GAT)3	F: TGGGGAGGTCAAGTTGTTCTG R: CAAAAATTATGTGGACCGATGC	8	0.5098	0.7760	0.7350**	0.3453	0.1666	5
LO520	196	58	(TTTG)4	F: CTGGTGGCTACAGAAGGTGCTC R: GCACTGGTTTAGTCTCTCCTTGGTA	5	0.3922	0.4628	0.4280	0.1540	0.0754	NA
LO55	153	61	(TGTT)4	F: TGGTGGGGATTTTTTGGTTG R: CACAGTGCTTCAAGTTGCTCTCTC	5	0.4118	0.6937	0.6390**	0.4088	0.1865	14
LO566	201	56	(CTG)3	F: AAAATCCTGCCCGACTACGA R: GCTTTATCTCCTCAGCATCTCCT	6	0.6863	0.7556	0.7070	0.0926	0.0413	9
LO579	135	56	(AAAC)4	F: ACAAAATGTAATCAGTTGAAGTCCA R: GGCAAAAGCATTCTGAAAAGGA	5	0.6078	0.6492	0.5890	0.0643	0.0182	3
LO61	156	53	(TCTA)4	F: GCATATACTTGTACCATGGAAAAC R: ACCACACAAACTCCAGAACACT	5	0.8431	0.7290	0.6730	-0.1584	0.0000	1A
LO63	186	56	(GGAT)4	F: GCAATGCAACTCCAAAACTCAG R: CCAGGAAGCAGAAGGACAAAGT	6	0.6863	0.7269	0.6730	0.0564	0.0210	1A
Na, number of alleles at each locus; Ho, observed heterozygosity; He, expected heterozygosity; PIC, polymorphic information content; FIS, inbreeding coefficient; r, frequency of null allele Asterisk means significantly different after Bonferroni correction for HWE (where * denotes deviation with Bonferroni correction p < 0.05 and ** denotes deviation with Bonferroni correction p < 0.0025)

 | Show Table

DownLoad: CSV

Genetic characters of the Emei Shan Liocichla samples

The N_a ranged from 4 to 10 with an average of 6.08. H_o and H_e varied from 0.3922 to 0.9020 (mean = 0.6617) and 0.4628 to 0.8464 (mean = 0.7047), respectively. PIC values ranged from 0.4280 to 0.8180 with an average of 0.6546. After the Bonferroni correction, five loci exhibited deviations from HWE (LO131, LO299, LO313, LO514 and LO55), and two pairwise loci exhibited linkage disequilibrium (LO131 and LO55, LO313 and LO55). Therefore, we removed these five potential problematic loci in following analyses. The overall F_IS value of the remaining 19 loci was -0.00976, significantly different from zero (p = 0.0020). Estimated frequencies of null alleles per locus ranged from 0 to 0.1739, and there is no allelic dropout in each locus. The detailed characteristics of all individuals at the 24 loci are summarized in Table 1.

The overall F_ST between populations from Daxiangling and Xiaoliangshan is 0.022 (p < 0.05). The recent migration rate calculated from BayesAss was 0.3128 from populations in Daxiangling to populations in Xiaoliangshan, and 0.0291 vice versa.

Discussion

Our results identified 24 novel microsatellites of the Emei Shan Liocichla using Illumina paired-end RAD-Seq strategy. The higher average depth per locus yielded by RAD-Seq ensured the accuracy of these microsatellites (Andrews et al. 2016).

PIC is the most widely applied formula to measure the information content of molecular markers (Nagy et al. 2012). According to the criteria of Botstein et al. (1980), 21 of these 24 microsatellite loci characterized by RAD sequencing showed a high level of polymorphism (PIC > 0.5). Other researchers like Huang et al. (2015) and Zhang et al. (2015) who have used other NGS methods to identify microsatellites within threatened bird species have observed similar polymorphism levels. Therefore, the use of RAD sequencing strategy to identify microsatellite is economically appealing and effective (Barchi et al. 2011).

The levels of microsatellite genetic diversity in the Emei Shan Liocichla shown in our results are lower than those in the non-threatened Liocichla steerii (McKay et al. 2010) (Table 2). Although the comparison between different microsatellite loci can be problematic, the N_a in L. steerii (9-51) (McKay et al. 2010) are much higher than in Emei Shan Liocichla (4-10, Table 2), which could still imply that the Emei Shan Liocichla might have lost some genetic diversity. Additionally, the frequency of null alleles was low for most loci (r < 0.05), suggesting that the effective population size for the Emei Shan Liocichla might be low (Chapuis and Estoup 2007), consistent with the low number of mature individuals estimated (1500‒7000) (BirdLife International 2017).

Table  2.  Genetic diversity of the Emei Shan Liocichla (Liocichla omeiensis) and Steere's Liocichla (Liocichla steerii)

Threatened species							Non-endangered species
	A	Ho	He	n	FIS		A	Ho	He	n	FIS
DXL Liocichla omeiensis	5.58	0.668	0.695	27	-0.03910		Liocichla steerii a	19.6	0.76	0.79	122	0.026
XLS Liocichla omeiensis	5.54	0.655	0.698	24	-0.00117
Total	6.08	0.662	0.705	51	-0.00976
DXL, Daxiangling Mountain; XLS, Xiaoliangshan Mountain; A, Average number of alleles per locus; Ho, observed heterozygosity; He, expected heterozygosity; n, sample sizeh aData from McKay et al. (2010)

 | Show Table

DownLoad: CSV

After removing HWE-deviated and LD loci, the overall F_IS showed negative values in all samples (Table 2), which means that the populations of the Emei Shan Liocichla do not exhibit inbreeding. It has been proposed that sexual dimorphism, which implies potential mating preference, can be a possible mechanism to avoid inbreeding (Tregenza and Wedell 2002). The males of Emei Shan Liocichla differ from females in having conspicuous orange-red on flight-feather bases and tips (Collar et al. 2016). In contrast, L. steerii show less sexual dimorphism between males and females, and their F_IS valve was 0.026 (Table 2) even in such large populations. Therefore, the low-level inbreeding in the Emei Shan Liocichla may be due to strong sexual selection, an explanation that needs to be further investigated.

The global F_ST between populations from Daxiangling and Xiaoliangshan was low (0.022), but highly significant, which indicate that there is a modest degree of genetic differentiation between populations from these two mountains. However, the recent migration rate analysis indicated there was relatively high contemporary gene flow from Daxiangling populations to Xiaoliangshan populations (0.3128), suggesting the populations are not strictly isolated yet. It is worth noting that the recent migration rate was highly asymmetrical, implying that a source-sink population structure exists in the Emei Shan Liocichla. Worse yet, the natural forests within its habitat have been cleared or replaced by tea plantations, and although China's Grain-for-Green Program has been in effect since 1999, the restored forests between these two mountains are mostly monoculture forest and semi-mixed forest (Hua et al. 2016), which may be unsuitable for bird dispersal (Twedt et al. 2002). Therefore, although the populations exhibit low levels of genetic differentiation and moderate unidirectional gene flow now, still the long-term separation between these populations and unsuitable habitat between the two mountains may increase the risk of genetic fragmentation.

Conclusions

Based on the aforementioned results, 19 neutral microsatellites, with considerably high polymorphisms, were explored in the Emei Shan Liocichla by Illumina paired-end RAD-Seq. This set of de novo nuclear markers, combined with other molecular markers, such as mitochondrial DNA and nuclear introns, would facilitate further studies on population demographics, phylogeography and paternity analysis in vulnerable Emei Shan Liocichla and might be employed across species within the Liocichla genus, such as the critically endangered Bugun Liocichla (Liocichla bugunorum). In our genetic diagnoses, the Emei Shan Liocichla showed decreased genetic diversity, and a modest degree of genetic differentiation between populations, although they showed no inbreeding yet. These results indicate that if we conserve this species properly, before the emergence of strict loss of genetic diversity and inbreeding, we may still be able to prevent the population from becoming fully genetically fragmented.

Authors' contributions

AY designed primers, conducted genetic analyses and led efforts to draft the manuscript. DC analyzed the data, interpreted the outcomes and helped draft the manuscript. PW and YF helped with acquired samples. ZZ and YF conceived the ideas, improved the manuscript and directed the research. All authors read and approved the final manuscript.

Acknowledgements

This work was supported by the National Natural Science Foundation of China (Grant No. 31272330) to YF and the National Key Programme of Research and Development, Ministry of Science and Technology (2016YFC0503200) to ZZ. We would like to thank Mr. Jianzhi Zhang for his help in sample collection, and Mr. Yinong Liu for his help in figure preparation. We appreciate Dr. Rowden for his help with the English corrections. We are grateful to Prof. Rebecca Kimball, Prof. Shou-Hsien Li and Dr. Lu Dong for their comments on the manuscript. We also thank the anonymous reviewers for their comments on the manuscript.

Competing interests

The authors declare that they have no competing interests.

Availability of supporting data

The data set supporting the results of this article is available in the NCBI repository. The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.