High Sequence Coverage Research Articles

Genomic history of the Italian population recapitulates key evolutionary dynamics of both Continental and Southern Europeans

BackgroundThe cline of human genetic diversity observable across Europe is recapitulated at a micro-geographic scale by variation within the Italian population. Besides resulting from extensive gene flow, this might be ascribable also to local adaptations to diverse ecological contexts evolved by people who anciently spread along the Italian Peninsula. Dissecting the evolutionary history of the ancestors of present-day Italians may thus improve the understanding of demographic and biological processes that contributed to shape the gene pool of European populations. However, previous SNP array-based studies failed to investigate the full spectrum of Italian variation, generally neglecting low-frequency genetic variants and examining a limited set of small effect size alleles, which may represent important determinants of population structure and complex adaptive traits. To overcome these issues, we analyzed 38 high-coverage whole-genome sequences representative of population clusters at the opposite ends of the cline of Italian variation, along with a large panel of modern and ancient Euro-Mediterranean genomes.ResultsWe provided evidence for the early divergence of Italian groups dating back to the Late Glacial and for Neolithic and distinct Bronze Age migrations having further differentiated their gene pools. We inferred adaptive evolution at insulin-related loci in people from Italian regions with a temperate climate, while possible adaptations to pathogens and ultraviolet radiation were observed in Mediterranean Italians. Some of these adaptive events may also have secondarily modulated population disease or longevity predisposition.ConclusionsWe disentangled the contribution of multiple migratory and adaptive events in shaping the heterogeneous Italian genomic background, which exemplify population dynamics and gene-environment interactions that played significant roles also in the formation of the Continental and Southern European genomic landscapes.

Transcriptomic Data Sets To Determine Gene Expression Changes Mediated by the Presence of PBT2 in Growth Medium of Multidrug-Resistant Neisseria gonorrhoeae WHO Z.

Neisseria gonorrhoeae causes the sexually transmitted infection gonorrhea. High-coverage (∼3,300-fold) transcriptome sequencing data have been collected from multidrug-resistant N. gonorrhoeae strain WHO Z grown in the presence and absence of PBT2.

Neandertal DNA in Africans

Sophisticated analyses of current human populations compared to a high-coverage Neandertal genome sequence indicate that, contrary to the previous consensus, African genomes carry a small but significant amount of Neandertal-specific DNA. This indicates back-migration into Africa of modern humans (carrying some Neandertal sequences) and underlines the complexity of ancient human migrations.

Insights into human genetic variation and population history from 929 diverse genomes.

Genome sequences from diverse human groups are needed to understand the structure of genetic variation in our species and the history of, and relationships between, different populations. We present 929 high-coverage genome sequences from 54 diverse human populations, 26 of which are physically phased using linked-read sequencing. Analyses of these genomes reveal an excess of previously undocumented common genetic variation private to southern Africa, central Africa, Oceania, and the Americas, but an absence of such variants fixed between major geographical regions. We also find deep and gradual population separations within Africa, contrasting population size histories between hunter-gatherer and agriculturalist groups in the past 10,000 years, and a contrast between single Neanderthal but multiple Denisovan source populations contributing to present-day human populations.

PyBSASeq: a simple and effective algorithm for bulked segregant analysis with whole-genome sequencing data

BackgroundBulked segregant analysis (BSA), coupled with next-generation sequencing, allows the rapid identification of both qualitative and quantitative trait loci (QTL), and this technique is referred to as BSA-Seq here. The current SNP index method and G-statistic method for BSA-Seq data analysis require relatively high sequencing coverage to detect significant single nucleotide polymorphism (SNP)-trait associations, which leads to high sequencing cost.ResultsWe developed a simple and effective algorithm for BSA-Seq data analysis and implemented it in Python; the program was named PyBSASeq. Using PyBSASeq, the significant SNPs (sSNPs), SNPs likely associated with the trait, were identified via Fisher’s exact test, and then the ratio of the sSNPs to total SNPs in a chromosomal interval was used to detect the genomic regions that condition the trait of interest. The results obtained this way are similar to those generated via the current methods, but with more than five times higher sensitivity. This approach was termed the significant SNP method here.ConclusionsThe significant SNP method allows the detection of SNP-trait associations at much lower sequencing coverage than the current methods, leading to ~ 80% lower sequencing cost and making BSA-Seq more accessible to the research community and more applicable to the species with a large genome.

A functional network of gastric-cancer-associated splicing events controlled by dysregulated splicing factors.

