Graph construction method impacts variation representation and analyses in a bovine super-pangenome

We built pangenomes from autosomal sequences of domestic cattle and three of their wild relatives using minigraph, cactus, and pggb to assess how different graph pangenomes represent the same underlying input sequences. Nineteen haplotype-resolved assemblies from eight breeds of taurine (Bos taurus taurus) and indicine (Bos taurus indicus) cattle, yak (Bos grunniens), bison (Bison bison bison), and gaur (Bos gaurus) as well as the current Hereford-based Bos taurus taurus reference genome sequence were considered (Fig. 1a), of which the reference genome and eleven haplotype-resolved assemblies were used for pangenome construction and eight intermediate-quality assemblies were held out for analysis. These genomes were assembled using different sequencing and algorithmic approaches, but this has limited effect on the pangenomes [10]. However, the larger amount of centromeric/telomeric sequence in the HiFi-based assemblies (Fig. 1b, c) does require additional consideration. Minigraph used ARS-UCD1.2 as a backbone, whereas pggb and cactus are reference-free methods, although cactus is guided by an approximate phylogenetic tree (black subset of Fig. 1d).

Fig. 1figure 1

Input assemblies considered for the pangenome analyses. a Autosomal length of the Bos taurus taurus reference sequence (HER, ARS-UCD1.2) and 19 haplotype resolved assemblies considered during pangenome construction and downstream analyses. Three pangenomes were created with ARS-UCD1.2 [HER] and eleven input assemblies (Brown Swiss [BSW], Piedmontese [PIE], Highland [HIG], Angus [ANG], Original Braunvieh [OBV], Simmental [SIM], Brahman [BRA], Nellore [NEL], gaur [GAU], bison [BIS], yak [YAK]), whereas eight additional Original Braunvieh or Brown Swiss assemblies indicated with red text (OD1, OD2, OS1, OS2, B11, B12, B21, B22) were only considered for downstream analyses. The colour of the bars indicates the primary sequencing technology used to construct the assemblies. The black dashed line indicates the length of ARS-UCD1.2. b Centromeric and c telomeric completeness is generally higher in the HiFi- and ONT- than CLR-based input assemblies. The marker is the sum over the autosomes and error bars indicate the 95% confidence interval. d A tree constructed with mash from the assemblies reveals the expected separation between taurine, indicine and non-cattle. The asterisk indicates the backbone genome used by minigraph

Pangenome construction and sequence content

The three pangenomes spanned between 427 k and 198 M nodes and contained between 2.6 and 3.0 gigabases (Table 1). Minigraph primarily incorporates larger (> 50 bp) DNA variation, and so builds a smaller graph compared to pggb and cactus which include all sizes of variation, which is also reflected in the number of path steps in the graph and final file sizes. The non-reference nodes contain nearly five times as much sequence in cactus (552 Mb) and pggb (523 Mb) than minigraph (109 Mb). Pangenome construction took 16 and 253 times more CPU hours and required 8 and 7 times more memory for cactus and pggb than minigraph.

Table 1 Profiles of three bovine pangenomes constructed from autosomal sequences of twelve assemblies. Path steps are the total number of steps needed to trace all 12 assemblies through the graph

Pggb and cactus also contain more repetitive sequence than minigraph, 46.4%, 44.3%, and 42.2% respectively, although this is largely due to including centromeric sequence. The ARS-UCD1.2 reference genome contains little centromeric sequence, preventing minigraph from anchoring centromeric sequence present in the HiFi-based (Brown Swiss, Nellore, Original Braunvieh, Piedmontese, and gaur) and to a lesser extent in ONT-based (bison and Simmental) assemblies into its pangenomes. Although pggb and cactus do include centromeric sequence into their pangenomes, they span a similar amount of bases as the sum of the input assemblies, suggesting either they are biologically distinct or cannot be effectively collapsed into a graph representation in a way similar to other repeat elements (e.g. SINE, LINE, LTR, etc., Additional file 1: Fig. S1).

Consensus of genomic variation

