Selection of a marker gene to construct a reference library for wetland plants, and the application of metabarcoding to analyze the diet of wintering herbivorous waterbirds

Ethics statement

Our research work did not involve capture or any direct manipulation or disturbances of animals. We collected samples of plants and feces for molecular analyses. We obtained access to the reserve under the permission of the Shengjin Lake National Nature Reserve Administration (Chizhou, Anhui, China), which is responsible for the management of the protected area and wildlife. We were forbidden to capture or disturb geese in the field.

Study area

Shengjin Lake (116°55′–117°15′E, 30°15′–30°30′N) was established as a National Nature Reserve in 1997, aiming to protect waterbirds including geese, cranes and storks. The water level fluctuates greatly in this lake, with maximal water level of 17 m during summer (flood season) but only 10 m during winter (dry season). Due to this fluctuation, receding waters expose two large Carex spp. meadows and provide suitable habitats for waterbirds. This makes Shengjin Lake one of the most important wintering sites for migratory waterbirds (Zhao et al., 2015). Greater white-fronted goose and bean goose are the dominant herbivores wintering (from October to April) in this area, accounting for 40% and 60% of populations along the East Asian–Australian Flyway Route, respectively (Zhao et al., 2015).

Field sampling

The most common plant species that these two geese may consume were collected in May 2014 and January 2015, especially species belonging to Carex and Poaceae. Fresh and intact leaves were carefully picked, tin-packaged in the field and stored at −80 °C in the laboratory before further treatment. Morphological identification was carried out with the assistance of two botanists (Profs. Zhenyu Li and Shuren Zhang from Institute of Botany, Chinese Academy of Sciences).

All feces were collected at the reserve (Fig. 2) in January 2015. Based on previous studies and the latest waterbird survey, sites with large flocks of geese (i.e., more than 200 individuals) were chosen (Zhang et al., 2011). As soon as geese finished feeding and feces were defecated, fresh droppings were picked and stored on dry ice. Droppings of bean geese were generally thicker than those of smaller greater white-fronted goose, to the degree that these could be reliably distinguished in the field (Zhao et al., 2015). Disposal gloves were changed for each sample to avoid cross contamination. To avoid repeated sampling and to make sure samples were from different individuals, each sample was collected with a separation of more than 2 m. In total, 21 feces were collected, including 11 for greater white-fronted goose and 10 for bean goose. All samples were transported to laboratory on dry ice and then stored at −80 °C until further analysis.

Selection of molecular markers and corresponding primers

