1. Introduction

Human mitochondrial DNA (mtDNA) has proven to be a useful tool for a variety of anthropological investigations such as forensics genetics, human evolutionary history, migration patterns, and population studies ^[1,2,3,4,5]. Sequence diversity within the mitochondrial D-loop hypervariable regions (HVR1 and HVR2) has been applied for this purpose since the level of polymorphism in these regions is high enough to permit its use as an important tool in population diversity studies ^[6,7]. However, most of these studies are based on an analysis of a controlled cohort of individuals which are randomly selected to be representative of the population of the geographical region of interest ^[8,9]. In this study, we present an alternate approach for the analysis of population diversity by targeting human mtDNA directly from environmental waters impacted by human contamination.

DNA is naturally shed into the environment by virtually all animal species through feces, urine, exudates, or tissue residues ^[10,11]. There are numerous sources of human mtDNA in environmental waters. These include fecal waste from combined sewer overflows (CSO), sanitary sewer overflows, household sewage treatment systems, and agriculture/urban runoff ^[12,13]. Human fecal waste has been shown to have large amounts of exfoliated epithelial cells, each cell harboring thousands of mitochondrial copies making mtDNA an adequate molecular target in environmental studies. Recently, several studies have taken advantage of human-specific mtDNA signature sequences to implicate human feces as the primary source of contamination in fecally-contaminated effluents ^[14,15,16]. Consequently, human mtDNA sequences obtained from environmental waters are reliable, quantitative and real-time indicators of diversity of the contributing populations.

With the exception of one study, the aforementioned studies have focused on the detection of mtDNA using qPCR assays to detect fecal pollution sources. Recently, Kapoor et al. ^[13] demonstrated the use of mtDNA sequence analysis to both determine the importance of specific human fecal pollution sources in an urban watershed (Cincinnati, OH), as well as the relative abundance of population haplogroups associated with the contributing populations. We hypothesize that human mtDNA sequences in sewage are reliable, quantitative, and real-time indicators of population diversity in a community. To describe, characterize, and track human population diversity in a watershed region, we used high-throughput DNA sequencing technology to profile the HVR2 sequences in water samples taken from a tropical watershed (Río Grande de Arecibo (RGA), Puerto Rico) impacted by human sewage. Like previous controlled population studies ^[6,9], the single-nucleotide polymorphisms (SNPs) present in HVR2 was used to differentiate populations on the basis of their frequencies of occurrence. Furthermore, we extracted haplotypes and assigned mitochondrial haplogroups to identify the mtDNA biological ancestry of the populations impacting the watershed. We demonstrate the potential of these data for surveying the distribution of population diversity in this region and their intersection with orthogonal data like U.S. Census data. These data establish a regional-scale, baseline population profile, which represents a unique metagenomics tool for studying population diversity, regional migration, and other anthropological investigations.

---

2. Materials and Methods

---

2.1. Study area and sampling sites

The Río Grande de Arecibo (RGA) watershed is located along the western-central part of Puerto Rico and has a catchment area of approximately 769 km², with water flowing northward from the central mountain range into a coastal valley before discharging into the Atlantic Ocean. Multiple point sources, including leaking septic and sewer systems and discharge from wastewater treatment plants (WWTPs) contribute to human fecal pollution in the watershed, in addition to several nonpoint sources associated with recreational activities. Three secondary sewage treatment plants discharge disinfected secondary effluents into the watershed: two drain into Río Cidra and Río Caunillas, tributaries of the RGA, while the third drains directly into the RGA. The water quality of the RGA watershed is a major concern as it is an important drinking water reservoir and some sections are used in recreational activities. Thus, fecal contamination of the RGA is a significant public health concern and has a negative economic impact. Most of the population in the RGA watershed is located in the coastal alluvial plain near the municipality of Arecibo ^[17]. The upper watershed is mostly forested, undeveloped land.

The sampling sites (Figure 1) were identified and assessed for the presence of human fecal contamination through PCR-based detection of human fecal markers as described in a previous study ^[18]. These sites had a high human density based on previously recorded fecal pollution levels and potential impact from human fecal pollution via sewage overflow and watershed runoff ^[18]. Three sites (4, 7, and 10) were located downstream of a wastewater treatment plant (WWTP) for the municipalities of Adjuntas, Utuado, and Jayuya, respectively. Sites 6 and 7 represent sites before and after a WWTP. Site 6 is located approximately 1.62 km upstream from site 7, and site 7 is located 120 m downstream from the sewage treatment plant. Site 8 is located at the mouth of the watershed right before the RGA drains into the Atlantic Ocean and close to the town center of Arecibo.