To assess the commonality of genomic variation across pggb, cactus, and minigraph, as well as the suitability of their representation for downstream analyses, we decomposed and normalized the graph pangenomes into VCF files with respect to the ARS-UCD1.2 reference. Decomposing default minigraph pangenomes failed to output many expected structural variant (SV) alleles, particularly deletions in multi-node bubbles (Additional file 1: Fig. S2). By adding path information (P-lines in GFA) through minigraph-based realignment, we recovered an additional 21,770 SVs (13.5% increase). Realignment rarely produces paths incongruent with graph topology, where an assembly does not trace through nodes originating from that same assembly, affecting approximately 0.03%, 0.03%, and 0.05% of taurine, indicine, and non-cattle nodes (Additional file 1: Fig. S3). Cactus had the most SV genotypes marked as CONFLICT (2.2%), indicating that a haploid assembly had multiple possible ambiguous alleles, with minigraph and pggb significantly lower, 0.5% and 0.4% respectively. Such conflicts occurred more often in SV than small variation alleles, and were more frequent in divergent assemblies for cactus (Additional file 1: Fig. S4). VCF decomposition took 142 and 50 times longer for cactus and pggb compared to minigraph respectively, using approximately 26 and 22 times more memory (Additional file 2: Table S1), although pggb and cactus also contain substantial small variation to process.

Merging SVs, allowing for variation in breakpoints up to the minimum of 50% of the SV length and 1 Kb, identified 221 k nonredundant SVs ≥ 50 bp from the three pangenomes of which 135 k (61.2%) were common (Fig. 2a). Cactus contained 184 k SVs of which 26 k (14.1%) were private, i.e. not present in minigraph or pggb. We identified fewer SVs with pggb (172 k) and minigraph (169 k), of which 8 k (4.7%) and 9 k (5.0%) were respectively private. We classified the ARS-UCD1.2 reference genome into five regions: centromeric satellites (0.13%), tandem repeats (6.21%), repetitive elements (35.78%), low mappability (0.44%), and “normal” which encapsulated all remaining sequence (57.44%). The overlap between the three pangenomes was highest in normal (76.8%, 53.6 k) and repetitive regions (70.3%, 50.7 k), while tandem repeat regions (52.5%, 35.8 k) had lower overlap. The most challenging regions for alignment, low mappability (23.0%, 566) and especially centromeric satellites (2.1%, 81), had substantially lower overlap between the pangenomes. This also reaffirms that tandem repeats are disproportionately large contributors to structural variation [17]. Applying stricter requirements to merge SVs across the pangenomes, requiring variant breakpoints to be within a 5 bp window, had limited effect, and even requiring exact basepair resolution had 92 k overlapping SVs, suggesting most SVs are precisely represented in the pangenomes (Additional file 1: Fig. S5).

Fig. 2figure 2

Consensus of structural variation between three pangenomes. a Overlap of SVs between three pangenomes and assembly-based SV discovery. Overall bar height represents the total number of SVs. Variants in different classifications of genomic regions are indicated by colour of the stacked bar. b Number of SVs identified in the input assemblies through assembly-based mapping and the three pangenomes. c Allele count of the SVs (colours described in panel b). d, e Number of and ratio between deletions and insertions recovered from the pangenomes and assembly-based mapping (colours described in panel b)