Comprehensive understanding of aberrant splicing in gastric cancer is lacking. We RNA-sequenced 19 gastric tumor–normal pairs and identified 118 high-confidence tumor-associated (TA) alternative splicing events (ASEs) based on high-coverage sequencing and stringent filtering, and also identified 8 differentially expressed splicing factors (SFs). The TA ASEs occurred in genes primarily involved in cytoskeletal organization. We constructed a correlative network between TA ASE splicing ratios and SF expression, replicated it in independent gastric cancer data from The Cancer Genome Atlas and experimentally validated it by knockdown of the nodal SFs (PTBP1, ESRP2 and MBNL1). Each SF knockdown drove splicing alterations in several corresponding TA ASEs and led to alterations in cellular migration consistent with the role of TA ASEs in cytoskeletal organization. We have therefore established a robust network of dysregulated splicing associated with tumor invasion in gastric cancer. Our work is a resource for identifying oncogenic splice forms, SFs and splicing-generated tumor antigens as biomarkers and therapeutic targets.

Ultrafast enzymatic digestion of proteins by microdroplet mass spectrometry

Enzymatic digestion for protein sequencing usually requires much time, and does not always result in high sequence coverage. Here we report the use of aqueous microdroplets to accelerate enzymatic reactions and, in particular, to improve protein sequencing. When a room temperature aqueous solution containing 10 µM myoglobin and 5 µg mL−1 trypsin is electrosonically sprayed (−3 kV) from a homemade setup to produce tiny (∼9 µm) microdroplets, we obtain 100% sequence coverage in less than 1 ms of digestion time, in sharp contrast to 60% coverage achieved by incubating the same solution at 37 °C for 14 h followed by analysis with a commercial electrospray ionization source that produces larger (∼60 µm) droplets. We also confirm the sequence of the therapeutic antibody trastuzumab (∼148 kDa), with a sequence coverage of 100% for light chains and 85% for heavy chains, demonstrating the practical utility of microdroplets in drug development.

Noninvasive chimeric DNA profiling identifies tumor-originated HBV integrants contributing to viral antigen expression in liver cancer.

Host genome integration of HBV sequence is considered to be significant in HBV antigen expression and the development of hepatocellular carcinoma (HCC). We developed a probe-based capture strategy to enrich integrated HBV DNA for deep-sequencing analysis of integration sites in paired patient samples derived from tumor, liver tissue adjacent to tumor, saliva and plasma, as a platform for exploring the correlation, significance and utility of detecting integrations in these sample types. Most significantly, alpha fetoprotein levels significantly correlated to the amounts of integrations detected in tumor. Viral-host chimeric DNA fragments were successfully detected at high sequencing coverage in plasma rather than saliva samples from HCC patients, and each fragment of this type was only seen once in plasma from chronic hepatitis B patients. Almost all plasma chimeric fragments were derived from integrations in tumor rather than in adjacent liver tissues. Over 50% of them may produce viral-host chimeric transcripts according to deep RNA sequencing in paired tissue samples. Particularly, in patients with low HBV DNA level (< 250 UI/ml), the seemingly normal HBsAg titers may be explained by larger amounts of integrations detected. Meanwhile, we developed a strategy to predict integrants by pairing breakpoints for each integration event. Among four resolved viral patterns, the majority of Pattern I events (81.2%) retained the complete opening reading frame for HBV surface proteins. We achieve the efficient enrichment of plasma cell-free chimeric DNA from integration site, and demonstrate that chimeric DNA profiling in plasma is a promising noninvasive approach to monitor HBV integration in liver cancer development and to determine the ability of integrated sequences to express viral proteins that can be targeted, e.g. by immunotherapies.

Chromosome-level assembly of wild Bactrian camel genome reveals organization of immune gene loci.

Camelids are characterized by their unique adaptive immune system that exhibits the generation of homodimeric heavy-chain immunoglobulins, somatic hypermutation of T-cell receptors, and low genetic diversity of major histocompatibility complex (MHC) genes. However, short-read assemblies are typically highly fragmented in these gene loci owing to their repetitive and polymorphic nature. Here, we constructed a chromosome-level assembly of wild Bactrian camel genome based on high-coverage long-read sequencing and chromatin interaction mapping. The assembly with a contig N50 of 5.37Mb and a scaffold N50 of 76.03Mb, represents the most contiguous camelid genome to date. The genomic organization of immunoglobulin heavy-chain locus was similar between the wild Bactrian camel and alpaca, and genes encoding for conventional and heavy-chain antibodies were intermixed. The organizations of two immunoglobulin light-chain loci and four T cell receptor loci were also fully deciphered using the new assembly. Additionally, the complete classical MHC region was resolved into a single contig. The high-quality assembly presented here provides an essential reference for future investigations examining the camelid immune system.