---

2.2. Sample collection and DNA extraction

Ten samples (Table 1) were chosen from the water samples collected within the RGA watershed sites. The water samples collected within the RGA watershed represented different degrees of human contamination. Water sample collection and DNA extraction was performed as described earlier ^[18,19]. Briefly, all samples were collected using sterile bottles and transported on ice to the laboratory at the University of Puerto Rico—Río Piedras Campus where the samples (100 mL) were filtered through polycarbonate membranes (0.4-µm pore size, 47-mm diameter; GE Water and Process Technologies, Trevose, PA) and stored at −80 ℃ until DNA extraction. The membranes were shipped overnight on dry ice to the EPA laboratory (Cincinnati, OH) for DNA extraction. Total DNA was extracted from filters samples using the PowerSoil DNA isolation kit, following the manufacturer's instructions (Mo Bio Laboratories, Inc.). DNA extracts were stored at −20 ℃ until further processing.

---

2.3. High throughput sequencing

The human mitochondrial hypervariable region Ⅱ sequences were elucidated via Illumina sequencing of HVR2 libraries generated with DNA extracts and barcoded primers HVR2-F (5′-GGTCTATCACCCTATTAACCAC-3′) and HVR2-R (5′-CTGTTAAAAGTGCATACCGCC-3′) ^[13]. We generated PCR amplicon libraries for each water DNA extracts. PCR reactions were performed in 25 μL volumes using the Ex Taq kit (Takara) with 200 nM each of the forward and reverse primer and 2 μL of template DNA. Cycling conditions involved an initial 5 min denaturing step at 94 ℃, followed by 35 cycles of 45 s at 94 ℃, 60 s at 56 ℃, and 90 s at 72 ℃ and a final elongation step of 10 min at 72 ℃. Prior to multiplexed sequencing, PCR products were visualized on an agarose gel to confirm product sizes. Sequencing of the pooled library was performed on an Illumina Miseq benchtop sequencer using pair-end 250 bp kits at the Cincinnati Children's Hospital DNA Core facility. The HVRII sequence of the operator was also determined through Sanger sequencing and confirmed that it did not contribute to experimental data.

---

2.4. Bioinformatics analyses

All HVR2 sequences were sorted according to barcodes and grouped under their respective sampling event. The sequences were processed and cleaned using the software MOTHUR v1.25.1 ^[20]. Briefly, fastq files for forward and reverse reads were used to form contigs which were first screened for sequence length (no greater than 420 bp). To compensate for potential sequencing errors, sequences having an average quality under 20, having ambiguous bases (Ns), or being shorter than 300 bp were discarded. The quality-filtered sequences were then aligned to the revised Cambridge Reference Sequence (rCRS) ^[21] for human mitochondrial DNA (NC_012920.1| Homo sapiens mitochondrion, complete genome); and analyzed by using custom scripts to detect the SNPs present in the sequences. All SNPs with frequency greater than 5% were used for further analyses. Additionally, the sequences were exported to CLC Genomics Workbench Version 6.5 (CLC Bio, Cambridge, MA) and aligned to the rCRS, after which the Quality-based Variant Detection was called to detect insertions and deletions (indels) as well as SNPs with reference to the rCRS as described previously ^[13]. The mitochondrial genome databases, including MITOMAP ^[22], mtDB ^[23] and Phylotree ^[24] were referred to validate the occurrence of detected variants. Haplotypes were extracted and submitted to MITOMASTER version Beta 1 ^[25] to assign mitochondrial haplogroups based on variants present in HVR2.

---

3. Results and Discussion

---

3.1. Variant detection and frequency

In total, more than 100, 000 sequence reads were retrieved with a mean output exceeding 20, 000 per barcoded sample, which were then filtered and grouped according to their respective sampling events. HVR2 DNA from ten samples was sequenced and screened producing an average read length of approximately 423 bp. Of this, a 300 bp portion (i.e., from base position 50 to 350) was used for variant detection since SNPs in this region have been well documented ^[22]. A total of 19 distinct variants were detected with frequency greater than 5% of the total number of unique reads, all of which are present in MITOMAP—database of mtDNA Control Region Sequence Variants ^[22]. We observed some SNPs that were common to all samples with varying frequencies, while other SNPs were sample-specific, allowing each sample to have a unique human mtDNA signature in the form SNP allelic frequencies (Figure 2). The variation in SNP frequencies could be the result of several factors including limited sample size, population changes related to migration, changes in sampling time and storm runoff volumes during wet weather events, or a combination of them. Variants 73G and 263G were detected in all samples with high frequency ( > 90%), while variant 150T was detected in all samples except the samples belonging to sites 4 and 8. Interestingly, variant 263G has been observed for mitochondrial genomes from European populations ^[26], and is compatible with the European ancestry that originated in the island over five centuries ago. All the variants detected at site 6 were also detected for samples belonging to site 7, except for 232G which was detected at site 6 with relatively low frequency. This is expected since sites 6 and 7 are located in close proximity to each other. Site 8 had two unique variants (67T, 81T) which may be attributed to the influx of water from several different tributaries, since site 8 is located right before the RGA drains into the Atlantic Ocean at sea level and it is the most downstream of all sampling sites.

