Characterization of the complete mitochondrial genome and a set of polymorphic microsatellite markers through next-generation sequencing for the brown brocket deer Mazama gouazoubira.

The complete mitochondrial genome of the brown brocket deer Mazama gouazoubira and a set of polymorphic microsatellite markers were identified by 454-pyrosequencing. De novo genome assembly recovered 98% of the mitochondrial genome with a mean coverage of 9-fold. The mitogenome consisted of 16,356 base pairs that included 13 protein-coding genes, two ribosomal subunit genes, 22 transfer RNAs and the control region, as found in other deer. The genetic divergence between the mitogenome described here and a previously published report was ∼0.5%, with the control region and ND5 gene showing the highest intraspecific variation. Seven polymorphic loci were characterized using 15 unrelated individuals; there was moderate genetic variation across most loci (mean of 5.6 alleles/locus, mean expected heterozygosity = 0.70), with only one locus deviating significantly from Hardy-Weinberg equilibrium, probably because of null alleles. Marker independence was confirmed with tests for linkage disequilibrium. The genetic variation of the mitogenome and characterization of microsatellite markers will provide useful tools for assessing the phylogeography and population genetic patterns in M. gouazoubira, particularly in the context of habitat fragmentation in South America.

The brown brocket deer (Mazama gouazoubira), one of the ten recognized species of brocket deer, has the greatest ecological plasticity among these deer and occurs in all biomes of Brazil (except the Amazon basin), the pre-Andean regions in Argentina and Bolivia, southern Peru, Paraguay and Uruguay (Black-Décima ). The relatively wide distribution of this species, its occurrence in a number of protected areas and its apparently large populations that do not appear to be declining at a rate sufficient to pose a threat to survival at this time mean that the M. gouazoubira is classified as being of "least concern" by the IUCN (2012). Although the brown brocket deer is not on the Brazilian national list of threatened species, it has nevertheless been classified as "vulnerable" in the list of threatened species in the southern Brazilian state of Rio Grande do Sul (Marques, 2002) and as "endangered" in the state of Rio de Janeiro (SEMA, 1998). The main threats to this species are poaching and habitat loss by the expansion of agriculture, which has destroyed remnants of forests and large portions of the Cerrado (Brazilian savanah); the latter has led to a reduction and fragmentation of deer populations (Black-Décima ). Despite these threats, the impact of poaching and deforestation/fragmentation in South America on the genetic variability of brown brocket deer populations has seldom been investigated, partly because of our limited knowledge of the biology, ecology and population genetic structure of most Mazama species. Population genetic data in particular can provide valuable information for monitoring species that are of interest for management and conservation (Schwartz ). Genetic markers such as microsatellites allow the population genetic structure to be easily inferred (Sunnucks, 2000). These markers can be developed quickly and at low cost for most species using next-generation sequencing (NGS) (Abdelkrim ; Castoe ; Gardner ; Ekblom and Galindo, 2011).

Six morphological forms of gray brockets have been considered as subspecies of M. gouazoubira (Rossi RV, 2000, MSc dissertation, Universidade de São Paulo, São Paulo, Brazil). However, no genetic study has assessed the validity of the current taxonomy nor has there been any attempt to understand the evolutionary relationships among these taxa. Recent studies have demonstrated that NGS can be used to effectively recover additional molecular markers, including whole mitochondrial genome (mitogenome) sequences (Seabury ; Miller ). Investigation of genetic variation at mitochondrial loci can be particularly useful for assessing patterns of sex-biased dispersal, species phylogeography and demographic histories (Avise, 2000). The number of reports of complete mitogenomes has increased steadily since the report of the first complete mitochondrial genomes based on work using the polymerase chain reaction (PCR) (Arnason ); however, only a few of these have dealt with Neotropical deer.

In this study, we developed a suite of molecular markers for the brown brocket deer M. gouazoubira, the most abundant deer species of the Neotropical region (Duarte, 1996). In order to characterize the genetic structure of M. gouazoubira and determine the impact of habitat fragmentation on the genetic variability of deer populations, we undertook a partial 454-pyrosequencing run based on which we designed a set of seven polymorphic microsatellite markers and recovered the full mitochondrial DNA sequence. Combined, these markers will provide a valuable resource for investigating the population genetics and demographic history of M. gouazoubira in response to habitat fragmentation in South America.