A benchmark dataset to validate SVs is not available for Bovinae, and so we assessed the SV representation accuracy in the three pangenomes comparing against SVs called directly from reference-alignment with the same set of assemblies. We identified between 13.6 k and 68.3 k SVs in the haplotype-resolved assemblies, with substantially more SVs in gaur, yak, and bison than the indicine and taurine haplotypes (Fig. 2b). More than three quarter (78.22%) of the SV alleles were observed in only one haplotype assembly. The eleven haplotype assemblies contained 163 k nonredundant SVs of which 151 k (92.30%), 150 k (91.90%) and 146 k (89.45%) were also recovered with minigraph, pggb and cactus, respectively. This approach validated the vast majority (N = 135 k, 96.26%) of the 141 k SVs that were found in all three pangenomes. However, it validated only 6.12%, 14.06% and 19.24% of the SVs that were respectively private to cactus, pggb and minigraph. Overall, the SV F-score for each minigraph, pggb, and cactus against the assembly truth was 0.908, 0.896, and 0.842 respectively, with normal and repetitive regions outperforming low mappability and satellite regions (Additional file 1: Fig. S6). Similar alleles are collapsed more frequently in minigraph than pggb and cactus, and so the allele frequency spectrum differs between the three tools with minigraph containing substantially less singleton SVs (Fig. 2c). While all tools recovered more insertions than deletions, minigraph contained proportionally less insertions than pggb and cactus (Fig. 2d, e).

We used optical mapping from two Nellore samples [27] unrelated to the NEL haplotype and from the BRA and ANG haplotypes [28], as well as ONT long reads from the OBV, BSW, PIE, NEL, and GAU haplotypes [10] to further validate SV consensus in the pangenomes. Since optical maps call larger SVs (typically minimum SV size of 1 Kb) than sequencing reads, and the ONT datasets had different read lengths, quality, and coverage, the overlap with the pangenome SVs was not expected to be complete. Both sets of orthogonal data substantiate that minigraph and pggb best represent SVs with nearly equal support (1355 and 1283 pangenome SVs overlapping optical map SVs, and ONT SV F-score of 0.799 and 0.802, respectively), while cactus had comparatively lower support (1045 optical map SV overlaps and ONT F-score of 0.743) (Additional file 3: Table S2).

We further assessed the commonality of small variation (SNPs and indels < 50 bp) in pggb and cactus, as minigraph primarily only represents SVs. We recovered 68.27 M nonredundant small variants from the two pangenomes. Pggb and cactus contained 63.02 M (55.37 M SNPs/7.65 M indels) and 64.27 M (56.61 M SNPs/7.66 M indels) small variants of which 59.02 M (53.28 M SNPs/5.74 M indels) were common. The overlap between pggb and cactus was substantially larger for SNPs than indels and multiallelic variants (Fig. 3a–c).

Fig. 3figure 3

Overlap of SNPs and indels between pangenome- and assembly-based discovery. ac Overlap of small variations (SNPs and indels smaller than 50 bp) between pggb, cactus, and assembly-based discovery presented for a biallelic SNPs, b biallelic indels, and c multiallelic SNPs/indels. Variants in different classifications of genomic regions are indicated by colour of the stacked bar. d Number of small variations recovered from each haplotype from the pangenome and the assembly-based mapping. The faded and dark areas of the bars represent indels and SNPs, respectively. e Precision and recall of pggb and cactus for SNPs, indels, and multiallelic variants assuming the assembly-based calls are truth. The black WGS point represents gold-standard accuracy for 30 × sequencing coverage

We assessed small variation representation accuracy in the two pangenomes again using assembly-derived calls as an approximate truth set. The eleven input assemblies contained 62.04 M (54.55 M SNPs / 7.49 M indels) nonredundant small variations. Each assembly had between 5.12 and 27.05 M small variations, again finding substantially fewer in the taurine assemblies than their indicine and non-cattle counterparts (Fig. 3d). Pggb and cactus contained 96.3% and 94.8% of the assembly-based small variants, and 94.8% and 91.5% of the pangenome variation respectively were found in the assembly calls (Fig. 3e). Pggb had a higher overall F-score than cactus, 0.96 and 0.93, respectively, suggesting it encapsulates small variation more accurately. As observed for SVs, accuracies were highest in normal and repetitive regions, although pggb and cactus had F-scores of 0.93 and 0.89, respectively, in tandem repeat regions, indicating a strong ability to resolve repeat motif variation which may only differ by a single or a few bases (Additional file 1: Fig. S6).