---

3.2. Haplotypes and haplogroup classification

Mitochondrial haplogroups have arisen from mutation and migration during human evolution and largely correspond to the geographic regions of their origin ^[23,24]. These mitochondrial haplogroups can be used to define ancestry based on the frequency of observation as a means of investigating population diversity ^[4,27]. For instance, there is broad correspondence between the L haplogroups and African ethnicity assignments, while the H haplogroups are most common among the Europeans. Consequently, we sought to use our human mtDNA sequences to extract haplotypes and classify them into haplogroups by comparing them to the Phylotree database ^[24]. We observed abundant diversity of haplotypes from HVR2 amplicons for all samples, which is consistent with the clear indication of human-associated pollution in the watershed ^[18]. The major haplotypes obtained from each sample are presented in Table 2.

The mitochondrial sequences were compared and assigned to haplogroups based on the differences in HVR2 sequence mutations with respect to the rCRS. Since most accurate haplogroup prediction is based on full mtDNA sequences, sequences were assigned to the closest haplogroup for which the HVR2 sequence contained all mutations that define the haplogroup. The most salient features of the haplogroup distribution (Figure 3) in the clustered sequences were the relatively high frequencies of haplogroup H (32%). This haplogroup is very common in Europe ^[28,29] and its presence in the mtDNA sequences from our study is in agreement with the presence of European population on the island. Other dominant haplogroups were T (25%), L (24%) and B (11%). As an additional verification step, HVR2 PCR products were cloned (TOPO TA Cloning Kit for Sequencing, Invitrogen, Carlsbad, CA) and 90 colonies were randomly picked and sent for Sanger sequencing. Nucleotide sequences were assembled and edited by using Sequencher 4.7 software (Gene Codes, Ann Arbor, MI) and analyzed for haplogroup prediction using MITOMASTER. Most of the sequences belonged to the haplogroup H (40%), followed by T (20%), L (12%) and B (10%). The results obtained with Sanger sequencing corresponded well with Illumina high-throughput sequencing supporting the reproducibility of the results by alternative sequencing methods.

---

3.3. Population diversity

To further explore the applicability of our HVR2-derived haplogroup data to local population diversity, several mitochondrial databases and studies were consulted to assign haplogroups to the general population groups found in Puerto Rico. We assigned haplogroups H, T and J to "West Eurasians"; haplogroup L to "Sub-Saharan African"; and haplogroup B to "American Indian" according to Martínez‐Cruzado et al. ^[30]. Based on the average distribution of HVR2-derived population groups, most mtDNA haplogroups were identified as of West Eurasian ancestry (57.6%), followed by those of African (23.9%) and American Indian (11%) ancestries (Figure 4). According to U.S. census data for 2010 ^[31], populations belonging to these groups live in and around the study area. Figure 5 presents the comparative analysis of the population data obtained through the two strategies— census data for population (by race) viz-a-viz the HVR2-derived population groups for three different locations in the watershed. There was a strong correlation between the federal census data and the mitochondrial haplogroups as an indicator of population composition (Pearson product-moment correlation coefficient, r = 0.9) demonstrating the suitability of human mitochondrial sequences to infer the population structure of the neighborhoods impacting the watershed.