Here, we aimed to select gene markers with adequate discriminating power for our study. We included eight chloroplast genes— rbcL, rpoC1, rpoB, matK, trnH-psbA, trnL (UAA), atpF-atpH, and psbK-psbI for estimation. Although Shengjin Lake included an array of plant species, we focused mainly on the most likely food resources (Xu et al., 2008; Zhao et al., 2015) that geese would consume for candidate gene tests. These covered 11 genera and the family Poaceae (Table S1). For tests of all candidate genes, we recovered sequences of representative species in the selected groups from GenBank (http://www.ncbi.nlm.nih.gov/nuccore). We calculated inter-specific divergence within every genus or family based on the Kiruma 2-parameter model (K2P) using MEGA version 6 (Tamura et al., 2013). We also constructed molecular trees based on UPGMA using MEGA and characterized the resolution of species by calculating the percentage of species recovered as monophyletic based on phylogenetic trees (Rf). Secondly, primers selected out of eight candidate genes were used to amplify all specimens collected in Shengjin Lake and to check their amplification efficiency and universality. Thirdly, we calculated inter-specific divergence based on sequences that we obtained from last step. Generally, a robust barcode gene is obtained when the minimal inter-specific distance exceeds the maximal intra-specific distance (e.g., existence of barcoding gaps). Finally, to allow the recognition of sequences after high-throughput sequencing, both of the forward and reverse primers of the selected marker gene were tagged specifically for each sample with 8nt nucleotide codes at the 5′end (Parameswaran et al., 2007).

DNA extraction, amplification and sequencing

Two hundred milligrams of leaf was used to extract the total DNA from each plant sample using a modified CTAB protocol (Cota-Sanchez, Remarchuk & Ubayasena, 2006). DNA extraction of feces was carried out using the same protocol with minor modification in incubation time (elongate to 12 h). Each fecal sample was crushed thoroughly and divided into four quarters. All quarters of DNA extracts were then pooled together. DNA extraction was carried out in a clean room used particularly for this study. For each batch of DNA extraction, negative controls (i.e., extraction without feces) were included to monitor possible contamination.

For plant DNA extracts, PCR amplifications were carried out in a volume of 25 µl with ∼100 ng total DNA as template, 1U of Taq Polymerase (Takara, Dalian, Liaoning Province, China), 1× PCR buffer, 2 mM of Mg^{2 +}, 0.25 mM of dNTPs, 0.1 µM of forward primer and 0.1 µM of reverse primer. After 4 min at 94 °C, the PCR cycles were as follows: 35 cycles of 30 s at 94 °C, 30 s at 56 °C and 45 s at 72 °C, and the final extension was 10 min at 72 °C. We applied the same PCR conditions for all primers. All the successful PCR products were sequenced with Genewiz (Suzhou, Jiangsu Province, China).

For fecal DNA extracts, PCR mixtures (25 µl) were prepared in six replicates for each sample to reduce biased amplification. Each replicate was subjected to the same amplification procedure used for plant extracts. The six replicates of each sample were pooled and purified using the Sangon PCR product purification kit (Sangon Biotech, Shanghai, China). Quantification was carried out to ensure equilibrium of contribution of each sample using the NanoDrop ND-2000 UV-Vis Spectrophotometer (NanoDrop Technologies, Wilmington, Delaware, USA). High-throughput sequencing was performed using Illumina MiSeq platform following manufacturer’s instructions by BGI (Shenzhen, Guangdong Province, China). Reads of high-throughput sequencing could be found at NCBI’s Sequence Read Archive (Accession number: SRP070470).

Data analysis for estimating diet composition

After high-throughput sequencing, pair-ended reads were merged with the fastq_mergepairs command using usearch (http://drive5.com/usearch, Edgar, 2010). Reads were then split into independent files according to unique tags using the initial process of RDP pipeline (https://pyro.cme.msu.edu/init/form.spr). We removed sequences (i) that didn’t perfectly match tags and primer sequences; (ii) that contained ambiguous nucleotide (N’s). Tags and primers were then trimmed using the initial process of RDP pipeline. Further quality filtering using usearch discarded sequences with (i) quality score less than 30 (<Q30) and (ii) length shorter than 100 bp. Unique sequences were clustered to operational taxonomy units (OTUs) at the similarity threshold of 98% (Edgar, 2013). All OTUs were assigned to unique taxonomy with local blast 2.2.30+ (Altschul et al., 1990). We detected a plant within the reference library for each sequence with the threshold of length coverage >98%, identity >98% and e-value <1.0e⁻⁵⁰. If a query sequence matched two or more taxa, it was assigned to a higher taxonomic level which included all taxa.

Microhistology analysis

We used the method described by Zhang et al. (2011) to perform microhistologic examination of fecal samples. Each sample was first washed with pure water and filtered with a 25-µm filter. Subsequently, the suspension was examined under a light microscope at 10× magnification for quantification statistics and at 40× magnification for species identification. We compared photos of visible fragments with an epidermis database of plants from Shengjin Lake to identify food items (Zhang et al., 2011).

Selection of genes and corresponding primers and reference library construction

A total of 3,296 representative sequences were recovered from GenBank, ranging from 0 to 345 sequences per gene per taxon (Table S1). For Eleocharis and Trapa, only sequences of rbcL gene and trnL gene were retained, which makes it unfair to compare the efficiency and suitability of eight candidate genes. For the other ten taxa, trnL, trnH-psbA, rbcL and psbK-psbI showed the largest inter-specific divergence in five, three, one, and one taxonomic groups, respectively. In addition, trnH-psbA, atpF-atpH, trnL and psbK-psbI showed the highest mean divergence in four, four, one and one taxonomic groups, respectively. However, given the small number of sequences and coverage of species, the suitability and efficiency of atpF-atpH and psbK-psbI seem to be less reliable than others. This comparison makes trnH-psbA, trnL and rbcL to be selected out of the eight candidate genes. As matK used to be recommended as the standard barcode gene for Carex species (Starr, Naczi & Chouinard, 2009), which happened to be the dominant food for herbivorous geese in our study (Zhao et al., 2015), we included matK as a supplement at last.

Primers for these four genes (Table 1) were used to amplify the plants that we collected in the field. In total, we collected 88 specimens in the field, belonging to 25 families, 53 genera and 70 species (Table 2). The selected primers for trnL and rbcL successfully amplified 100% and 91% of all species, respectively, while primers for trnH-psbA and matK amplified only 71% and 43%, respectively. Therefore, we chose trnL and rbcL to test their discriminating power in our target plants.

We calculated the inter-specific divergence within genera and families with at least two species to compare their discriminating power. Maximal, minimal and mean inter-specific distances were calculated for seven dominant genera and six dominant families (Table 3). Neither gene could differentiate species of Vallisneria (mean =0.000 ± 0.000%) or Artemisia (mean =0.000 ± 0.000%). But trnL showed a larger divergence range for the other six genera and five families. Hence, we chose trnL as the barcoding gene for reference library constructing and high-throughput sequencing for our study. The discriminating power of trnL was strong for most species (Table 4). However, some species could only be identified at genus-level or family-level with trnL. For instance, five species of Potamogetonaceae shared the same sequences and this made them to be identified at genus-level. Species could be identified easily to genus and family, except for three grasses (Poaceae) Beckmannia syzigachne, Phalaris arundinacea, and Polypogon fugax which shared identical sequences.

Data processing for estimating diet composition

In total, 0.21 and 0.18 million reads were generated for greater white-fronted goose (GWFG) and bean goose (BG), respectively (Table 5). The number of recovered OTUs ranged from 8 to 123 for GWFG and BG samples. We used local BLAST to compare these sequences with the Shengjin Lake reference database. Finally, with DNA metabarocoding, 12 items were discovered in the feces of GWFG, including one at family-level, three at genus-level and eight at species-level (Table 6). Four items were discovered in the feces of BG, including one at genus-level and three at species-level. In total, this method identified 15 taxa in feces of these geese.

However, the sequence percentage of each food item varied greatly (Table 6). For GWFG, the majority of sequences (96.36%) were composed of only five items—Poaceae spp. (47.98%, except Poa annua), Poa annua (21.86%), Carex heterolepis (17.51%), Carex spp. (9.01%, except Carex heterolepis), and Alopecurus aequalis (3.21%). For BG, almost all the sequences belonged to Carex heterolepis (99.49%). Other items only occupied a relatively small proportion of sequences. In addition, the presence of each item per sample was also unequal (Table S2). For example in GWFG, Carex heterolepis, Carexspp., Poa annua and Potentilla supina were present in almost all the samples, while Stellaria media, Asteraceae sp. and Lapsana apogonoldes occurred in only about one third of samples.

When the microhistologic examination was performed using the same samples, eight items were found in the feces of greater white-fronted goose, including one at family-level, four at genus-level and three at species-level (Table 6). Dominant items were Poaceae spp. (45.68%), Alopecurus Linn. (30.93%) and Carex heterolepis (16.39%). Seven items were found in the feces of bean goose, including four at genus-level and three at species-level (Table 6). Dominant items were Carex heterolepis (62.85%), Asteraceae sp. (14.55%), and Alopecurus Linn. (13.18%).

Marker selection and reference library constructing for diet analysis

With greatly reduced cost, extremely high throughput and information content, metabarcoding has revolutionized the exploration and quantification of dietary analysis with noninvasive samples containing degraded DNA (Fonseca et al., 2010; Shokralla et al., 2014). Despite enormous potential to boost data acquisition, successful application of this technology relies greatly on the power and efficiency of genetic markers and corresponding primers (Bik et al., 2012; Zhan et al., 2014). In order to select the most appropriate marker gene for our study, we compared the performance of eight commonly used chloroplast genes (rbcL, rpoB, rpoC1, matK, trnL, trnH-psbA, atpF-atpH, and psbK-psbI) and their corresponding primers. Although a higher level of discriminating power was shown in several studies, atpF-atpH and psbK-psbI were not as commonly used as other barcoding genes (Hollingsworth, Graham & Little, 2011). As one of the most rapidly evolving coding genes of plastid genomes, matK was considered as the closest plant analogue to the animal barcode COI (Hilu & Liang, 1997). However, matK was difficult to amplify using available primer sets, with only 43% of successful amplification in this study. In spite of the higher species discrimination success of trnH-psbA than rbcL + matK in some groups, the presence of duplicated loci, microinversions and premature termination of reads by mononucleotide repeats lead to considerable proportion of low-quality sequences and over-estimation of genetic difference when using trnH-psbA (Graham et al., 2000; Whitlock, Hale & Groff, 2010). In contrast, the barcode region of rbcL is easy to amplify, sequence, and align in most plants and was recommended as the standard barcode for land plants (Chase et al., 2007). The relatively modest discriminating power (compared to trnL) precludes its application for our study aiming to recover high resolution of food items. Consequently, trnL was selected out of eight candidate markers, with 100% amplification success, more than 90% of high quality sequences, and relatively large inter-specific divergence.

One of the biggest obstacles in biodiversity assessment and dietary analysis is the lack of a comprehensive reference library, without which it is impossible to accurately interpret and assign sequences generated from high-throughput sequencing (Valentini, Pompanon & Taberlet, 2009; Barco et al., 2015). In this study, we constructed a local reference library by amplifying the most common species (70 morpho-species in total) during the wintering period with the trnL gene. Although not all of them could be identified at species-level with trnL due to relatively low inter-specific divergence, many species could be separated with distinctive sequences. Previous studies have recommended group-specific barcodes to differentiate closely related plants at the species level (Li et al., 2015). For instance, matK has been proved to be more efficient for the discrimination of Carex spp. (Starr, Naczi & Chouinard, 2009). However, the primer set of matK failed to amplify species of Carex spp. in our study, suggesting the universality of selected primer pairs should be tested in each study (Zhan et al., 2014).

Applications of metabarcoding for geese diet analysis

A variety of recent studies have demonstrated the great potential of metabarcoding for dietary analysis, mainly owing to the high throughput, high discriminating power, and the ability to process large-scale samples simultaneously (Creer et al., 2010; Taberlet et al., 2012; Shehzad et al., 2012). In this study, we applied this method to recover diets of herbivorous geese and provided standard protocols for dietary analysis of these two ecologically important waterbirds. Our results further proved the more objective, less experience-dependent and more time-efficient character of DNA metabarcoding. However, not all the species in the reference library could be identified at species-level, owing to low inter-specific divergence. We suggest that multiple group-specific markers to be incorporated in the future, as in De Barba et al. (2014). Two species, Carex thunbergii and Fabaceae sp., were only discovered via microhistologic analysis rather than metabarcoding. This failure might reflect the biased fragment amplification of current technology, of which dominant templates could act as inhibitors of less dominant species (Piñol et al., 2015). However, three species of Poaceae were only discovered using metabarcoding. In total, more taxa and higher resolution were attained using metabarcoding. But microhistology still proved a powerful supplementary. Previous studies using metabarcoding usually detected dozens of food items, even as many as more than one hundred species. For instance, 18 taxa prey were identified for leopard cat (Prionailurus bengalensis) (Shehzad et al., 2012); 44 plant taxa were recovered in feces of red-headed wood pigeon (Columba janthina nitens) (Ando et al., 2013); while more than 100 taxa were found in diet studies of brown bear (Ursus arctos) (De Barba et al., 2014). The relatively narrow diet spectrum of herbivorous geese may lead to misunderstanding that this result of our study is merely an artefact due to small sampling effort. However, this result is credible since these two geese species only feed on Carex meadow, where the dominant vegetation is Carex spp., with other species such as Poaceae and dicots (Zhao et al., 2015). Even though other wetland plants exist, they usually composed only a small proportion of the geese diets.

Quantification of food composition is another key concern in dietary analysis. Although the relative percentage of sequences was not truly a quantitative estimate of diet, taxa of the majority sequences in this study were in accord with microhistologic observations, which was considered an efficient way to provide quantitative results (Wang et al., 2013). Discrepancies might come from the semi-quantitative nature of metabarcoding methods (Sun et al., 2015). This is likely derived from PCR amplification, which always entails biases caused by universal primer-template mismatches, annealing temperature or number of PCR cycles (Zhan et al., 2014; Piñol et al., 2015). Other methods such as shotgun sequencing or metagenomic sequencing could be incorporated in the future to give information on abundances of food items (Srivathsan et al., 2015).

Implications for waterbird conservation and wetland management

For long-distance migratory waterbirds, such as the wild geese in this study, their abundance and distribution are greatly influenced by diet availability and habitat use (Wang et al., 2013). For example, waterbirds may be restricted at (forced to leave) certain areas due to favoring (loss) of particular food (Wang et al., 2013), while the recovery of such food may contribute to return of bird populations (Noordhuis, Van der Molen & Van den Berg, 2002). Results of both metabarcoding and microhistologic analysis in this study revealed that Carex and Poaceae were dominant food components which is in accordance with previous studies. The increasing number of these two geese wintering at the Shengjin Lake may be attributed to the expansion of Carex meadow, which offers access to abundant food resources (Zhao et al., 2015). Considering the long-distance migratory character of these birds, it is important to maintain energy balances and good body conditions in wintering areas because this might further influence their departure dates and reproductive success after arriving at breeding areas (Prop, Black & Shimmings, 2003). Based on this, it is important for wetland managers to maintain the suitable habitats and food resources for sustainable conservation of waterbirds, which highlights the significance of diet information. Our study also indicated that overlap and dissimilarity existed between the diets of these two geese. Animals foraging in the same habitats may compete for limited food resources (Madsen & Mortensen, 1987). This discrepancy of food composition may arise from the avoidance of inter-specific competition (Zhao et al., 2015). However, with the increase of these two species in Shengjin Lake, further research is needed to investigate the mechanisms of food resource partitioning and spatial distribution.

Shengjin Lake is one of the most important wintering sites for tens of thousands of migratory watebirds, while annual life cycles of these birds depend on the whole migratory route, including breeding sites, stop-over sites and wintering sites (Kear, 2005). Thus, a molecular reference library covering all the potential food items along the whole migratory route will be useful both for understanding of wetland connections and waterbird conservation. In addition, the ability of DNA metabarcoding to process lots of samples simultaneously enables rapid analyses and makes this method helpful for waterbird studies.