Contribution of common and rare variants to bipolar disorder susceptibility in extended pedigrees from population isolates

Current evidence from case/control studies indicates that genetic risk for psychiatric disorders derives primarily from numerous common variants, each with a small phenotypic impact. The literature describing apparent segregation of bipolar disorder (BP) in numerous multigenerational pedigrees suggests that, in such families, large-effect inherited variants might play a greater role. To identify roles of rare and common variants on BP, we conducted genetic analyses in 26 Colombia and Costa Rica pedigrees ascertained for bipolar disorder 1 (BP1), the most severe and heritable form of BP. In these pedigrees, we performed microarray SNP genotyping of 838 individuals and high-coverage whole-genome sequencing of 449 individuals. We compared polygenic risk scores (PRS), estimated using the latest BP1 genome-wide association study (GWAS) summary statistics, between BP1 individuals and related controls. We also evaluated whether BP1 individuals had a higher burden of rare deleterious single-nucleotide variants (SNVs) and rare copy number variants (CNVs) in a set of genes related to BP1. We found that compared with unaffected relatives, BP1 individuals had higher PRS estimated from BP1 GWAS statistics (P = 0.001 ~ 0.007) and displayed modest increase in burdens of rare deleterious SNVs (P = 0.047) and rare CNVs (P = 0.002 ~ 0.033) in genes related to BP1. We did not observe rare variants segregating in the pedigrees. These results suggest that small-to-moderate effect rare and common variants are more likely to contribute to BP1 risk in these extended pedigrees than a few large-effect rare variants.

Abstract PD8-07: A pro-metastatic alternative splicing program that drives breast cancer progression

Abstract Background: Breast cancer progression accompanies broad reprogramming of the alternative splicing landscape with functional implications for tumor growth and metastasis. However, the underlying regulatory pathways that are hijacked by cancer cells to drive pathological splicing programs remain largely unknown. Here, we have described a novel regulatory pathway, mediated by the splicing factor SNRPA1, that promotes breast cancer progression. Results: We used high-coverage RNA sequencing to compare alternative splicing patterns in highly metastatic breast cancer cells (MDA-LM2) to their less metastatic parental line (MDA-231). We identified a novel structural splicing enhancer (SSE) element that is enriched among the exons up-regulated in highly metastatic cells and showed that it is bound by the RNA-binding protein SNRPA1, which in turn acts as a splicing factor. We used loss-of-function and gain-of-function both in vitro and in vivo xenograft mouse models to establish SNRPA1 as a promoter of breast cancer cell invasion (~2-fold change in two independent cell lines, P=0.002) and metastasis (~5-10 fold change in lung colonization in two independent cell lines, P&amp;lt;0.01). We also identified specific splicing events downstream of SNRPA1 that are required for its pro-metastatic activity. Single and multivariate analyses of publicly available datasets established that the higher expression of SNRPA1 is significantly associated with metastatic progression (METABRIC: n=1456, P&amp;lt;1e-5; EXPO: n=243, P=0.002), as well as reduced metastasis-free (kmplot: n=1746, HR=1.4, P&amp;lt;1e-3) and disease-free survival (METABRIC: n=1979, HR=1.7, P&amp;lt;1e-6; multivariate: HR=1.3, P=0.01). Finally, we measured the expression of SNRPA1 and its target exon in 96 samples from normal tissue biopsy to stage IV metastatic disease and observed a concordant increase in the activity of this novel pro-metastatic splicing program (ANOVA: P&amp;lt;1e-10; U test stage IV to N: P&amp;lt;1e-3). Methods: We used a framework called ‘mixture of isoforms’ (MISO) to quantify changes in alternative splicing. We then applied iTEISER, a computational tool we have developed, to find cis-regulatory elements that drive alternative splicing patterns. We used bait RNA to co-precipitate and mass-spectrometry to identify the RNA-binding proteins that interact with SSEs. Lentiviral delivery of shRNAs (two independent hairpins) or transgenes to modulate SNRPA1 expression and RNA-seq was used to measure their regulatory consequences on splicing. For the resulting lines (both MDA-231 and HCC-1806 backgrounds), in vitro cancer cell invasion and in vivo lung colonization assays were also performed (n=5 mice in each arm). Log-rank test (univariate) and Cox Proportional Hazard Models (multivariate) were used to perform survival analyses in METABRIC and the kmplot aggregate dataset. Mann-Whitney or ANOVA was used to compare expression of SNRPA1 across samples stratified based on sub-type and tumor grade/stage in public datasets as well as our own measurements in clinical samples (n=96; 5 healthy, 23 stage I, 30 stage II, 29 stage III, and 9 stage IV). Conclusion: In sum, we have discovered a previously unknown regulatory pathway of RNA splicing that acts as a promoter of metastasis in breast, which may inform new therapeutic targets among the SNRPA1 regulon for adjuvant therapy. Citation Format: Hani Goodarzi, Lisa Fish, Bruce Culbertson. A pro-metastatic alternative splicing program that drives breast cancer progression [abstract]. In: Proceedings of the 2019 San Antonio Breast Cancer Symposium; 2019 Dec 10-14; San Antonio, TX. Philadelphia (PA): AACR; Cancer Res 2020;80(4 Suppl):Abstract nr PD8-07.