We also simulated assemblies for chromosome 29 by introducing known variation per sample to the ARS-UCD1.2 reference genome, removing uncertainty arising from variation that was not called from the real assemblies when aligned to ARS-UCD1.2. Pangenomes constructed from these simulated assemblies had minor improvements to precision and recall for minigraph, pggb, and cactus, as well as for all types of variants (mean F1 improvement for SNP: 0.014, Indel: 0.027, multi-allelic: 0.023, SV: 0.039, Additional file 1: Fig. S7) compared to pangenomes constructed from the real assemblies. The slight improvement suggests that hard-to-call regions or centromeric sequence (not included in simulated assemblies) hurt pangenome accuracy, but that there is still some loss of accuracy in pangenome construction or downstream conversion to VCF even with simulated assemblies.

Alignment of assemblies to pangenomes

We realigned all twelve input assemblies (including ARS-UCD1.2) against the pangenomes to calculate edit distances and quantify sequence content of the pangenomes. Since centromeric, and to a lesser extent telomeric, sequence is challenging to align, we analysed these “centro-/telomeric” regions separately to the rest of the “bulk” sequence (Additional file 4: Table S3). Minigraph had a substantially higher edit rate (0.494%) than pggb (0.010%) and cactus (0.010%) in bulk regions and also in centro-/telomeric regions, 3.43%, 0.396%, and 0.719%, respectively (Fig. 4a). Even though cactus includes more centromeric sequence in the graph (Additional file 1: Fig. S1) and had overall comparable edit rates, its centro-/telomeric edit rates were nearly double that of pggb. Minigraph also had a significantly lower edit rate for ARS-UCD1.2, as all reference sequence is included by definition of minigraph’s algorithm. Aligning to only the linear reference ARS-UCD1.2 resulted in higher edit rates (bulk: 0.527%, centro/telo: 3.56%) and lower query coverage compared to pangenomic alignments.

Fig. 4figure 4

Realignment of included and held-out assemblies to the pangenomes. a Edit rate of the twelve assemblies used in pangenome construction. Each faded dot represents one autosome of each assembly in either the bulk or centro-/telomeric ranges. Mean and 95% confidence intervals are also indicated for each category of sequence per assembly. b Similar plot to a but for eight additional assemblies held-out from pangenome construction. c Query coverage for the bulk sequence for each autosome of all assembly. Values above the dashed line indicate there were more aligned bases than query bases, suggesting multiple ambiguous alignments that were equally scored. d Similar to c, but for the centro-/telomeric sequence. e Bar plots of for the number of alignment errors reported by GraphAligner when using finite “tangle effort”. Bar heights reflect the mean across 348 (included) and 232 (not included) autosomes across the assemblies, while error bars represent 95% confidence intervals

Edit rates were higher for assemblies more diverged from the reference backbone in minigraph pangenomes, indicating they do not incorporate highly diverged segments containing multiple small variations but were not large enough to form bubbles. Contrastingly, pggb and cactus pangenomes did not have any divergence bias for edit rate, suggesting they may be more suited to include more divergent assemblies into super-pangenomes [29].

We also aligned eight additional assemblies held out from the pangenomes (2 × haplotypes from each parent of the OBV and 2 × haplotypes for two unrelated Brown Swiss cattle (Fig. 1)), and found highly similar edit rates for minigraph, but significantly higher edit rates for pggb and cactus (Fig. 4b). This is expected and reflects true genomic variation in the held-out assemblies not present in the pggb/cactus pangenomes, while minigraph never captured small variation anyway and had higher edit rates to begin with. The OD1/2 (dam of OBV sample) assemblies had a lower edit rate in pggb and cactus, as expected since the OBV haplotype is the maternal haplotype.