DNA sequencing was done with 454-pyrosequencing on a Genome Sequencer FLX-454 System (GS FLX sequencer). Sample preparation and DNA sequencing were done according to the manufacturer’s instructions (Roche Diagnostics, Mannheim, Germany). Approximately 100 μg of genomic DNA was extracted from a blood sample of a single female (from the Brazilian Pantanal) using digestion with proteinase K followed by purification with a standard phenol-chloroform-isoamyl alcohol procedure (Sambrook ). Sequencing was done using half of a 70 × 75 mm Pico Titer Plate in conjunction with a Roche GS FLX (454) system at the Genomic Sciences and Molecular Biotechnology Laboratory at the Catholic University of Brasília, Brasília, DF, Brazil.

Unique sequence contigs possessing microsatellite motifs were identified using MSATCOMMANDER (Faircloth, 2008). Subsequently, each unique contig with enough flanking sequence was selected for primer design using PRIMER3 (Rozen and Skaletsky, 1999). A selection of contigs including di-, tri-, tetra-, penta- and hexa-nucleotide repeats was used for subsequent analysis. The forward primer for each locus was appended with a universal M13 sequence (5′CACGACGTTGTAAAACGAC) at the 5′ end. To enable fluorescent labeling, the complement sequence to the M13 tail was modified with a 5′ fluorophore (FAM or HEX). Amplification reactions were done with 45 ng of M. gouazoubira DNA in 15 μL reaction volumes containing 50 mM KCl, 10 mM Tris-HCl (pH 8.4), 2 mM MgCl2, 0.12 mM each dNTP, 0.02 μM forward primer, 0.08 μM fluorescent labelled M13 primer, 0.08 μM reverse primer and 0.4 U of Taq DNA polymerase. Amplification was done in a Biometra T-1 Thermoblock with the following conditions: 95 °C for 5 min, 35 cycles at 94 °C for 1 min, 50–59 °C for 1 min, 72 °C for 1 min, and a final extension step at 72 °C for 10 min. The PCR products were run on 2% agarose gels to confirm the amplification and check the fragment size; 1 kb plus DNA ladder markers were included in all runs (Invitrogen). Genotyping was subsequently done using an Applied Biosystems 3130 capillary analyzer and product lengths were scored manually and assessed for polymorphisms using GENEMARKER version 2.4.2 (Softgenetics).

Loci were pooled into groups for multiplexing based on observed locus-specific allele size ranges and screened for polymorphisms using template DNA from 15 individuals captured in the wild, in the central region of the Brazilian Pantanal (state of Mato Grosso do Sul) known as Nhecolândia. Individuals were captured and sampled in July, 2011 by the field staff of the Center for Research and Conservation of Cervidae (NUPECCE), at the College of Agricultural and Veterinary Sciences at Jaboticabal (FCAVJ/UNESP). IDENTITY version 1.0 software (Wagner and Sefc, 1999) was then used to estimate the expected (HE) and observed (HO) heterozygosities, number of alleles (NA), paternity exclusion probability (Q) and probability of identity (I), whereas conformation to Hardy-Weinberg equilibrium and linkage disequilibrium estimates between all pairs of loci were examined using open-source GENEPOP (Web software version 4; Raymond and Rousset, 1995). The significance values were adjusted for multiple comparisons using Bonferroni corrections (Rice, 1989). All loci were assessed using MICRO-CHECKER to check for null alleles and scoring errors (Van Oosterhout ).

Sequence reads in SFF format were edited by trimming 454 adaptor sequences and converted to fasta and quality file format using the 454 Sequencing system software. Genomic sequence contigs were assembled using de novo default 454 parameters and coverage statistics were obtained with GENEIOUS version 5.6.6 (Drummond ). Gene positions, codon usage, transcriptional orientations and secondary structures were obtained with open-source DOGMA software (Wyman ), using 60% and 80% identity cutoffs for protein coding genes and RNAs, respectively. Annotations and reading frames were confirmed by visual inspection in GENEIOUS using the previously published mitogenome of M. gouazoubira (Genbank accession number JN632658) as a reference. The ribosomal subunit (rRNA) gene boundaries were estimated with alignments implemented in GENEIOUS using the M. gouazoubira genomes, with a high degree of conservatism at the beginning and end of the respective genes across the mitogenome. GENEIOUS was also used to generate the finalized annotated mitochondrial genome map (Genebank accession number KJ772514). Gene arrangement, composition and transcriptional polarity were compared with the mitogenome of other Neotropical deer, including Hippocamelus antisensis (GenBank accession number JN632646), Ozotoceros bezoarticus (JN632681), Blastocerus dichotomus (JN632603) and Mazama nemorivaga (JN632660). These species were chosen due to close phylogenetic relationship with M. goauzoubira (Hassanin ). Intraspecific variation in the complete mitochondrial genome of M. gouazoubira was quantified by comparing the mitogenome described here with that previously published for M. gouazoubira from Colombia (GenBank accession number JN632658).

