{"paper_id":"432639c6-50a0-4e8d-bd79-560436de2af7","body_text":"Structural variants are enriched for molecular QTL in a \ncattle long read cohort \n \nXena Marie Mapel1,*, Alexander S. Leonard1,*, Hubert Pausch1,# \n1 Animal Genomics, ETH Zurich, Universitaetstrasse 2, 8092 Zurich, Switzerland \n* equal contribution \n# Corresponding author: hubert.pausch@usys.ethz.ch \n \n \nAbstract \nSequencing mapping cohorts with long-read technology is crucial to understand the impact \nof structural variants (SVs) on complex traits. Here, we obtained 4.86 terabases of HiFi \nreads with an average read N50 of 16.3 Kb from 120 Bos taurus taurus bulls, yielding a \nmean coverage depth of 13.5-fold. We genotyped 23.8 M small variants and 79.3 k SVs to \nperform association testing with molecular phenotypes derived from a subset of 117 bulls \nwith total RNA sequencing data from testis tissue. We identified 27.3 k molecular QTL \nincluding 316 for which SVs were the top variants. This corresponds to a 2.1- and 5.6-fold \nenrichment of SVs among expression and splicing QTL, respectively. When considering SVs \nin perfect LD with the lead small variant, the enrichment increased to 6.1- and 12-fold for \nexpression and splicing QTL, respectively. Imperfect genotyping for large SVs and other \nvariants limited our ability to detect all SV top variants, suggesting that the true enrichment of \nSVs among molecular QTL may be even higher. These results demonstrate that SVs have \nprofound impacts on gene expression and splicing variation in cattle but highlight the \nnecessity of improved SV genotyping to fully leverage long-read sequencing cohorts for \ndissecting complex traits. \n \nIntroduction \nUnderstanding the genetic basis of complex traits and diseases requires comprehensive \nexploration of genomic variation within and between individuals1,2. Genome-wide association \nstudies (GWAS) have identified numerous loci linked to specific traits, with single nucleotide \npolymorphisms (SNPs) derived from microarrays or short-read sequencing serving as the \npredominant source of genomic data3–5. Structural variants (SVs), specifically, insertions, \ndeletions, inversions, and duplications larger than 50 bp, contribute more to genomic \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nvariation between individuals than SNPs and other small variants6. Since short-read \nsequencing has inherent limitations in detecting and genotyping SVs, until recently, a \nsubstantial portion of genomic variation remained undercharacterized7–10. \nLong-read sequencing enables the analysis of DNA fragments orders of magnitudes longer \nthan short reads, offering the resolution necessary to identify and genotype SVs11–13. Studies \nin humans and plants have begun to explore SV diversity at the cohort level using long-read \nsequencing14–16. While these efforts remain limited in number, they have provided compelling \nevidence that SVs disproportionately contribute to phenotypic variation. These studies also \ndemonstrated that accurately genotyping large and complex SVs remains challenging, even \nwith long and highly precise reads17,18. \nResearch on SVs in cattle has primarily relied on pangenome analyses19–21, short-read \nsequencing22–24, and quantitative trait loci (QTL) fine-mapping25,26. These studies have \nrevealed several SVs with large impacts on phenotypic traits including polledness27, coat \ncolour variation26,28–30, milk production31, mastitis susceptibility32, and fertility31. The discovery \nof such impactful SVs highlights their biological and economical relevance underscoring the \ncritical need to systematically investigate them in genome-wide analyses. However, long-\nread sequencing at the cohort level has yet to be conducted in cattle. \nHere, we assess genetic variation in a cohort of 120 bulls from moderately covered HiFi \nreads. We identify and genotype 23.8 M small variants and 79.3 k SVs larger than 50 bp \nrepresenting the vast majority of genetic variants segregating in that population. We perform \nmolecular QTL mapping with deeply sequenced total RNA from testis tissue of the same \nindividuals, identifying more than 27 k gene expression and splicing QTL. Our findings show \nthat molecular QTL are enriched for SVs but also highlight that obtaining accurate genotypes \nfor structural variants remains challenging. \n \nResults \nHiFi sequencing, alignment, and variant calling of 120 cattle generates an exhaustive \nset of variants \nWe resequenced a biobanked cohort of 120 Bos taurus taurus samples of primarily \nBraunvieh (BV) ancestry with existing short sequencing reads33 with long reads to assess \nboth small variant and SV diversity. We collected 4.86 terabases of HiFi reads from 49 8M \nSMRT cells sequenced on Sequel IIe and 41 25M SMRT cells sequenced on Revio. The \nmean depth of coverage was 13.8-fold and 13.5-fold for the existing Illumina and newly \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\ncreated HiFi reads, respectively (Figure 1a), enabling a fair comparison for variant discovery. \nThe mean read N50 was 16.3 Kb with a mean Phred quality score of 33.7, demonstrating \nthe expected strong negative correlation between read length and quality for HiFi reads \n(Pearson’s r: -0.80, p=5.87 × 10-37; Supplementary Figure 1).  \n \nFigure 1. Comparison of alignment coverage and variant accuracy between Illumina and HiFi reads. (a) HiFi and \nIllumina sequencing depth was roughly comparable across bulls with breeds assigned as Original Braunvieh \n(OBV), Brown Swiss (BSW), ambiguous Braunvieh ancestry (BV), crosses between Braunvieh and a non-\nBraunvieh breed (Cross), or animals without Braunvieh ancestry (non-BV). (b) Fraction of chromosomes covered \nby at least two reads with MAPQ≥5 to enable variant calling. F1 score for SNPs (c) and indels (d), taking the HiFi \nvariants as truth. Variants are stratified by regions annotated as tandem repeats, transposable elements \n(SINE/LINE/LTR), or neither (normal). \n \nWe aligned the 120 HiFi and Illumina samples to the ARS-UCD2.0 Bos taurus reference \ngenome, which included the new T2T assembly of a Y chromosome from a Wagyu bull, \nwhereas the remaining genome is from a Hereford cow34. Alignment depth of both the short \nand long read sequencing was comparable across the autosomes, with only minor increase \nin the fraction of autosomal sequence covered by alignments suitable for variant calling for \nHiFi reads (99.4% versus 99.0%; Figure 1b). Conversely, the alignment improvement was \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nnoticeable for the X chromosome (90.0% versus 86.0%) and substantial for the Y \nchromosome (58.9% versus 32.0%). In the former case, this was largely due to the improved \nmapping of HiFi reads distinguishing high sequence similarity of X:19309231-28499819 to \nthe unplaced contig NKLS02002208.1, while the latter was due to a large fraction of \nsequence containing higher order repeats, which was almost completely inaccessible with \nshort read sequencing (Supplementary Figure 2). However, the male-specific region of the Y \nchromosome (MSY) still displayed uneven and inconsistent alignment with the HiFi reads \ndue to the repeat structure typically exceeding the length of the sequenced long reads which \nposes substantial challenges to unambiguous read alignment. \nWe called small variants with DeepVariant separately for HiFi and Illumina reads, again \nfinding minor differences in the number of detected variants and called genotypes on \nautosomes and larger differences on sex chromosomes (Table 1). The increased number of \nvariants called on X from the HiFi alignments was evenly distributed across the \nchromosome, while the larger number of variants discovered on Y primarily resulted from the \nincrease in bases covered with HiFi reads in the MSY, with a substantial increase in variants \ncalled in newly accessible regions that were almost entirely uncalled from short read \nsequencing. Given all the samples are male, we specified that variants could have \nheterozygous genotypes in the pseudo-autosomal region (PAR) of the Y chromosome, while \nthey could only be hemizygous in the X or MSY regions. However, there were “real” signals \nof heterozygous variants in the hemizygous region, primarily due to copy number variants of \namplicon genes (e.g., TSPY, HSFY2, and RBMY) relative to the reference, where the variant \nallele frequency could estimate the copy number in addition to coverage-based estimates \n(Supplementary Figure 3). \nWe assessed potential functional consequences of small variants with the Ensembl Variant \nEffect Predictor (Table 1), finding 247 and 166 biallelic SNPs with “HIGH” impacts \nrespectively private to HiFi- and Illumina-based variant calls, with 2,013 common to both. \nAfter manual inspection of 15 HIGH impact variants private to each group (Supplementary \nTable 1), we identified poor mapping of Illumina reads as the primary cause for the \ndiscrepancy between the variants called from short and long reads. All but one Illumina-only \nHIGH impact variant were likely genotyping artifacts as evidenced by implausible insert sizes \nand interchromosomal-mapping reads, with the sole exception of a singleton variant missed \nby the HiFi reads due to an uneven balance of alleles (7:1 read depth for the alleles). \nConversely, all HiFi-only HIGH impact variants appeared correct, mostly missed by Illumina \nreads due to mapping quality of 0 or highly diverged regions potentially leading to local \nreference bias (Supplementary Figure 4). \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nTable 1. Variants called by DeepVariant from the HiFi or Illumina alignments. The Y chromosome is split into the \nY PAR (where heterozygous genotypes are allowed) and MSY (where only hemizygous genotypes are allowed). \nThe number of variants identified by VEP as HIGH impact are given in parentheses. \n Autosomes X Y PAR MSY MT \nIllumina \n22,636,961 \n(2,198) \n383,835 \n(20) \n102,570 \n(4) \n20,177 \n(0) \n379 \n(8) \nHiFi \n23,163,376 \n(2,273) \n444,306 \n(24) \n112,420 \n(6) \n69,402 \n(1) \n386 \n(8) \nIncrease \n2.3% \n(3.4%) \n15.8% \n(20.0%) \n9.6% \n(50.0%) \n244.0% \n(NA) \n1.8% \n(0.0%) \n \nWe found generally high genotyping concordance for SNPs on the autosomes \n(F1=0.93±0.02), taking the HiFi variants as truth and short read variants as query, although \nthere was a moderate drop (F1=0.86±0.03) in genomic regions classified as tandem repeats \n(Figure 1c). There was no such drop for regions containing transposable elements, \nsuggesting that there was sufficient variation within different transposable elements for short \nreads to correctly align, while they struggle in long tandem repeats with little motif variation. \nVariants genotyped on the X chromosome (F1=0.84±0.04) and the PAR of the Y \nchromosome (F1=0.85±0.02) were less concordant across all types of regions, while the \nlimited number of short read-based variant calls in the MSY region led to a substantial drop \nin F1 score (F1=0.27±0.02). Indels behaved similarly, with an even more pronounced drop in \naccuracy around tandem repeats (Figure 1d). Given the observed improvements in small \nvariant genotyping from long reads, we only retained the HiFi-based small variants for \ndownstream use. \nWe also called structural variants from the HiFi read alignments with sniffles2, finding 75,164 \nand 4,111 SVs on the autosomes and sex chromosomes, respectively. We found \napproximately 23.9±0.6 k SVs per genome, with only on the order of 100s of novel SVs for \neach additional individual after 100 samples (Supplementary Figure 5). Given the \npredominant Braunvieh breed ancestry of the cohort, which has an estimated effective \npopulation size of 70, we likely captured the bulk of non-rare SVs present in this breed. SVs \nwere overrepresented in regions annotated as tandem repeats, totalling just 3.2% of the \ngenome, with 30.2% of SVs having at least half of their sequence contained within these \nregions. We observed a disproportionately large number of SVs of length ~300 bp, ~1.3 Kb, \n~5.5 Kb, and ~8.4 Kb (Figure 2a), corresponding to known transposable elements like Bov-\nA2 (SINE), BTLTR1B (LTR), and BovB (LINE). Given the average read N50 of 16 Kb for this \ncohort, there is likely a small number of longer SVs that remains unidentified (Supplementary \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nFigure 6). We further found that 7,944 SVs (10.02%) were present in only one sample (6,988 \nsingletons and 956 doubletons), ranging between 6 and 573 SVs unique to each sample. \nWe found 37,380 novel autosomal SVs compared to a previous autosome-only, pangenome-\nderived SV panel constructed from 15 Braunvieh haplotypes, including 19,507 SVs with a \nminor allele frequency (MAF) between 1% and 15% (Figure 2b), as well as 376 duplications \nand 422 inversions which cannot be easily genotyped through the previously applied k-mer \napproach35.  \n \nFigure 2. Structural variant calling in a long-read cohort. (a) SV size distribution across the four SV classes (DEL \n– deletion, INS – insertion, INV – inversion, DUP – duplication) examined. Large spikes in insertions and \ndeletions correspond to known transposable elements. (b) The majority of SVs which were only discovered \nthrough the long read cohort compared to a pangenome panel have low minor allele frequency, whereas those \npreviously discovered are typically common SVs. Allele frequency was calculated from the cohort, and so the \npanel-only SVs have no assigned allele frequency. (c) Forced genotyping of SVs substantially reduces \nmissingness to similar levels identified in SNP variant calling. \n \nDuring sniffles2 joint-calling, “similar” SVs were merged into single alleles, collapsing \nmultiallelic SVs. This overwhelmingly homogenised SVs that might only differ by small \ndifferences in length or position (Supplementary Figure 7). Approximately 82% of the merged \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nalleles had allele length standard deviations below 10, while only approximately 6% were \nabove 100. However, SVs overlapping tandem repeats were particularly affected by merging \nalleles, either incorrectly treating different tandem repeat allele lengths as distinct SVs \n(rather than distinct alleles) or treating different tandem repeat allele lengths as a single \nconsensus allele (Supplementary Figure 8). Several large tandem repeats contained multiple \nSVs, primarily differing only by the start coordinate (corresponding to inserting/removing \nidentical motifs at different locations within the tandemly duplicated sequence), which were \nnot merged. On the other hand, smaller tandem repeats tended to be merged too \naggressively, picking a consensus allele length which did not accurately reflect the copy-\nnumber diversity present in individual, pre-merged alleles. \nJoint-calling SVs across the entire cohort resulted in a relatively high proportion of missing \ngenotypes (mean of 6.8%), with duplications having the highest missingness (mean of \n9.9%). The force-calling of candidate SVs from the original alignments substantially reduced \nthe missingness to a mean of 1.2% (with the lowest mean missingness now observed in \nduplications at 0.7%), comparable to the missingness of joint-called SNPs (Figure 2c). \nGenotypes that remained missing after force-calling were significantly associated with low \nalignment coverage (p=5.0 × 10-23) and, to a much lesser extent, mean read length (p=2.9 × \n10-6) after conducting a Type II ANOVA. Missingness was elevated at the start and end of \nchromosomes, roughly corresponding to centromeric and telomeric regions in the \nacrocentric cattle autosomes (Supplementary Figure 9), with other peaks corresponding to \nregions known to be highly polymorphic (e.g., BoLA on BTA23) or challenging for alignment \n(e.g., large segmental duplication on BTA10). Approximately 4% of the SV genotypes \n(380,911 genotypes across 54,561 unique SVs) changed during force-calling. The \noverwhelming majority of these changes (89.7%) were newly filled missing values, while \n9.6% were changes between non-missing genotypes, and the remaining 0.7% changed from \nnon-missing to missing genotypes (Supplementary Figure 10). A disproportionately high \nnumber of non-missing genotype changes were within centromeric-like regions (the first 200 \nKb of chromosomes, corresponding to 0.2% of the genome), accounting for 11% of these \nunexpected changes. Other non-missing genotype changes occurred in regions where even \nmanual assignment of genotypes was not obvious from the alignments, with instances of \nforce-calling either improving or worsening genotype accuracy (Supplementary Figure 11).  \nGiven the near-complete catalogue of SVs established for the cohort, we further examined \nthe linkage disequilibrium between small variants and SVs. Even with the improved \nresolution of small variants within tandem repeats, we still found SVs overlapping tandem \nrepeats as the most poorly tagged class of SVs. Over half of SVs (52%) containing tandem \nrepeats were not in high linkage disequilibrium (r2>0.8) with any small variant within a ±1 Mb \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nwindow, compared to only 17% of SVs which did not contain tandem repeats or transposable \nelements. \n \nmolQTL mapping with HiFi data reveals the influence of structural variation on gene \nexpression and splicing \nThe contribution of SVs detected from long read sequencing to complex trait variability has \nnot been investigated in large cohorts of cattle. To do so, we assessed impacts of small \nvariants and SVs on molecular phenotypes through cis-molQTL mapping using deeply \nsequenced total RNA from testis tissue of the same bulls that was available from previous \nwork33. After quality control on the RNA sequencing data, we considered 117 bulls for \nexpression QTL (eQTL) and splicing QTL (sQTL) identification. Molecular phenotypes were \nestimated for 24,281 expressed genes and 16,975 spliced genes (corresponding to 50,853 \nintron clusters and 207,773 splice junctions), enabling association testing with 20,308,853 \nsmall variants (SNPs and small INDELs) and 61,861 SVs that were within 1 Mb of a \nmolecular feature's start site and had a MAF of at least 1% (Table 2). \nTable 2. Results for molecular QTL mapping with small variants and SVs.  \nVariant type \nTested \nvariants \neGenes eQTL \neVariants \n(unique) \nTop \nvariants \n(unique) \nsGenes sQTL \nsVariants \n(unique) \nTop \nvariants \n(unique) \nSmall \nvariants \n20,308,853 \n12,584 16,664 \n5,866,009 \n(3,612,237) \n16,539 \n(16,334) \n7,567 10,611 \n6,363,647 \n(2,357,818) \n10,400 \n(10,359) \nSVs 61,861 \n15,218 \n(9,618) \n105 \n(102) \n14,351 \n(5,535) \n211* \n(176) \n* 136 SV sQTL contained an SV with the same p-value as the top small variant \n \nHalf of the expressed genes had at least one eQTL (12,584 genes; eGenes), corresponding \nto 16,644 independent-acting eQTL and 3,621,855 variants that passed the nominal \nsignificant threshold (eVariants; Table 2). This included 261 eQTL (218 eGenes) on the X \nchromosome and 15 eQTL (13 eGenes) on the Y chromosome. Approximately 40% of \neGenes had at least one eVariant that was an SV. SVs were lead variants for 105 eQTL \n(hereafter referred to as SV eQTL; Table 2), i.e. they either had the smallest p-value or the \nsame p-value as the top small variant for an eGene. Three SVs were lead variants for \nmultiple eGenes, including an 85,550 bp duplication on chromosome 8 (8:103486033–\n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\n103571583) that was associated with ORM1 and COL27A1 expression. This SV was \npreviously described as a likely causal variant for an increased expression of ORM1 in liver \ntissue23. SV eQTL were depleted in intergenic regions and enriched in exons, CpG islands, \npromoters, and enhancers (Figure 3a; Supplementary Table 2). We identified 23 SV eQTL \nthat were not in strong LD (r2 >= 0.8) with a small variant, including a deletion that was \nassociated with the expression of ATXN7L3B on chromosome 5. Notably, this 218 bp \ndeletion was seven orders of magnitude more significant than the next eVariant (SV nominal \np-value = 1.55 × 10-17; Figure 3b,c). \nNearly half of the alternatively spliced genes possessed sQTL (Table 2). We detected at \nleast one sQTL for 10,388 intron clusters (23,292 splice junctions) within 7,567 genes \n(sGenes), which totalled to 2,363,353 significant variants (sVariants) and 10,611 \nindependent-acting sQTL. The X chromosome had 80 sQTL for 77 genes while the Y \nchromosome had 9 sQTL for 7 genes. Approximately 9% of the SVs considered for \nassociation testing were sVariants, totalling to 5,535 SV sVariants. SVs were lead variants \n(hereafter referred to as SV sQTL; Table 2) for 211 independent-acting sQTL for 173 genes \nand included 176 unique variants. One SV was associated with multiple intron clusters within \nthe same gene, while 30 SVs were associated with multiple junctions from the same cluster. \nOnly one SVs was associated with multiple sGenes; a 783 bp deletion on chromosome 25 \n(17,283,696 bp) was an SV sQTL for VPS35L, where the SV was located within an intron, \nand an sQTL for KNOP, which was approximately 23 Kb upstream of the variant. SV sQTL \nwere enriched in introns, exons, and splice sites, but depleted in intergenic regions \n(Supplementary Table 2). Approximately 16% of SV sQTL were not in strong LD with a small \nvariant (r2 >= 0.8), though some of these SVs were much more significant than the next \nsVariant, such as a 514 bp deletion on chromosome 23 that was associated with alternative \nsplicing of ADGRF1 (Figure 3d,e). One SV sQTL was not tagged (r2 < 0.2) by any nearby \nsmall variants—a 120 bp insertion on chromosome 21 which was the only variant that \npassed the significance threshold for LOC112443211. The significant SV and gene both \nresided within a complex region near the beginning of the chromosome with over 150-fold \ncoverage (at least 10 times higher than average) (Supplemental Figure 12). The lack of LD \nwith nearby SNPs is almost certainly driven by the SV and gene residing in a poorly resolved \nregion of the reference genome, although the sQTL may still reflect a valid association. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\n \nFigure 3. SV molQTL characteristics. (a) Enrichment of SV molQTL in functional elements. Odds ratios and p-\nvalues were inferred with a Fisher’s test, with significant (p < 0.01) observations circled in red. (b) Manhattan plot \nfor a poorly tagged SV eQTL that was associated with expression of ATXN7L3B. (c) Boxplot of ATXN7L3 \nexpression (TPM) across genotypes of the SV eQTL. Median TPM values for each genotype are reported above, \nand number of samples belonging to each genotype are reported below. (d) Manhattan plot for a poorly tagged \nSV sQTL that was associated with splicing of ADGRF1. (e) Boxplot of ADGRF1 intron usage (PSI) across \ngenotypes of the SV sQTL. Median PSI values for each genotype are reported above, and number of samples \nbelonging to each genotype are reported below. (f/g) MAF of each SV eQTL (f) and sQTL (g), and effect size \nmagnitude, coloured by variant type. (h) molQTL are depleted for rare variants, particularly for SVs. Colours \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\ncorrespond to variant type (SV or small variant) and whether the variant was a top variant. (i/j) Effect size \nmagnitude for eQTL (i) and sQTL (j) across variant types. Left panel shows all SV (>= 50 bp) and small variant (< \n50 bp) molQTL, while the right panel is separated by specific variant types. Median effect size (Beta) and \nstandard deviation (SD) is reported above, and number of variants within each category is below within \nparentheses. (k) Number of SV molQTL with TEs. Proportion is coloured by number of different TE classes for \neach SV eQTL or sQTL. (l) Proportion of TE SV molQTL across repeat classes.    \n \nWe found that SV eQTL and SV sQTL exhibited similar characteristics. Specifically, low \nfrequency variants had larger effect sizes (Figure 3f,g), as did variants that were within the \nfeature’s boundary (gene for eQTL, intron cluster and gene for sQTL; Supplemental Figure \n13). Both eQTL and sQTL were underrepresented among low-frequency variants, reflecting \nthe reduced statistical power to detect rare molecular QTLs (Figure 3h). This \nunderrepresentation was stronger for SVs than small variants. We observed no statistical \ndifference in effect size between SVs and small variants for eQTL (p = 0.23; Figure 3j), but \ninsertions and deletions had slightly smaller effect sizes than small variants for sQTL (p = 2.7 \n× 10-8; Figure 3j). Regardless, top variants for eQTL and sQTL were 2.1 (p = 9.29 × 10-11) \nand 5.6-fold (p = 1.92 × 10-70) enriched for SVs, respectively, highlighting the functional \nrelevance of SVs. Manual inspection of molQTL regions and HiFi alignments confirmed that \nthe SV top variants were genuine compelling causal candidates. \nApproximately 55% of the SV molQTL contained the full or partial sequence of at least one \ntransposable element (Figure 3k). LINEs were the most common transposable element class \namong SV molQTL; however, they were underrepresented compared to their overall \ngenome-wide abundance (Figure 3l; Supplementary Table 3). LTRs were also depleted \nwhen considering their overall abundance, while DNA and RNA transposable elements were \nslightly enriched among SV molQTL, though not significantly (Supplementary Table 3). There \nwas no statistical difference in effect size across transposable element classes for both SV \nmolQTL; though, when considering classes with more than two observations, median effect \nsize was slightly larger for RNA and LTR SV eQTL and LTR and SINE SV sQTL. Only 8% \nand 12% of SV eQTL and SV sQTL contained tandemly repeated motifs, which was lower \nthan the overall proportion of SVs containing such repetitive sequence. However, we \nencountered cases where incorrectly merged multi-allelic tandem repeats were omitted as \npotential SV molQTL (Supplemental Figure 14). \n \nImperfect genotyping led to an underestimation of SV molQTL \nInspection of the molQTL summary statistics revealed many SVs that appeared as strong \nputative candidate causal variants but had slightly larger p-values than the top small variant. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nThere were 167 eQTL and 168 sQTL where a small top variant was in perfect LD (r2 = 1.0) \nwith at least one SV that had a larger p-value (Figure 4a,c). This difference in p-values \nbetween two variants in perfect LD was often a result of an unbalanced proportion of missing \ngenotypes. Nearly all (97%) of these 335 molQTL had at least one “missing” call for the most \nsignificantly associated small variant, while 12 of the 289 SVs in LD had at least one \n“missing” call. Only 17 of the 335 molQTL had a “missing” call for both the small variant and \nthe SV in LD.  \n \nFigure 4. Misgenotyping of SVs impacts SV molQTL identification (a/c) Nominal p-value of top small variant and \nsignificant SVs for large effect eQTL (a) and sQTL (c). Colours correspond to whether the SV was in perfect LD \nwith the top small variant, had evidence of misgenotyping or if we did not identify an obvious error.  (b/d) Nominal \np-value of top small variant and significant SVs for large effect eQTL (b) and sQTL (d). Colours correspond to \nwhether the SV was in perfect LD with the top small variant, had evidence of misgenotyping or if we did not \nidentify an obvious error. (e) Comparing sample-wise average coverage across chromosome 2 and coverage of \nthe 681,722 bp duplication clearly separates heterozygous carriers of the duplication (black symbols) from non-\ncarriers (grey symbols). The red dotted line is an identity line. Four heterozygous samples (yellow stars) were \nmisgenotyped as homozygous reference by sniffles2. (f) Boxplots representing the expression of STK39 and \nCERS6 in carriers (0/1) and non-carries (0/0) of the 681,722 bp duplication the was misgenotyped for 4 samples \n(red symbols). Median TPM values for each group are listed above. \n \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nThe enrichment of SVs among molQTL was substantially stronger (6.1-fold enrichment for \neQTL, p = 9.38 × 10-130; 12.0-fold enrichment for sQTL, p = 8.04 × 10-256) when considering \nSVs in perfect LD with a top small variant. Furthermore, the proportion of SV eQTL and \nsQTL increased to 1.4% and 3.6%, respectively, when considering these molQTL. This \nincluded an MAGI3 eQTL for which a SNP (3:29677355_A_G) was a lead variant. This SNP \nhad a “missing” genotype in one sample and was in perfect LD with a 6,966 bp insertion that \nwas within an intron of MAGI3 but had non-missing genotypes for all samples. The nominal \np-value (p-value = 1.03 × 10-14) of the SNP was slightly lower than that of the SV (nominal p-\nvalue = 1.14 × 10-14; Supplemental Figure 15), demonstrating that an unequal number of \nsamples with missing genotypes across variants introduces bias in association tests. We \nalso identified instances where “missing” genotypes for SVs resulted in missed SV molQTL. \nFor instance, we observed an 8,790 bp insertion that had “missing” genotypes in four \nsamples, which we identified as homozygous for the reference allele upon manual inspection \nof the overlapping read alignments (Supplemental Figure 16). This insertion was in complete \nLD with a small variant sQTL for LOC768028. Exclusion of the “missed” genotype samples \nfor this SV during association testing had a substantial impact on the estimated p-value \n(nominal p-value = 4.97 × 10-12), which was consequentially fifteen orders of magnitude \nlarger than the p-value of the top SNP (nominal p-value = 3.60 × 10-27). \nThe effect size distributions for both sQTL and eQTL revealed that the most extreme \nmolecular QTL effects were primarily driven by small variants rather than SVs (Figure 3h,i). \nWe manually examined the top 1% of large effect molQTL with at least one SV passing the \nfeature’s significance threshold to determine whether erroneously genotyped SVs prevented \ntheir identification as molQTL. Among the 33 eQTL and 31 sQTL that passed these criteria, \nwe identified genotyping errors for SVs that were strong causal candidates for 9 eQTL and \n10 sQTL (27.3% and 32.3% of the manually examined large effect eQTL and sQTL, \nrespectively; Figure 4a,b,c,d). Most genotyping errors stemmed from low overall coverage, \nreads failing to span the length of the variant, and insufficient coverage of one haplotype \n(often leading to heterozygous genotypes being miscalled as homozygous; Supplemental \nFigure 17). Large insertions were particularly susceptible to misgenotyping, and we \nobserved multiple instances of insertions exceeding 8 Kb being missed as SV molQTL. \nNotably, the set of eVariants and sVariants were depleted for insertions between 8–9Kb, \nhighlighting the broader effects of misgenotyping variants of this size (eVariants: OR = 0.41, \np = 5.82 × 10-18; sVariants: OR = 0.27, p = 1.08 × 10-18). For example, we identified an 8,573 \nbp insertion that contained an L1 LINE element and was associated with an alternative \nsplicing event of PPM1H. This variant was nearly ten orders of magnitude less significant \nthan the top SNP (SV nominal p-value = 5.20 × 10-24; SNP nominal p-value = 5.13 × 10-33; \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nFigure 4c,d), despite being located within the intron that was affected by differential splicing. \nInspection of the HiFi alignments revealed that two samples were misgenotyped as \nhomozygous insertion carriers due to low coverage of the reference haplotype \n(Supplemental Figure 18). Correcting the two misgenotyped samples resulted in a stronger \nassociation between the splicing event and the SV (corrected SV nominal p-value = 1.31 × \n10-26). \nOur SV catalogue contained 43 duplications spanning from 3 to 956 Kb that completely \noverlapped 135 testis-expressed genes, yet only four were identified as molQTL. This low \nproportion is unexpected, considering that previous studies have highlighted the substantial \nimpact of such duplications on gene expression23. Nineteen of these 43 duplications could \nbe confidently verified though coverage-based genotyping, corroborating that sniffles can \nidentify duplications that are much larger than the HiFi reads. Taking the coverage-based \ngenotypes as truth and the sniffles-called genotypes as query, we found that at least one \nsample was misgenotyped for all nineteen variants examined. Overall, we found that a \nsubstantial fraction of the sniffles-called genotypes for these long duplications were \nerroneous in both lowly and highly covered samples (Supplemental Figure 19). We identified \na relatively common 682 Kb duplication on chromosome 2 (2:27313972–27995694), which \noverlaps the entire coding sequences of CERS6 and STK39. This duplication was called as \nheterozygous by sniffles in 16 samples; however, the coverage-based genotyping indicated \nthat 4 heterozygous samples were incorrectly genotyped as homozygous reference (i.e., did \nnot have the duplication; Figure 4e; Supplemental Figure 20). These genotyping errors \nresulted in this duplication being missed as an SV eQTL for both overlapped genes. We \nobserved that the top small variants for these genes were one to three orders of magnitude \nmore significantly associated than the erroneously genotyped duplication (top small variant \np-value CERS6 = 3.31 × 10-15, STK39 = 1.90 × 10-14; duplication p-value CERS6 = 8.34 × \n10-12, STK39 = 7.17 × 10-13; Figure 4f). However, if genotyped correctly, the duplication is the \ntop variant for both genes (corrected duplication p-value CERS6 = 1.33 × 10-15, STK39 = \n3.96 × 10-15). \n \n \nDiscussion \nWe have generated and analysed a cohort of Braunvieh cattle that exceeded the effective \npopulation size with long and accurate HiFi reads, representing one of the largest long-read \ndatasets available for a non-model organism. Notably, this is the first bovine long-read DNA \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\ncohort which also has paired functional data, providing the opportunity to investigate the \nconsequences of small and structural variation on molecular phenotypes. This cohort \nenabled us to detect and characterize 23,789,890 small variants and 79,275 SVs, including \ntens of thousands of previously unidentified SVs. Our collection of SVs is almost twice as \nlarge as that identified in a previously built 16-sample cattle pangenome35 and substantially \nlarger than previously reported in similar-sized, or larger, short-read sequenced cattle \ncohorts22–24. While short-read based SV cohorts are either depleted for insertions or lack \nthem completely10,22,36, we identified—as expected19,37—more insertions than deletions from \nthe long read alignments, thereby providing comprehensive access to a vastly understudied \nvariant type. A saturation analysis corroborated that our SV set is a near complete catalogue \nof common SVs segregating within the Braunvieh population. However, we suspect that \nlonger SVs (>15,000 bp) were missed for samples with shorter average read lengths. Long-\nread sequence data also improved alignments to the X and Y chromosomes—which are \noften neglected in genome-wide association studies due to their repetitiveness—and thereby \nfacilitated the identification of over one hundred thousand small variants that were previously \ninaccessible with short-read sequencing, including 7 SNPs predicted as HIGH impact. We \nalso demonstrated improved small variant calling on the autosomes, particularly in tandem \nrepeat regions, when using long reads over short reads, thereby creating a comprehensive \ngenetic resource for genome-wide association testing. \nOur expression and splicing QTL analyses revealed that molQTL are enriched for SVs, \nemphasizing their crucial role in complex trait variation, as previously reported in cattle and \nother species14,38–40. We observed that 0.6% of eQTL and 1.9% of sQTL had an SV as the \ntop variant. These proportions of SV molQTL were between 2.1-times and 14-times greater \nthan those reported in short-read-based studies in bulls22,35. More than half of the SV \nmolQTL were due to insertions that were largely inaccessible with short reads. Notably, the \nnumber of the SV molQTL nearly doubled when accounting for SVs in complete LD with a \nlead small variant, increasing to 1.4% of eQTL and 3.6% of sQTL. These estimates are more \ncomparable to short- and long-read based SV molQTL studies that were conducted in \nhumans with higher coverage sequence data and larger sample sizes40,41. Contrary to other \nSV molQTL studies, we did not observe a statistical difference in effect size magnitude \nbetween small variants and SVs40. However, we identified multiple large effect molQTL \nwhere small variants were prioritized as candidate causal variants despite the presence of \ncompelling candidate causal SVs. We suspect that several large allelic substitution effects \nwere erroneously attributed to the small variants due to imperfect SV and small variant \ngenotyping, which introduces some bias in this assessment.   \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nThe forced genotyping approach substantially reduced the proportion of missing genotypes \nfor all types of SVs. Yet, many SVs were missed as molQTL due to erroneous SV genotypes \nand an unequal distribution of missing genotypes between SVs and small variants. Almost \none-third of the manually examined molQTL and all manually examined large duplications \nwere affected by imperfectly genotyped small variants and SVs. Using the manually revised \ngenotypes, we identified compelling candidate causal SVs in many of the examined molQTL \nwhere the summary statistics from the association tests initially prioritized small variants. We \nwere able to attribute some of the genotyping errors to insufficient coverage. This suggests \nthat an average HiFi coverage higher than the 13.5-fold obtained for our cohort will produce \nmore accurate genotypes for both small variants and SVs which will benefit association \ntesting. Such findings also indicate that SVs are even stronger enriched among molQTLs \nthan reflected by the results from the statistical association testing. Imputation can enhance \nstatistical power and mitigate some of the missing genotype biases in association testing \nthat we observed42,43. However, the accuracy of imputation can be much lower for SVs than \nsmall variants, which negatively impacts association testing and causal variant \nidentification44. Since it remains unclear whether the benefits outweigh the drawbacks, we \ndid not attempt to impute missing genotypes neither for SVs nor for small variants. \nGenerating high coverage long-read data at the sample-sizes required for detecting QTL is \nbecoming increasingly feasible16,45, which will reduce the missingness, thereby improving the \nstatistical power of association studies utilising long-reads for variant discovery and \ngenotyping.  \nGenotyping errors within complex regions, such as those containing multiallelic SVs or \ntandem repeats, and an inappropriate merging of SVs within these regions likely caused the \nobserved underrepresentation of SV molQTL annotated as tandem repeats. In addition, the \nlack of reads spanning large variants resulted in erroneously genotyped large SVs, including \nsome molecular phenotype-associated variants, which remained undetected by our \nassociation analyses. This contributed to the lack of association of gene-containing \nduplications and likely contributed to the observed underrepresentation of LINE and LTR \ntransposable element molQTL, which have previously been associated with large impacts on \ngene expression46,47. Further research is warranted to develop refined genotyping and \nassociation testing methods to fully capitalize on long-read sequencing data and explore the \nassociation between functionally relevant SVs and complex traits 48. \n \n \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nMethods \nEthics statement \nThe Bos taurus taurus tissues used in this study were previously collected from a \ncommercial abattoir in Zürich, Switzerland after regular slaughter. None of the authors were \ninvolved in the decision to slaughter the bulls. No ethics approval was required for this study. \n \nSample extraction and DNA sequencing \nReproductive tissues (including testis and the caput of the epididymis) from 120 mature bulls \nwere sampled by Mapel et al33. Tissue samples were flash frozen in liquid nitrogen and \nstored at -80 C. Sampling was conducted between 2019 and 2021 and HiFi sequencing was \nperformed from December 2022 to May 2023.  \nWe extracted HMW DNA from flash-frozen testis (N=105) or epididymis (N=15) tissue \nsamples with the Monarch HMW DNA Extraction Kit for Tissue (New England Biolabs). We \nfollowed the manufacturer recommended protocol for both tissue types. DNA was shipped \non dry ice to PacBio (Rolling Stock Yard, London) for fragment length analysis, library \npreparation, and sequencing. \n \nDNA read alignment \nWe aligned the HiFi reads with minimap2 (v2.24)49 to ARS-UCD2.0 with the X chromosomal \npseudo-autosomal region (PAR) hard masked (X:133300518-139009144), using the “map-\nhifi” preset. Alignments were coordinate sorted with SAMtools (v1.19.2)50. The previously \ncollected short reads were aligned with bwa-mem2 (v2.2.1)51 to the same reference, before \ncollating by name, marking duplicates, and coordinate sorting with SAMtools. \nWe assessed alignment coverage using SAMtools bedcov and BEDtools (v2.30)52 coverage, \nusing 100 Kb windows generated by BEDtools makewindows. Regions with a minimum \nprimary read depth of 2 and MAPQ of 5 were considered suitable for variant calling. \n \nSmall variant calling \nWe called small variants with DeepVariant (v1.6.0)53, using the “PACBIO” model for the HiFi \nreads and “WGS” for the short Illumina paired-end reads, specifying --haploid_contigs \"X,Y” \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nand the PAR bed file, producing a gvcf for each sample. Samples were merged into \npopulation vcfs using GLnexus (v1.4.1)53 with the DeepVariant config and \nrevise_genotypes=False. DeepVariant variants larger than 45 bp were removed due to avoid \nunintentional overlapping with the dedicated structural variant calling. \nWe classified regions as tandem repeats using the pbsv (v2.9.0) utility findTandemRepeats \n(https://github.com/PacificBiosciences/pbsv/blob/master/annotations/findTandemRepeats) or \nas transposable elements using the RepeatMasker (v.4.1.5) (https://www.repeatmasker.org/) \nutility RM2Bed.py on the USCS browser repeat file \n(https://hgdownload.soe.ucsc.edu/hubs/GCF/002/263/795/GCF_002263795.3/GCF_002263\n795.3.repeatMasker.out.gz) and keeping all SINE/LINE/LTR elements. Regions overlapping \nin tandem repeats and transposable elements were prioritized into the tandem repeat bed. \nAll remaining regions were assigned as “normal” using BEDtools complement. \nWe assessed variant calling and genotyping accuracy per sample with hap.py (v0.3.15), \nremoving all reference calls, left-shifting indels, and stratifying with the genomic regions \nspecified above. \n \nStructural variant calling \nWe called SVs with sniffles2 (v2.2)54, providing the tandem repeat locations from above. We \nmerged per-sample snf files across the cohort with sniffles2, removing any breakend SVs or \nSVs larger than 1 Mb, before force-genotyping each sample again with the merged set of SV \ncandidates. We finally merged force-called VCFs with BCFtools merge and removed SVs \nwith more than 10% missingness. \nWe assessed the repetitive element and tandem repeat content of the SV sequence (the \nreference for deletions and alternate for insertions) using Repeatmasker and TRF \n(v4.09.1)55, respectively. \n \nVariant analyses \nWe used VEP (release 113)56 to assess functional impacts and consequences of small \nvariants using the RefSeq (release 106) cattle annotation. We used plink (v1.90b6.26)57 to \ncalculate the LD between SVs and small variants, using a 1 Mb window and a minimum r2-\nthreshold of 0.2 and 0.8 to assess weak and strong tagging, respectively. We tested for \nmutually present SVs between our cohort and the SV pangenome panel from Leonard et \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nal.35 using Jasmine (v1.1.5)58, with the flags “--pre_normalize max_dist_linear=0.5 \nmax_dist=250” to allow slight mismatches in SV position or length between the two callsets. \nWe used BEDtools intersect to examine the overlap between SVs and the previously \ndescribed tandem repeat file. \n \nCoverage analysis \nWe used mosdepth (version 0.2.2)59 to calculate the average coverage in 500 bp windows \nfor the HiFi alignments. Heterozygous and homozygous genotypes were assigned to the \nsniffles2-called duplications when the average coverage over the duplicated sequence was \nat least 1.3-fold and 1.8-fold higher, respectively, than the average coverage outside the \nduplication on the chromosome containing the duplication. \n \nMolecular phenotype preparation \nTotal RNA from the HiFi-sequenced individuals was sequenced previously and made publicly \navailable in Mapel et al.33 (Supplementary Table 4). We filtered RNA reads with fastp \n(v0.23.1)60 to remove adapter sequences, poly-A-tails, ploy-G-tails, and low-quality bases.  \nWe used Kallisto61 (v0.50.0) and the RefSeq annotation for cattle to quantify transcript \nexpression (TPM) and counts, which were aggregated to gene-level with tximport61 \n(v1.34.0). We considered genes with TPM >= 0.1 and at least 6 supporting reads in >=10% \nof samples for eQTL mapping. Expression values were inverse normal transformed and \nquantile normalized.  \nTo identify and quantify splicing events, we first aligned cleaned reads to the ARS-UCD2.0 \nreference and RefSeq annotation with STAR (version 2.7.11a)62 and included WASP filtering \nto account for allelic bias, using heterozygous sites called from the HiFi reads63. We \nextracted exon-exon junctions with Regtools (v0.5.2)64, then used LeafCutter (v0.2.9)65 to \ncluster introns, calculate intron excision ratios, perform filtering, and normalize splicing \nphenotypes for sQTL mapping.  \n \nmolQTL mapping \nCis-molQTL mapping was conducted with QTLtools (v1.3.1)66. We included variants that \nwere within a 1 Mb window of the feature’s start site and had MAF >= 1%. We selected \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\ncovariates for eQTL and sQTL testing with PCAForQTL67 which implemented an elbow test \nto estimate the number of principal components (PCs) from the expression or splicing matrix \nto include as hidden confounders (hereafter, “RNA PCs”). We used the \nfilterKnownCovariates function (R2 thresholds = 0.5) to remove known covariates that were \ncaptured by the RNA PCs. The covariates considered for eQTL mapping included PCs 2-10 \nof the genome relationship matrix (constructed from 381,111 unlinked variants with plink \n(v2.00a3.6LM)57 and --indep-pairwise 1000 5 0.2, RIN, age, and RNA PCs 1-12. The sQTL \ncovariates included PCs 2-10 of the genome relationship matrix, RIN, age, and RNA PCs 1-\n11. We conducted 1,000 permutations with --permute function in QTLtools to infer beta \ncorrected p-values, then used the qtltools_runFDR_cis.R script to apply a 5% false discovery \nrate (FDR) and estimate significance thresholds for each gene. sQTL testing include the --\ngrp-best flag to account for multiple intron clusters within a gene. We performed a \nconditional analysis (as described in Delaneau et al.66) to identify independent-acting signals \nfor each gene or intron cluster and the corresponding most significant variant. \nWe defined “SV molQTL” as QTL for which either an SV had the smallest p-value, or an SV \nhad a p-value that was identical with a lead small variant. Enrichment was inferred with a \nFisher’s exact test. We used BEDtools to identify SVs that overlapped functional elements. \nPromoters and enhancer annotations for cattle testis tissue were obtained from Salavati et \nal.68, while CpG islands were obtained from the UCSC genome browser69. \n \nData availability  \nDNA and RNA sequencing data of the analysed cohort are available in the ENA database at \nthe study accessions PRJEB42335 (Long-read sequencing data from cattle for the purpose \nof de-novo genome assembly), PRJEB28191 (Short read sequencing of cattle) and \nPRJEB46995 (Testis transcriptome of mature bulls). Accession identifiers for all samples are \navailable as Supplementary Table 4. Gene expression and splicing matrices, a VCF file of \ngenome-wide small and structural variant genotypes used for e/sQTL mapping, a cross-table \nto link genotype and transcriptome data as well as results from e/sQTL mapping have been \narchived at zenodo (https://zenodo.org/records/15431127). \n \nCode availability  \nAll workflows are available through https://github.com/AnimalGenomicsETH/HiFi_cohort  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\nAcknowledgements \nThis study was supported by an ETH Research Grant and a grant from the Swiss National \nScience Foundation (SNSF, grant-ID 204654). The funding bodies were neither involved in \nthe design of the study and collection, analysis, and interpretation of data nor in writing the \nmanuscript. We thank Eirini Lampraki from Pacific Biosciences for DNA fragment analysis \nand sequencing. We thank Audald Lloret-Villas and Qiongyu He for valuable discussions. \n \nAuthor contributions \nX.M.M. sampled tissue and purified HMW DNA, aligned RNA reads against the reference, \ndeveloped and applied workflows to quantify gene expression and splicing variation, \nconducted molecular QTL mapping, interpreted results, and drafted the manuscript; A.S.L. \naligned DNA reads against the reference, called variants from short and long read \nalignments, interpreted results, and drafted the manuscript; H.P. conceived the study, \ncontributed to the analysis of expression and splicing QTL, interpreted results, and \ncontributed to the writing of the manuscript. All authors approved the final version of the \nmanuscript. \n  \nCompeting interests \nThe authors declare no competing interests. \n \nReferences \n1. Daetwyler, H. D. et al. Whole-genome sequencing of 234 bulls facilitates mapping of \nmonogenic and complex traits in cattle. Nat Genet 46, 858–865 (2014). \n2. Altshuler, D. M. et al. An integrated map of genetic variation from 1,092 human \ngenomes. Nature 491, 56–65 (2012). \n3. Sollis, E. et al. The NHGRI-EBI GWAS Catalog: knowledgebase and deposition \nresource. Nucleic Acids Res 51, D977–D985 (2023). \n4. Hu, Z. L., Park, C. A. & Reecy, J. M. Bringing the Animal QTLdb and CorrDB into the \nfuture: Meeting new challenges and providing updated services. Nucleic Acids Res \n50, D956–D961 (2022). \n5. Canela-Xandri, O., Rawlik, K. & Tenesa, A. An atlas of genetic associations in UK \nBiobank. Nat Genet 50, 1593–1599 (2018). \n6. Weischenfeldt, J., Symmons, O., Spitz, F. & Korbel, J. O. Phenotypic impact of \ngenomic structural variation: Insights from and for human disease. Nature Reviews \nGenetics vol. 14 125–138 Preprint at https://doi.org/10.1038/nrg3373 (2013). \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\n7. Sedlazeck, F. J. et al. Accurate detection of complex structural variations using single-\nmolecule sequencing. Nat Methods 15, 461–468 (2018). \n8. Huddleston, J. et al. Discovery and genotyping of structural variation from long-read \nhaploid genome sequence data. Genome Res 27, 677–685 (2017). \n9. Kosugi, S. & Terao, C. Comparative evaluation of SNVs, indels, and structural \nvariations detected with short- and long-read sequencing data. Hum Genome Var 11, \n(2024). \n10. Cameron, D. L., Di Stefano, L. & Papenfuss, A. T. Comprehensive evaluation and \ncharacterisation of short read general-purpose structural variant calling software. Nat \nCommun 10, (2019). \n11. Ahsan, M. U., Liu, Q., Perdomo, J. E., Fang, L. & Wang, K. A survey of algorithms for \nthe detection of genomic structural variants from long-read sequencing data. Nature \nMethods vol. 20 1143–1158 Preprint at https://doi.org/10.1038/s41592-023-01932-w \n(2023). \n12. Mikheyev, A. S. & Tin, M. M. Y. A first look at the Oxford Nanopore MinION sequencer. \nMol Ecol Resour 14, 1097–1102 (2014). \n13. Eid, J. et al. Real-Time DNA Sequencing from Single Polymerase Molecules. Science \n(1979) 323, 133–138 (2009). \n14. Alonge, M. et al. Major Impacts of Widespread Structural Variation on Gene \nExpression and Crop Improvement in Tomato. Cell 182, 145-161.e23 (2020). \n15. Beyter, D. et al. Long-read sequencing of 3,622 Icelanders provides insight into the \nrole of structural variants in human diseases and other traits. Nat Genet 53, 779–786 \n(2021). \n16. Schloissnig, S. et al. Long-read sequencing and structural variant characterization in \n1,019 samples from the 1000 Genomes Project. bioRxiv (2024) \ndoi:10.1101/2024.04.18.590093. \n17. Duan, X., Pan, M. & Fan, S. Comprehensive evaluation of structural variant \ngenotyping methods based on long-read sequencing data. BMC Genomics 23, \n(2022). \n18. Yang, Q. et al. SVLearn: a dual-reference machine learning approach enables \naccurate cross-species genotyping of structural variants. Nature Communications 16, \n(2025). \n19. Leonard, A. S., Crysnanto, D., Mapel, X. M., Bhati, M. & Pausch, H. Graph \nconstruction method impacts variation representation and analyses in a bovine super-\npangenome. Genome Biol 24, (2023). \n20. Talenti, A. et al. A cattle graph genome incorporating global breed diversity. Nat \nCommun 13, (2022). \n21. Leonard, A. S. et al. Structural variant-based pangenome construction has low \nsensitivity to variability of haplotype-resolved bovine assemblies. Nat Commun 13, \n(2022). \n22. Bhati, M., Mapel, X. M., Lloret-Villas, A. & Pausch, H. Structural variants and short \ntandem repeats impact gene expression and splicing in bovine testis tissue. Genetics \n225, (2023). \n23. Lee, Y. L. et al. High-resolution structural variants catalogue in a large-scale whole \ngenome sequenced bovine family cohort data. BMC Genomics 24, (2023). \n24. Grant, J. R. et al. A large structural variant collection in Holstein cattle and associated \ndatabase for variant discovery, characterization, and application. BMC Genomics 25, \n903 (2024). \n25. Lee, Y. L. et al. A 12 kb multi-allelic copy number variation encompassing a GC gene \nenhancer is associated with mastitis resistance in dairy cattle. PLoS Genet 17, (2021). \n26. Trigo, B. B. et al. Variants at the ASIP locus contribute to coat color darkening in \nNellore cattle. Genetics Selection Evolution 53, (2021). \n27. Rothammer, S. et al. The 80-kb DNA duplication on BTA1 is the only remaining \ncandidate mutation for the polled phenotype of Friesian origin. Genetics Selection \nEvolution 46, (2014). \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\n28. Milia, S. et al. Taurine pangenome uncovers a segmental duplication upstream of KIT \nassociated with depigmentation in white-headed cattle. Genome Res (2024) \ndoi:10.1101/gr.279064.124. \n29. Durkin, K. et al. Serial translocation by means of circular intermediates underlies \ncolour sidedness in cattle. Nature 482, 81–84 (2012). \n30. Küttel, L. et al. A complex structural variant at the KIT locus in cattle with the \nPinzgauer spotting pattern. Anim Genet 50, 423–429 (2019). \n31. Kadri, N. K. et al. A 660-Kb Deletion with Antagonistic Effects on Fertility and Milk \nProduction Segregates at High Frequency in Nordic Red Cattle: Additional Evidence \nfor the Common Occurrence of Balancing Selection in Livestock. PLoS Genet 10, \n(2014). \n32. Venhoranta, H. et al. Ectopic KIT Copy Number Variation Underlies Impaired \nMigration of Primordial Germ Cells Associated with Gonadal Hypoplasia in Cattle (Bos \ntaurus). PLoS One 8, (2013). \n33. Mapel, X. M. et al. Molecular quantitative trait loci in reproductive tissues impact male \nfertility in cattle. Nat Commun 15, (2024). \n34. Olagunju, T. A. et al. Telomere-to-telomere assemblies of cattle and sheep Y-\nchromosomes uncover divergent structure and gene content. Nat Commun 15, 8277 \n(2024). \n35. Leonard, A. S., Mapel, X. M. & Pausch, H. Pangenome-genotyped structural variation \nimproves molecular phenotype mapping in cattle. Genome Res 34, 300–309 (2024). \n36. Sudmant, P. H. et al. An integrated map of structural variation in 2,504 human \ngenomes. Nature 526, 75–81 (2015). \n37. Ebert, P. et al. Haplotype-resolved diverse human genomes and integrated analysis of \nstructural variation. Science (1979) 372, (2021). \n38. Scott, A. J., Chiang, C. & Hall, I. M. Structural variants are a major source of gene \nexpression differences in humans and often affect multiple nearby genes. Genome \nRes 31, 2249–2258 (2021). \n39. Zhang, Y. et al. Structural variation reshapes population gene expression and trait \nvariation in 2,105 Brassica napus accessions. Nat Genet (2024) doi:10.1038/s41588-\n024-01957-7. \n40. Chiang, C. et al. The impact of structural variation on human gene expression. Nat \nGenet 49, 692–699 (2017). \n41. Billingsley, K. J. et al. Long-read sequencing of hundreds of diverse brains provides \ninsight into the impact of structural variation on gene expression and DNA \nmethylation. bioRxiv (2024) doi:10.1101/2024.12.16.628723. \n42. Huang, L., Wang, C. & Rosenberg, N. A. The Relationship between Imputation Error \nand Statistical Power in Genetic Association Studies in Diverse Populations. Am J \nHum Genet 85, 692–698 (2009). \n43. Marchini, J. & Howie, B. Genotype imputation for genome-wide association studies. \nNature Reviews Genetics vol. 11 499–511 Preprint at https://doi.org/10.1038/nrg2796 \n(2010). \n44. Noyvert, B. et al. Imputation of structural variants using a multi-ancestry long-read \nsequencing panel enables identification of disease associations. medRxiv (2023) \ndoi:10.1101/2023.12.20.23300308. \n45. Gong, J. et al. Long-read sequencing of 945 Han individuals identifies structural \nvariants associated with phenotypic diversity and disease susceptibility. Nat Commun \n16, 1494 (2025). \n46. Tang, L. et al. GWAS reveals determinants of mobilization rate and dynamics of an \nactive endogenous retrovirus of cattle. Nat Commun 15, (2024). \n47. Adelson, D. L., Raison, J. M. & Edgar, R. C. Characterization and Distribution of \nRetrotransposons and Simple Sequence Repeats in the Bovine Genome. PNAS \nAugust vol. 4 (2009). \n48. Cui, Y. et al. Multi-omic quantitative trait loci link tandem repeat size variation to gene \nregulation in human brain. Nat Genet (2025) doi:10.1038/s41588-024-02057-2. \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint \n\n49. Li, H. Minimap2: Pairwise alignment for nucleotide sequences. Bioinformatics 34, \n3094–3100 (2018). \n50. Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, \n2078–2079 (2009). \n51. Li, H. Aligning sequence reads, clone sequences and assembly contigs with BWA-\nMEM. (2013). \n52. Quinlan, A. R. & Hall, I. M. BEDTools: A flexible suite of utilities for comparing genomic \nfeatures. Bioinformatics 26, 841–842 (2010). \n53. Yun, T. et al. Accurate, scalable cohort variant calls using DeepVariant and GLnexus. \nBioinformatics 36, 5582–5589 (2020). \n54. Smolka, M. et al. Detection of mosaic and population-level structural variants with \nSniffles2. Nat Biotechnol (2024) doi:10.1038/s41587-023-02024-y. \n55. Benson, G. Tandem Repeats Finder: A Program to Analyze DNA Sequences. Nucleic \nAcids Research vol. 27 https://academic.oup.com/nar/article/27/2/573/1061099 \n(1999). \n56. McLaren, W. et al. The Ensembl Variant Effect Predictor. Genome Biol 17, (2016). \n57. Chang, C. C. et al. Second-generation PLINK: Rising to the challenge of larger and \nricher datasets. Gigascience 4, (2015). \n58. Kirsche, M. et al. Jasmine and Iris: population-scale structural variant comparison and \nanalysis. Nat Methods 20, 408–417 (2023). \n59. Pedersen, B. S. & Quinlan, A. R. Mosdepth: Quick coverage calculation for genomes \nand exomes. Bioinformatics 34, 867–868 (2018). \n60. Chen, S., Zhou, Y., Chen, Y. & Gu, J. Fastp: An ultra-fast all-in-one FASTQ \npreprocessor. in Bioinformatics vol. 34 i884–i890 (Oxford University Press, 2018). \n61. Soneson, C., Love, M. I. & Robinson, M. D. Differential analyses for RNA-seq: \ntranscript-level estimates improve gene-level inferences. F1000Res 4, 1521 (2015). \n62. Dobin, A. et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21 \n(2013). \n63. Van De Geijn, B., Mcvicker, G., Gilad, Y. & Pritchard, J. K. WASP: Allele-specific \nsoftware for robust molecular quantitative trait locus discovery. Nature Methods vol. \n12 1061–1063 Preprint at https://doi.org/10.1038/nmeth.3582 (2015). \n64. Feng, Y.-Y. et al. RegTools: Integrated analysis of genomic and transcriptomic data for \ndiscovery of splicing variants in cancer. bioRxiv (2018) doi:10.1101/436634. \n65. Li, Y. I. et al. Annotation-free quantification of RNA splicing using LeafCutter. Nat \nGenet 50, 151–158 (2018). \n66. Delaneau, O. et al. A complete tool set for molecular QTL discovery and analysis. Nat \nCommun 8, (2017). \n67. Zhou, H. J., Li, L., Li, Y., Li, W. & Li, J. J. PCA outperforms popular hidden variable \ninference methods for molecular QTL mapping. Genome Biol 23, (2022). \n68. Salavati, M. et al. Improving the annotation of the cattle genome by annotating \ntranscription start sites in a diverse set of tissues and populations using Cap Analysis \nGene Expression sequencing. G3: Genes, Genomes, Genetics 13, (2023). \n69. Karolchik, D. et al. The UCSC Genome Browser Database. Nucleic Acids Research \nvol. 31 51–54 Preprint at https://doi.org/10.1093/nar/gkg129 (2003). \n  \n.CC-BY 4.0 International licensemade available under a \n(which was not certified by peer review) is the author/funder, who has granted bioRxiv a license to display the preprint in perpetuity. It is \nThe copyright holder for this preprintthis version posted May 16, 2025. ; https://doi.org/10.1101/2025.05.16.654493doi: bioRxiv preprint","source_license":"CC-BY-4.0","license_restricted":false}