While their relative abundance is different, the average census abundance patterns (White > African American > American Indian) are similar to our findings suggesting that the results correspond with the census data for population (by race). The mtDNA sequencing analysis suggests that American Indian ancestry is more prevalent than that the census data reports. Similarly, results from studies using HVR1 and other mtDNA-restriction profiles have also suggested that the presence of American Indian signals in Puerto Rico is more prevalent than previously considered ^[30], which is in agreement with our findings. The HVR2 motifs that are characteristic of the 'American Indian' haplogroups detected in this study are 73G, 143A and 263G (32) whereas Martínez‐Cruzado (30) used predetermined restriction motifs as defining markers, along with HVR1 sequences to resolve inconclusive results. The latter approach to define haplogroups is more exact since it is based on haplogroup-defining markers for the entire mtDNA and not only just HVR2. However, classification of haplogroups based on analysis of small mtDNA regions with maximal discriminative power is useful for environmental studies due to concerns related to DNA damage in the environment. This approach has proven useful in past anthropological studies involving analysis of Neanderthal-type specimen to sequence small regions (300–350 bp) of Neanderthal mtDNA ^[33,34]. Deducing population diversity from mtDNA sequences retrieved from waste streams may be more accurate than census data since these are limited to people who respond to surveys and are subject to misclassification of self-declared racial/ethnic background while waste streams are impacted by everyone connected to the public sewer system. However, it is also possible that certain groups are overrepresented using the current approach either because they disproportionately use the water resources, are not connected to the sewer pipelines (e.g., use of septic tanks) and/or live in close proximity to sampling sites than others. Signature sequences from areas impacted by leaky septic tanks and combined sewer overflows will also be reflected in these types of molecular surveys. While further studies are needed to better understand how all these different sources may impact haplogroup distribution, we suggest that the use of these methods could provide complementary information in epidemiological studies.

The overall bioinformatics strategy in this study included the following steps: (ⅰ) trim/clean sequencing reads and group them according to sites, (ⅱ) map sample specific reads to the rCRS, (ⅲ) annotate the mapped sequences to detect variants in HVR2 region, and (ⅳ) extract haplotypes from individual reads and assign haplogroups based on HVR2 sequence motifs. As reported here, next-generation sequencing technology of the mitochondrial hypervariable sequences enabled the identification of a great number of mtDNA variants and at varied allele frequencies. The methods used in our study for haplogrouping uses only HVR2 sequences which may result in coarse haplogroup assignments. Since most accurate haplogroup prediction is based on full mtDNA sequences, sequences were assigned to the closest haplogroup for which the HVR2 sequences contain all SNPs that define the haplogroup. We believe that future global sequencing efforts associated with distinct populations will provide improved phylogenetic resolution of the human mtDNA hypervariable regions as a tool for defining genetic ancestry. It has not escaped our attention that extending our technique to include other mtDNA regions and/or assembling full mitochondrial genomes through metagenomics approaches on a massively parallel scale would allow for tracking humans through public waste streams, thus ethical concerns remain an important consideration in future work.

Mitochondrial DNA analysis has been applied in several biomedical investigations of human evolution, for example, studies tracing the origin of modern humans or of certain human populations. In addition, mtDNA analysis is extremely effective in a forensic setting for the identification of criminals and victims of crimes or accidents. Although our study was confined to analysis of HVRII region of human mtDNA for samples collected in a limited number of geographic locations, it can be inferred that by targeting specific regions of mtDNA, we can estimate cancer rates, occurrence of diseases, and population diversity in watershed regions impacted by human contamination. Moreover, we envision that a similar approach could be used to study the population diversity of different animal species in natural settings, such as local versus migratory birds.

---

4. Conclusions

We investigated the occurrence of HVR2 allelic frequencies of human mtDNA derived from water samples taken within a fecally impacted tropical watershed. The SNPs within the human HVR2 sequences represented a unique molecular signature for evaluating anthropogenic site-specific inputs. We observed several HVR2 haplotypes linked to these samples, and used this haplogroup data to derive human population diversity within the different sites of the watershed. There was a strong correspondence between the demographic census data and the population composition based on mitochondrial haplogroups, demonstrating the suitability of human mitochondrial sequences to infer the population structure of the neighborhoods impacting the watershed. As the levels of human mtDNA is significantly high in point and non-point sources of fecal pollution, detecting mtDNA allelic signatures in environmental waters provides a unique approach for simultaneously studying fecal waste source tracking, human population diversity and other many anthropological investigations.

---

Acknowledgements

We would like to thank Mehdi Keddache for help in data analysis, and David Wendell for providing access to the CLC Genomics program. VK was supported by U. S. Environmental Protection Agency (EPA) via a post-doctoral appointment administered by the Oak Ridge Institute for Science and Education through an interagency agreement between the U.S. Department of Energy and the U.S. EPA. The manuscript has been subjected to the EPA's peer review and has been approved as an EPA publication. Mention of trade names or commercial products does not constitute endorsement or recommendation by the EPA for use. The views expressed in this article are those of the authors and do not necessarily represent the views or policies of the U.S. EPA.

---

Conflict of interest

All authors declare no conflicts of interest in this paper.

---