A total of 576,646 reads with an average length of 430 bp (after trimming 454 adapter sequences) and covering up to 247 Mb of the M. gouazoubira genome was obtained by NGS. These data likely represent ∼7.4% of the genome, based on previous deer genome size estimates, e.g., ∼3,364 Mb for Muntiacus muntjak (Kent ). However, these figures are probably overestimated because of expected read redundancy.

Only singleton reads (364,782) were considered in the search for microsatellite loci. A total of 18,560 reads possessing microsatellite motifs were identified by MSATCOMMANDER analysis, of which 8,856 were found to possess optimal priming sites. Twenty loci were screened for polymorphism, with nine containing tri-nucleotide repeats and 11 containing tetranucleotide repeat motifs. The screening analysis showed that seven loci were polymorphic, three were monomorphic, five yielded alleles with a standard impossible to interpret and five failed to amplify.

All polymorphic loci were characterized by low to moderate genetic variation, with an average of 5.6 alleles per locus (range: 3–8 alleles) and heterozygosity estimates ranging from 0.27 to 1.00 (mean: 0.55; Table 1). Linkage disequilibrium analyses confirmed marker independence since none of the 21 pairwise tests was significant after Bonferroni correction. All loci conformed to Hardy-Weinberg expectations, except Mgoua7, which showed an excess of homozygotes (Table 1). MICRO-CHECKER detected the presence of null alleles only at this locus.

Table 1

Characterization of 10 microsatellite loci developed for the brown brocket deer (Mazama gouazoubira) by NGS technology.

Locus (GenBank number)	Repeat motif	Primer pair sequence (5′-3′)	Ta (°C)	N	NA	Size-range (bp)	HE	HO	Q	I
Mgoua7 (KJ772506)	(CAT)15	F - GTCTGATCCAAATCTGAGGGTCR - TCATAGCATCCAAGGCAAACTA	55	11	4	159–168	0.66	0.27*	0.406	0.169
Mgoua10 (KJ777821)	(TTG)13	F - TAGTGGGACGTTTGTTGTTGTTR - TGGATCTTTGGAGAGGGTCTAA	55	11	1	134	-	-	-	-
Mgoua16 (KJ772507)	(TACA)12	F - GGGACAGTGATAAACTAGGTGTR - CTAATGAGATAGCAAAGTACGC	55	9	8	209–241	0.84	1.00	0.690	0.043
Mgoua17 (KJ772508)	(ACA)11	F - GCAATCCCTATCAAGCTACCAAR - TGATTCCTCCAGTTCCATTCTT	55	11	1	175	-	-	-	-
Mgoua18 (KJ777822)	(ATT)11	F - TTCCAGGCAAGAATACAGGAGTR - GTAACTCGTTGAGCATAAGGGC	52	12	1	205	-	-	-	-
Mgoua19 (KJ772509)	(TCTA)11	F - GGCGATAGCTTGGTGTCATTAR - GCAGAGCATCAGTAGAGGTTCA	50	12	5	154–190	0.34	0.38	0.192	0.453
Mgoua20 (KJ772510)	(ATAA)11	F - ACAACTGGAGAAAACCCTTGTGR - AGCCTTTAGAGATGTTCTGTTTGG	55	14	3	238–274	0.42	0.54	0.214	0.380
Mgoua21 (KJ772511)	(CATA)11	F - GAGTACAACAGCCATGCAGAGAR - CATTGGGGTTCACCTAGAGAAG	52	15	7	164–188	0.74	0.60	0.536	0.099
Mgoua24 (KJ772512)	(ACA)10	F - AAGAAGCTCAAACTTGCCTGTCR - TCTTATTTCCACCTCTTTCCCA	55	12	6	171–186	0.80	0.83	0.610	0.068
Mgoua25 (KJ772513)	(CAT)10	F - AGGACAACCATGCACCTACTTTR - ATCCCAGCTCCTTTTAACACAA	55	12	6	191–209	0.71	0.62	0.472	0.132
All loci				12	5.6		0.70	0.55	0.989	1.083×10−6