Selective sweeps on novel and introgressed variation shape mimicry loci in a butterfly adaptive radiation.

Natural selection leaves distinct signatures in the genome that can reveal the targets and history of adaptive evolution. By analysing high-coverage genome sequence data from 4 major colour pattern loci sampled from nearly 600 individuals in 53 populations, we show pervasive selection on wing patterns in the Heliconius adaptive radiation. The strongest signatures correspond to loci with the greatest phenotypic effects, consistent with visual selection by predators, and are found in colour patterns with geographically restricted distributions. These recent sweeps are similar between co-mimics and indicate colour pattern turn-over events despite strong stabilising selection. Using simulations, we compare sweep signatures expected under classic hard sweeps with those resulting from adaptive introgression, an important aspect of mimicry evolution in Heliconius butterflies. Simulated recipient populations show a distinct ‘volcano’ pattern with peaks of increased genetic diversity around the selected target, characteristic of sweeps of introgressed variation and consistent with diversity patterns found in some populations. Our genomic data reveal a surprisingly dynamic history of colour pattern selection and co-evolution in this adaptive radiation.

De novo mutations across 1,465 diverse genomes reveal mutational insights and reductions in the Amish founder population

De novo mutations (DNMs), or mutations that appear in an individual despite not being seen in their parents, are an important source of genetic variation whose impact is relevant to studies of human evolution, genetics, and disease. Utilizing high-coverage whole-genome sequencing data as part of the Trans-Omics for Precision Medicine (TOPMed) Program, we called 93,325 single-nucleotide DNMs across 1,465 trios from an array of diverse human populations, and used them to directly estimate and analyze DNM counts, rates, and spectra. We find a significant positive correlation between local recombination rate and local DNM rate, and that DNM rate explains a substantial portion (8.98 to 34.92%, depending on the model) of the genome-wide variation in population-level genetic variation from 41K unrelated TOPMed samples. Genome-wide heterozygosity does correlate with DNM rate, but only explains <1% of variation. While we are underpowered to see small differences, we do not find significant differences in DNM rate between individuals of European, African, and Latino ancestry, nor across ancestrally distinct segments within admixed individuals. However, we did find significantly fewer DNMs in Amish individuals, even when compared with other Europeans, and even after accounting for parental age and sequencing center. Specifically, we found significant reductions in the number of C→A and T→C mutations in the Amish, which seem to underpin their overall reduction in DNMs. Finally, we calculated near-zero estimates of narrow sense heritability (h2), which suggest that variation in DNM rate is significantly shaped by nonadditive genetic effects and the environment.

Spitz melanoma is a distinct subset of spitzoid melanoma.