In addition to substantially lower edit rates, pggb and cactus also covered more of the bulk query sequence (98.4% and 98.6%, respectively) compared to minigraph (93.8%). The differences were more pronounced in the centro-/telomeric query sequence alignments (98.0%, 93.9%, and 65.8%, respectively), as expected given the relative lack of centromeric sequence in minigraph pangenomes (Fig. 4c). Some chromosomes (e.g. 6, 14, and 16) had multiple ambiguous but equally scored alignments in cactus pangenomes, leading to > 100% query coverage. More query sequences can be aligned with more sensitive alignment parameters, but this requires > 300 GB of memory per chromosome in cactus graphs (Additional file 5: Table S4). Even with relaxed alignment parameters, aligning to cactus and pggb pangenomes took 2.2 and 1.4 times more CPU time and 37 and 7 times more memory than to minigraph pangenomes (Additional file 6: Table S5), and was especially prounced for aligning assemblies not included in the pangenomes. Failed alignment of 500 Kb segments was more common in pggb and cactus than minigraph, generally resulting from exceeding GraphAligner’s “tangle effort” in highly complex regions (Fig. 4e). This was especially apparent when aligning not included assemblies to the cactus pangenomes.

Pangenome resolution of VNTR

Variable number tandem repeats (VNTR) account for substantial gene expression and complex trait variation, but due to their repetitive nature and high mutation rates, these loci are difficult to resolve from linear alignments [30, 31]. A catalogue of bovine VNTR had not been established and so we identified 9568 tandem repeats (TRs) in non-masked regions of the ARS-UCD1.2 reference sequence (Additional file 7: Table S6) and investigated their prevalence and variability in the other assemblies through the three pangenomes.

Pggb and cactus contain all assembly sequences, and so can arbitrarily convert between pangenome position and assembly coordinate with tools like odgi. As such, we can liftover the TR positions directly into each assembly through the pangenomes. On the other hand, minigraph contains information on assembly coordinates at graph bubbles, and so we can only estimate TR positions that effectively overlap with SVs. The former approach is substantially more complete and accurate, but incurs a much larger compute cost (Table 2). Approximately 95% of TRs had all twelve assemblies associated with a pangenome path, while the remaining TRs suffered from reduced mappability (Table 2). All pangenomes had instances where coordinates were erroneously translated, but this was most apparent in cactus with several huge outliers (Additional file 8: Table S7). Nearly all TRs examined by minigraph were variable between the assemblies, while approximately 12% of TRs examined by pggb and cactus had zero variation between samples (Table 2).

Table 2 All TRs were examined from pggb and cactus, while TRs overlapping SVs were examined for minigraph. Genotypable TRs are TRs where all 12 assemblies had a known path through the local graph structure. VNTRs are genotypable TRs that had at least one sample with a different number of TR counts. CPU and Memory indicate the compute resources needed to identify the paths of all assemblies through all examined TRs

Tandem repeat counts were similar between the three pangenomes, with 5293 commonly genotyped VNTRs. There were 3332 (63%) VNTRs with identical counts across all twelve assemblies in all three pangenomes and 4084 (77%) identical in at least two pangenomes (Fig. 5a). The remaining VNTR counts still broadly agreed with minor variability (Additional file 1: Fig. S8). The average squared Spearman correlation was 0.92 between the pangenomes, with several outliers in cactus and minigraph skewing the squared Pearson correlation, which was otherwise around an average of 0.8 (Fig. 5b, all significant). TRs present in only pggb and cactus had low variability between cattle and non-cattle (Fig. 5c), suggesting minigraph efficiently captures nearly all VNTRs of interest for further investigation. We also genotyped 465 TRs with adVNTR [31] using HiFi reads from the gaur, Nellore, Piedmontese, Brown Swiss, and Original Braunvieh samples as an approximate truth set, finding good concordance (Fig. 5d, Additional file 7: Table S6), with a median difference of counts of 7, 8, and 11 for pggb, cactus, and minigraph to adVNTR. Minigraph occasionally over- or underestimated TR counts if the overlapping SV was significantly larger or smaller, as that determines the translated coordinates. There were also instances where all three pangenomes agreed on counts different to adVNTR, indicating even HiFi reads may not be powerful enough to genotype all VNTRs (Fig. 5e). Clustering based on TR counts produced trees that broadly grouped the taurine and indicine cattle as well as the non-cattle (Fig. 5f), although the minigraph and cactus trees had slight inconsistencies with the topology of the mash-derived tree (Fig. 1d).