H - expected heterozygosity, H - observed heterozygosity, I - probability of genetic identity, N - number of individual analyzed, NA - number of alleles, Q - paternity exclusion probability, and Ta - annealing temperature.

Loci with Hardy-Weinberg disequilibrium.

Approximately 0.06% of the total NGS reads (358 reads) were of mitochondrial origin, and de novo assembly of mtDNA sequence contigs revealed about 98% of mitogenome coverage with a mean coverage of nine-fold. To complete the mitogenome sequence, two fragments were amplified (540 bp and 322 bp) via PCR with primers designed based on the previously obtained genome and later sequenced by the Sanger method on an ABI 3130 sequencer. The absence of internal stop codons, ambiguous base calls and evidence of heteroplasmy suggested that our sequences were authentic mitochondrial targets rather than nuclear mitochondrial-like sequences (numts) which are common in mammals (Bensasson ). The nucleotide composition of the a-strand was 33.9% adenine, 22.6% cytosine, 12.9% guanine and 30.6% thymine.

The M. gouazoubira mitochondrial genome is a 16,356 bp circular molecule with a typical metazoan-gene composition that includes 13 protein-coding genes (ND1, ND2, COI, COII, ATP8, ATP6, COIII, ND3, ND4L, ND4, ND5, ND6 and Cytb), two ribosomal subunits RNA (srRNA and lrRNA), 22 transfer RNAs (trnX) and the control region (Figure 1; Table 2). The overall gene arrangement, including respective transcriptional polarities of genes, was similar to that previously described for M. gouazoubira and the other deer that we have compared, and typical of mammals in general (Gibson ). The main difference between the genome described here and those of other deer species was in relation to size, i.e., 16,357 bp in O. bezoarticus, 16,359 bp in B.dichotomus and M. nemorivaga, and 16,410 bp in H. antisensis. The larger size of the latter resulted from a 45 bp insertion (from 16,160 to 16,205) in the terminal portion of the control region.

Figure 1

Map of the brown brocket deer mitochondrial genome. ATP6-8: ATPase subunits 6 and 8, COX1-3: cytochrome c oxidase subunits I–III, CYTB: cytochrome b, ND1-6/4L: NADH dehydrogenase subunits 1-6/4L. All 22 transfer RNA genes are designated by tRNA-X. Arrowheads indicate the direction of transcription.

Table 2

Characteristics of the brown brocket deer (Mazama gouazoubira) mitochondrial genome.

Gene	Position		Size (bp)	Strand	Codon
	From	To			Initial	Terminal
trnFPhe (GAA)	1	69	69	H
srRNA (12S)	70	1023	954	H
trnVVal (UAC)	1024	1090	67	H
lrRNA (16S)	1091	2657	1567	H
trnLLeu (UAA)	2658	2732	75	H
NADH dehydrogenase subunit 1 (ND1)	2735	3691	957	H	ATG	TAA
trnIIle (GAU)	3691	3759	69	H
trnQGln (UUG)	3757	3828	72	L
trnMMet (CAU)	3831	3899	69	H
NADH dehydrogenase subunit 2 (ND2)	3900	4940	1041	H	ATA	TAG
trnWTrp (UGA)	4942	5008	67	H
trnAAla (UGC)	5010	5078	69	L
trnNAsn (GUU)	5080	5152	73	L
trnCCys (GCA)	5185	5250	66	L
trnYTyr (GUA)	5251	5318	68	L
Cytochrome c oxidase subunit 1 (COI)	5320	6864	1545	H	ATG	TAA
trnSSer (UCN)	6862	6930	69	L
trnDAsp (GUC)	6938	7005	68	H
Cytochrome c oxidase subunit 2 (COII)	7007	7690	684	H	ATG	TAA
trnKLys (AAA)	7694	7762	69	H
ATPase subunit 8 (ATPase8)	7764	7964	201	H	ATG	TAA
ATPase subunit 6 (ATPase6)	7925	8605	681	H	ATG	TAA
Cytochrome c oxidase subunit 3 (COIII)	8605	9388	784	H	ATG	TA-
trnGGly (UCC)	9389	9457	69	H
NADH dehydrogenase subunit3(ND3)	9458	9804	347	H	ATA	TA-
trnRArg (UCG)	9805	9873	69	H
NADH dehydrogenase subunit 4L (ND4L)	9874	10170	297	H	GTG	TAA
NADH dehydrogenase subunit 4 (ND4)	10164	11541	1378	H	ATG	T-
trnHHis (GUG)	11542	11610	69	H
trnSSer (AGY)	11611	11670	60	H
trnLLeu (CUN)	11672	11741	70	H
NADH dehydrogenase subunit 5 (ND5)	11742	13562	1821	H	ATA	TAA
NADH dehydrogenase subunit 6 (ND6)	13549	14073	525	L	ATG	TAA
trnEGlu (UUC)	14074	14142	69	L
Cytochrome b (cytb)	14147	15286	1140	H	ATG	T-
trnTThr (UGU)	15290	15359	70	H
trnPPr° (UGG)	15360	15424	65	L
Control region (D-loop)	15425	16356	932	H