Melanomas that have histopathologic features that overlap with those of Spitz nevus are referred to as spitzoid melanomas. However, the diagnostic concept is used inconsistently and genomic analyses suggest it is a heterogeneous category. Spitz tumors, the spectrum of melanocytic neoplasms extending from Spitz nevi to their malignant counterpart Spitz melanoma, are defined in the 2018 WHO classification of skin tumors by the presence of specific genetic alterations such a kinase fusions or HRAS mutations. It is unclear what fraction of ‘spitzoid melanomas’ defined solely by their histopathologic features belong to the category of Spitz melanoma or to other melanoma subtypes. We assembled a cohort of 25 spitzoid melanomas diagnosed at a single institution over an eight-year period and performed high coverage DNA sequencing of 480 cancer related genes. Transcriptome wide RNA sequencing was performed for select cases. Only 9 cases (36%) had genetic alterations characteristic of Spitz melanoma, including HRAS mutation or fusion involving BRAF, ALK, NTRK1, or MAP3K8. The remaining cases were divided into those with a MAPK activating mutation and those without. Both Spitz melanoma and spitzoid melanomas in which a MAPK activating mutation could not be identified tended to occur in younger patients on skin with little solar elastosis, infrequently harbored TERT promoter mutations, and had a lower burden of pathogenic mutations than spitzoid melanomas with non-Spitz MAPK activating mutations. The MAPK activating mutations identified affected non-V600 residues of BRAF as well as NRAS, MAP2K1/2, NF1 and KIT while BRAF V600 mutations, the most common mutations in melanomas of the WHO low-CSD category, were entirely absent.While the ‘spitzoid melanomas’ comprising our cohort were enriched for bona fide Spitz melanomas, the majority of melanomas fell outside of the genetically defined category of Spitz melanomas, indicating that histomorphology is an unreliable predictor of Spitz lineage.

Establishment of a highly precise multi-attribute method for the characterization and quality control of therapeutic monoclonal antibodies

ABSTRACT The multi-attribute method (MAM) has garnered attention as a new quality control method of therapeutic monoclonal antibodies (mAbs). MAM analysis allows multiple relative quantifications of several structural attributes of therapeutic mAbs; however, some issues remain to be addressed in its procedures especially for sample preparation. The goal of this study was to optimize the sample preparation method for MAM analysis of mAbs. Using a model mAb, we compared five sample preparation methods based on sequence coverage, peptide redundancy, missed cleavage and chemical deamidation. It was found that low pH buffer and short digestion time reduced artificial deamidation. The desalting process after carboxymethylation was essential to obtaining high sequence coverage by a short digestion time. The generation of missed cleavage peptides was also improved by using a trypsin/lysyl endopeptidase (Lys-C) mixture. Next, we evaluated the usefulness of our method as a part of MAM analysis. Finally, 17 glycopeptides, 2 deamidated peptides and N- and C-terminal peptides of the heavy chain were successfully monitored with acceptable mass accuracy and coefficient of variation (CV, %) of the relative peak area. On the other hand, 4 oxidated peptides indicated the unavoidable slightly higher inter-assay CV (%) of the peak area ratio due to the instability in the MS sample solution. Collectively, we demonstrated that our method was applicable as an easy and reliable sample preparation method for MAM analysis, and the variation in the relative peak area could be influenced by the modification type rather than by the amount of each peptide.

Predicted Archaic 3D Genome Organization Reveals Genes Related to Head and Spinal Cord Separating Modern from Archaic Humans.

High coverage sequences of archaic humans enabled the reconstruction of their DNA methylation patterns. This allowed comparing gene regulation between human groups, and linking such regulatory changes to phenotypic differences. In a previous work, a detailed comparison of DNA methylation in modern humans, archaic humans, and chimpanzees revealed 873 modern human-derived differentially methylated regions (DMRs). To understand the regulatory implications of these DMRs, we defined differentially methylated genes (DMGs) as genes that harbor DMRs in their promoter or gene body. While most of the modern human-derived DMRs could be linked to DMGs, many others remained unassigned. Here, we used information on 3D genome organization to link ~70 out of the remaining 288 unassigned DMRs to genes. Combined with the previously identified DMGs, we reinforce the enrichment of these genes with vocal and facial anatomy, and additionally find significant enrichment with the spinal column, chin, hair, and scalp. These results reveal the importance of 3D genomic organization in understanding gene regulation by DNA methylation.

Accurate Allele Frequencies from Ultra-low Coverage Pool-Seq Samples in Evolve-and-Resequence Experiments.