Fig. 5figure 5

VNTR concordance in three pangenomes. a TRs with identical counts in at least two pangenomes, with the median count for the cattle and non-cattle groups. A particular VNTR with substantially more repeats in non-cattle compared to cattle that we investigated further is circled in red. b Pearson and Spearman squared correlation coefficients across the TR counts for pggb-cactus (p–c), pggb-minigraph (p-m), and cactus-minigraph (c-m). Each point is one assembly, with box plots over the 12 assemblies. c Similar to a, except TRs with identical counts in pggb and cactus that were not present in minigraph. d adVNTR-derived genotypes for five HiFi samples in the three pangenomes. The black dashed line indicates the expected count using adVNTR as a ground truth. e VNTRs where all three pangenomes agreed to a different count than adVNTR, suggesting adVNTR may sometimes over-/underestimate assembly-based counts. f Trees derived from the TR counts across different input assemblies, with colours representing clusters of taurine cattle, indicine cattle, and non-cattle

An eVNTR mediates expression of neighbouring genes and non-coding RNA

The pangenomes identified a highly variable VNTR locus on chromosome 12 (86,161,072–86,162,351 bp), containing substantially more copies of a degenerate 12 bp motif in non-cattle and fewer in Nellore than in the taurine assemblies (Fig. 6a, b), prompting a detailed investigation. Although the genotyping was similar, bandage plots revealed significantly different graph structures of the VNTR across the three pangenomes (Fig. 6c). Minigraph added additional nodes to incorporate the (primarily non-cattle) additional tandem repeats, while pggb, and cactus to a lesser extent, incorporated the additional repeats by looping through the nodes containing the tandem repeat sequence multiple times (Additional file 1: Fig. S9). This different representation was observed generally in VNTRs with high repeat counts (Additional file 1: Fig. S10).

Fig. 6figure 6

A polymorphic eVNTR discovered from the pangenomes. a VNTR repeat motif and b number of repeat units in the twelve input assemblies. c Bandage plots of a VNTR upstream SPACA7. BLAST hits for the top 5 most common VNTR motifs are coloured per motif. d Genes (black) and (long) non-coding RNAs (dark grey) nearby the VNTR. Blue colour indicates the position of the VNTR. Arrows indicate the orientation of genes and lncRNA. Red colour indicates two genes and five lncRNAs whose expression is associated with the VNTR. ek Expression (quantified in transcripts per million (TPM)) of the associated genes and ncRNAs in testis tissues of 83 Brown Swiss bulls that are either homozygous for hap1 (H1/H1), homozygous for hap2 (H2/H2) or heterozygous (H1/H2). The number of animals per diplotype is below the boxplots. l Cis-expression QTL mapping for LOC112449094. Different colours indicate the pairwise linkage disequilibrium (r2) between the VNTR haplotype (violet) and all other variants. m Correlation between the abundance of mRNA and lncRNA for the associated genes and lncRNA. n Allelic imbalance of nine heterozygous exonic SNPs in LOC112449094 in testis tissue from a Nellore x Brown Swiss crossbred bull. Orange and blue colours represent paternal (Nellore) and maternal alleles (Brown Swiss)

We observed two VNTR “haplotypes” (referred to as hap1 and hap2) in long-read alignments of additional Original Braunvieh cattle demonstrating within-breed variability. While the alignments indicated several insertions and deletions in both haplotypes with respect to ARS-UCD1.2, hap2 is 24 bp longer than hap1 due to two additional copies of the repeat motif (Additional file 1: Fig. S11). There are 63 genes annotated within the cis-regulatory range (± 1 Mb) of the VNTR, of which 44 (Fig. 6d) were expressed at > 0.2 TPM in testis tissues of 83 mature Brown Swiss and Original Braunvieh bulls.