The H-strand encoded 12 protein-coding genes, while the L-strand encoded only one gene (ND6; Table 2, Figure 1). Nucleotide overlap, a common feature in deer mitochondrial genomes (Wada ), was observed between ATP8 and 6 (39 nucleotides; Table 2), ATP6 and COIII (one nucleotide), ND4L and ND4 (six nucleotides), and between ND5 and ND6 (17 nucleotides). The standard methionine (ATN) initiation codon was inferred for 12 of the 13 genes while the ND4L gene appeared to use GTG codon, as in the other deer that were compared (except for B. dichotomus) and also Sika deer (Cervus nippon; Wada ) (Table 2). Open reading frames terminated with the typical TAA and TAG codons for nine genes (Table 2). The termination codon for COIII, ND3, ND4 and Cytb ended with poly-adenylation (TA- or T-), as found in other Neotropical deer and also in the genus Cervus(Wada ).

Both ribosomal subunit genes were encoded by the H-strand, with the srRNA (12S) separated from the lrRNA (16S) by trnV in M. gouazoubira and in the other Neotropical deer examined. The genomic position and transcriptional polarity of the rRNA genes were typical of the "commom" vertebrate arrangement (Boore, 1999).

Twenty-two trns corresponding to the standard metazoan gene set were identified on the basis of their respective anticodons and secondary structures (Table 2). Gene lengths were largely congruent with other deer species. All trns could be folded into the canonical cloverleaf structure except for trnS Ser (ucu), which had four nucleotides in the un-paired loop but lacked the DHU arm. This feature is common among deer and metazoans in general (Wolstenholme, 1992).

A total of 992 noncoding nucleotides were identified, with 60 bp spread across 14 intergenic regions and a large contiguous 932 bp noncoding region (Table 2). The large noncoding region probably represented the putative control region (D-loop) based on its relative position between the trnP° and trnF (Boore, 1999).

A high similarity of 99.5% was found between the complete mitogenomic sequences of M. gouazoubira from Brazil and Colombia. Sequence similarity of the two mitogenomes was 99.5% in the other gene regions except for the control region (D-loop). A sequence similarity of 99.4% was found in the 13 protein-coding gene regions, whereas 99.8% was identical in the mitogenome regions covering all the 22 trns genes. There was 0.21 and 0.32 variation in the 12S and 16S rRNAs, respectively. The control region and the ND5gene were the regions with the highest percentage of variation, with the exception of some trns genes (Figure 2).

Figure 2

Percentage of variation per locus obtained by comparison of the complete mitogenome sequences of M. gouazoubira described here with that previously described from Colombia (GenBank acession number JN632658).

The 454 NGS platform has become a commonly used tool for the development of genetic markers for systematic research (Mardis, 2008). In this study, we successfully isolated 8,856 microsatellite-containing contigs for M. gouazoubira from a total of 576,646 reads covering ∼7.4% of the genome. From these contigs, seven polymorphic microsatellite markers were successfully characterized. We also used NGS in combination with bioinformatic tools to assemble and annotate the complete mitochondrial DNA sequence of M. gouazoubira. The low combined probability of genetic identity and the high power of paternity exclusion indicated that the battery of microsatellite developed for M. gouazoubira will allow detailed studies of parentage and population differentiation in this species; these topics have rarely been addressed using high-resolution microsatellite markers in Neotropical deer. Furthermore, combining these markers with mitochondrial DNA sequence analysis of the hypervariable regions, e.g., control region and ND5, will provide a valuable resource for investigating the demographic history of M. gouazoubira in response to habitat fragmentation, as well as assist in resolving taxonomic uncertainties in this taxon.