Evolve-and-resequence (E+R) experiments leverage next-generation sequencing technology to track the allele frequency dynamics of populations as they evolve. While previous work has shown that adaptive alleles can be detected by comparing frequency trajectories from many replicate populations, this power comes at the expense of high-coverage (>100x) sequencing of many pooled samples, which can be cost-prohibitive. Here, we show that accurate estimates of allele frequencies can be achieved with very shallow sequencing depths (<5x) via inference of known founder haplotypes in small genomic windows. This technique can be used to efficiently estimate frequencies for any number of bi-allelic SNPs in populations of any model organism founded with sequenced homozygous strains. Using both experimentally-pooled and simulated samples of Drosophila melanogaster, we show that haplotype inference can improve allele frequency accuracy by orders of magnitude for up to 50 generations of recombination, and is robust to moderate levels of missing data, as well as different selection regimes. Finally, we show that a simple linear model generated from these simulations can predict the accuracy of haplotype-derived allele frequencies in other model organisms and experimental designs. To make these results broadly accessible for use in E+R experiments, we introduce HAF-pipe, an open-source software tool for calculating haplotype-derived allele frequencies from raw sequencing data. Ultimately, by reducing sequencing costs without sacrificing accuracy, our method facilitates E+R designs with higher replication and resolution, and thereby, increased power to detect adaptive alleles.

The Genome Sequence of the Eastern Woodchuck (Marmota monax) – A Preclinical Animal Model for Chronic Hepatitis B

The Eastern woodchuck (Marmota monax) has been extensively used in research of chronic hepatitis B and liver cancer because its infection with the woodchuck hepatitis virus closely resembles a human hepatitis B virus infection. Development of novel immunotherapeutic approaches requires genetic information on immune pathway genes in this animal model. The woodchuck genome was assembled with a combination of high-coverage whole-genome shotgun sequencing of Illumina paired-end, mate-pair libraries and fosmid pool sequencing. The result is a 2.63 Gigabase (Gb) assembly with a contig N50 of 74.5 kilobases (kb), scaffold N50 of 892 kb, and genome completeness of 99.2%. RNA sequencing (RNA-seq) from seven different tissues aided in the annotation of 30,873 protein-coding genes, which in turn encode 41,826 unique protein products. More than 90% of the genes have been functionally annotated, with 82% of them containing open reading frames. This genome sequence and its annotation will enable further research in chronic hepatitis B and hepatocellular carcinoma and contribute to the understanding of immunological responses in the woodchuck.

LinkedSV for detection of mosaic structural variants from linked-read exome and genome sequencing data

Linked-read sequencing provides long-range information on short-read sequencing data by barcoding reads originating from the same DNA molecule, and can improve detection and breakpoint identification for structural variants (SVs). Here we present LinkedSV for SV detection on linked-read sequencing data. LinkedSV considers barcode overlapping and enriched fragment endpoints as signals to detect large SVs, while it leverages read depth, paired-end signals and local assembly to detect small SVs. Benchmarking studies demonstrate that LinkedSV outperforms existing tools, especially on exome data and on somatic SVs with low variant allele frequencies. We demonstrate clinical cases where LinkedSV identifies disease-causal SVs from linked-read exome sequencing data missed by conventional exome sequencing, and show examples where LinkedSV identifies SVs missed by high-coverage long-read sequencing. In summary, LinkedSV can detect SVs missed by conventional short-read and long-read sequencing approaches, and may resolve negative cases from clinical genome/exome sequencing studies.

Neutrophil myeloperoxidase harbors distinct site-specific peculiarities in its glycosylation

Anti-neutrophil cytoplasmic autoantibodies (ANCAs) are directed against lysosomal components of neutrophils. ANCAs directed to proteinase 3 and myeloperoxidase (MPO) in particular are associated with distinct forms of small vessel vasculitides. MPO is an abundant neutrophil-derived heme protein that is part of the antimicrobial defense system. The protein is typically present in the azurophilic granules of neutrophils, but a large portion may also enter the extracellular space. It remains unclear why MPO is frequently the target of antibody-mediated autoimmune responses. MPO is a homodimeric glycoprotein, posttranslationally modified with complex sugars at specific sites. Glycosylation can strongly influence protein function, affecting its folding, receptor interaction, and backbone accessibility. MPO potentially can be heavily modified as it harbors 5 putative N-glycosylation sites (10 in the mature dimer). Although considered important for MPO structure and function, the full scope and relative abundance of the glycans attached to MPO is unknown. Here, combining bottom-up glycoproteomics and native MS approaches, we structurally characterized MPO from neutrophils of healthy human donors. We quantified the relative occupancy levels of the glycans at each of the five sites and observed complex heterogeneity and site-specific glycosylation. In particular, we detected glycosylation phenotypes uncommon for glycoproteins in the extracellular space, such as a high abundance of phosphorylated high-mannose species and severely truncated small glycans having the size of paucimannose or smaller. We hypothesize that the atypical glycosylation pattern found on MPO might contribute to its specific processing and presentation as a self-antigen by antigen-presenting cells.