A degenerate repeat motif and an overall length > 1300 bp precluded the short sequencing read-based genotyping of the VNTR in the eQTL cohort with adVNTR. Instead, we genotyped the VNTR through two SNPs (Chr12:86,160,984 and Chr12:86,160,971) and one indel (Chr12:86,161,000) that tagged the two VNTR haplotypes. This enabled us to investigate putative cis- and trans-regulatory impacts of the VNTR on the expression of genes and long non-coding RNAs (lncRNA). Two genes (TUBGCP3, SPACA7) and five non-coding RNAs (LOC112449093, LOC112449094, LOC112449095, LOC112449100, LOC112449104) within ± 1 Mb of the VNTR were differentially expressed (P < 1.13 × 10–3, Bonferroni-corrected significance threshold) between the diplotypes indicating a putative cis-regulatory role (Fig. 6e–k). We did not detect trans-regulatory effects for the haplotype (Additional file 1: Fig. S12).

We mapped eQTL within the cis-regulatory range (± 1 Mb of the transcription start site) of the two significant genes and five significant non-coding RNAs to investigate if variants in cis other than the VNTR haplotype are associated with transcript abundance. The VNTR haplotype was amongst the top variants at the LOC112449094 eQTL, but eleven variants in strong linkage disequilibrium (r2 = 0.93) were more significant (P = 2.00 × 10–22vs. P = 5.18 × 10–20) (Fig. 6l). The eQTL peak was absent when the VNTR haplotype was fitted as a covariate. We made similar observations for the TUBGCP3 eQTL (Additional file 1: Fig. S13). Variants other than the VNTR haplotype were more strongly associated with the expression of SPACA7, LOC112449093, LOC112449095, LOC112449100, and LOC112449104 (Additional file 1: Fig. S13). These findings suggest that the VNTR haplotype primarily mediates the expression of LOC112449094 and TUBGCP3.

Hap2 containing more copies of the 12 bp motif increases abundance (βnormalized = 1.26 ± 0.10, P = 5.18 × 10–20) of LOC112449094 which is lowly-expressed (0.86 ± 0.31 TPM) in 83 bull testis transcriptomes. Conversely, hap2 is associated (βnormalized = –1.08 ± 0.14, P = 3.81 × 10–11) with lower TUBGCP3 mRNA. Negative correlation (–0.37, Fig. 6m) between LOC112449094 and TUBGCP3 abundance suggests that LOC112449094 could be a cis-acting lncRNA that represses expression of TUBGCP3. The VNTR haplotype is in strong LD with two SNPs (Chr12:86,173,431 (r2 = 0.87) and Chr12:86,173,509 (r2 = 0.93)) in LOC112449094 exons. While hap2 primarily segregates with the alternate alleles of these two SNPs, hap1 segregates with the respective reference alleles. Alternate allele support in RNA sequencing reads overlapping Chr12:86,173,431 and Chr12:86,173,509, respectively, was 71% (96 out of 135 alleles, Pbinom = 1.01 × 10–6) and 69% (94 out of 137 alleles, Pbinom = 1.57 × 10–5) in 22 animals that are heterozygous both at the VNTR haplotype and the SNPs confirming allelic imbalance due to lower LOC112449094 expression in hap1.

We confirmed a putative cis-regulatory role of the VNTR in an individual with indicine ancestry. Expression of LOC112449094 was relatively low (TPM = 0.66) in the testis tissue sampled from a Nellore (Bos taurus indicus) x Brown Swiss (Bos taurus taurus) crossbred bull (NxB), although it inherited hap2 from its Brown Swiss dam (Additional file 1: Fig. S11). The paternal (Nellore) haplotype is diverged from hap2 and hap1 and contains less VNTR repeat units (Fig. 6b). Alignment of parental-binned HiFi reads against ARS-UCD1.2 readily distinguished paternal from maternal alleles at nine heterozygous SNPs overlapping LOC112449094 exons. Inspection of these variants in the RNA sequencing alignments of NxB F1 showed allelic imbalance due to a disproportionally low number or even complete absence of paternal (= Nellore) alleles (Fig. 6n) suggesting that the low VNTR repeat count of the Nellore haplotype represses expression of LOC112449094.

留言 (0)

沒有登入
gif