Mitochondrial DNA (mtDNA) has an important, yet often overlooked, role in health and disease. Constraint models quantify depletion of genetic variation by negative selection, representing a powerful tool for identifying deleterious variation underlying human phenotypes1-4. However, a constraint model for the mtDNA had not been developed, due to the unique features of this genome. Here we describe the development of a mitochondrial constraint model and its application to the Genome Aggregation Database (gnomAD), a large-scale population dataset reporting mtDNA variation across 56,434 humans5. Our results demonstrate strong depletion of expected variation, showing most pathogenic mtDNA variants remain undiscovered. To aid their identification, we compute constraint metrics for every protein, transfer RNA, and ribosomal RNA gene in the human mtDNA, as well as for non-coding regions. This includes assessment of regional or positional constraint within every gene, and local constraint for every base position. We identify enrichment of pathogenic variation in the most constrained sites, which include loci typically overlooked in mtDNA analyses, and show that the most constrained gene regions and non-coding elements encode functionally-critical sites. Lastly, we demonstrate how these metrics can improve the discovery of mtDNA variation underlying rare and common human phenotypes.