Long read sequencing characterises a novel structural variant, revealing underactive AKR1C1 with overactive AKR1C2 as a possible cause of unexplained severe fatigue

doi:10.21203/rs.3.rs-3218228/v2

Download PDF

Research Article

Long read sequencing characterises a novel structural variant, revealing underactive AKR1C1 with overactive AKR1C2 as a possible cause of unexplained severe fatigue

https://doi.org/10.21203/rs.3.rs-3218228/v2

This work is licensed under a CC BY 4.0 License

Journal Publication

published 16 Nov, 2023

Read the published version in Journal of Translational Medicine →

Version 2

posted

You are reading this older preprint version

Read the latest preprint version →

Background

Causative genetic variants cannot yet be found for many disorders with a clear heritable component, including chronic fatigue disorders like myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS). These conditions may involve genes in difficult-to-align genomic regions that are refractory to short read approaches. Structural variants in these regions can be particularly hard to detect or define with short reads, yet may account for a significant number of cases. Long read sequencing can overcome these difficulties but so far little data is available regarding the specific analytical challenges inherent in such regions, which need to be taken into account to ensure that variants are correctly identified. Research into chronic fatigue disorders faces the additional challenge that the heterogeneous patient population likely encompasses multiple aetiologies with overlapping symptoms, rather than a single disease entity, such that each individual abnormality may lack statistical significance within a larger sample. Better delineation of patient subgroups is needed to target research and treatment.

Methods

We use nanopore sequencing in a case of unexplained severe fatigue to identify and fully characterise a large inversion in a highly homologous region spanning the AKR1C gene locus, which was indicated but could not be resolved by short-read sequencing. We then use GC-MS/MS serum steroid analysis to investigate the functional consequences.

Results

Several commonly used bioinformatics tools are confounded by the homology but a combined approach including visual inspection allows the variant to be accurately resolved. The DNA inversion appears to increase the expression of AKR1C2 while limiting AKR1C1 activity, resulting in a relative increase of inhibitory neurosteroids and impaired progesterone metabolism.

Conclusions

This study provides an example of how long read sequencing can improve diagnostic yield in research and clinical care, and highlights some of the analytical challenges presented by regions containing tandem arrays of genes. It also proposes a novel gene associated with a specific disease aetiology that may be an underlying cause of complex chronic fatigue and possibly other conditions too. It reveals biomarkers that could be assessed in a larger cohort, potentially identifying a subset of patients who might respond to treatments suggested by the aetiology.

AKR1C1

AKR1C2

Fatigue

Neurosteroids

Allopregnanolone

Nanopore Sequencing

Long Read Sequencing

Structural Variants

The dramatic reduction in the cost of high-throughput sequencing over recent years has made short-read whole genome sequencing (WGS) viable as a go-to option for clinical genomic analysis. The high accuracy at base level combined with increasingly sophisticated bioinformatic analytical tools has enabled diagnosis in many cases that were previously unexplained. Many other patients, however, are left still searching for answers.

Disorders of unexplained chronic fatigue (CF) such as myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) show strong evidence of heritability but research is hampered by the lack of clear diagnostic biomarkers [1, 2]. The heterogeneous nature of symptoms, triggers and disease progression suggests that the broad diagnostic criteria likely encompass multiple different aetiologies with overlapping symptoms which may confound any investigation of the population as a whole [1, 2]. These debilitating, multi-system conditions can leave people bedbound for years, unable to speak or feed themselves, often in constant pain and unable to tolerate any light, sound or touch [3]. There is an urgent need to distinguish between subsets of patients in order to better target research and treatment.

There is hope that the rise of next-generation WGS and large-scale genome-wide association studies (GWAS) will bring about a step change in understanding. However, financial constraints mean that GWAS must generally rely on a single nucleotide variant (SNV) array based approach, which by definition cannot identify novel SNVs or structural variants. Next-generation WGS may also fall short because some parts of the genome are inherently “dark” to short-read alignment due to high homology or low complexity, where entire reads can map equally well to multiple loci. These regions may be more prone to structural variation and are thought to account for a significant proportion of pathogenic variants [4].

One such region is a stretch of segmental duplication on human chromosome 10p15 to 10p14 that contains the four Aldo-Keto Reductase Family 1 Member C genes. AKR1C1-4 are arranged in a head-to-head, tail-to-tail manner with >86% amino acid sequence identity overall. AKR1C1 and AKR1C2 in particular share nearly 98% amino acid sequence identity and differ by only 7 amino acids [5]. The AKR1C enzymes are well known for their wide-ranging roles in steroid metabolism, including deactivation of progesterone (P) and synthesis of the neuroactive steroids allopregnanolone (allo) and pregnanolone (preg), which positively modulate the gamma-aminobutyric acid type A receptor (GABA_AR) in a similar way to benzodiazepine sedatives [6]. To date, however, germline variants in AKR1C1 and 1C2 have been surprisingly little studied in the context of disease [7].

Long read sequencing technologies such as those from Oxford Nanopore (ONT) and Pacific Biosciences (PacBio) have the potential to shed light on these dark regions. Smaller structural variants (SVs) of up to around 50 kb may be covered in their entirety by single reads, while with larger variants the reads spanning breakpoints are likely to cover enough unique bases to enable accurate alignment, eliminating the multiple-mapping problems of short reads. Early use of clinical long read sequencing has successfully identified a number of pathogenic variants missed by standard testing, but as yet there are few reports detailing the specific challenges that repetitive regions present for algorithmic analysis or how well the most-used bioinformatic tools perform in practice [8]. A better understanding of where automated pipelines may struggle could increase diagnostic yield by ensuring that fewer variants are missed.

Here, we demonstrate the power of long read sequencing to fully characterise a large inversion in a difficult-to-align region that is all but invisible to short reads. We highlight some of the challenges that can arise when analysing both long read and short read data from such regions, and show that even when SV callers do indicate structural variation they may not correctly resolve the SV, potentially suggesting a functional impact that can differ significantly from the true effect.

We also present biochemical evidence that the variant described here is functional in our patient. We propose that it represents a novel aetiology combining underactive AKR1C1 with overactive AKR1C2, which limits the clearance of intracellular progesterone while suppressing neuronal activity via the GABA_AR, resulting in a complex, multi-system syndrome of severe fatigue. We further propose that the genetic and biochemical biomarkers discovered in our patient could be assessed in a larger cohort, to see if they constitute a subset of CF patients who might respond to treatments suggested by the aetiology.

Clinical Details

The 42 year old patient has a long and unusually complex history, of which the most debilitating aspect is increasingly severe mental and physical fatigue. Onset occurred after a viral infection age 16, initially mild but gradually progressing. At present, the patient has been completely housebound for many years and is unable to sit up or speak for much of the day.

Clinical diagnoses include central adrenal insufficiency and partial central diabetes insipidus, both thought to be hypothalamic rather than pituitary, along with renal salt wasting and functional pan-enteric dysmotility. There was previously a diagnosis of ME/CFS for some years but that label can be a barrier to getting appropriate treatment so was removed after the later diagnosis of other conditions known to cause persistent fatigue. Mental health assessment concluded that there was no underlying psychiatric disease. None of the current diagnoses, however, fully explains the severity of illness or the full range of symptoms. It is also highly unusual to see this constellation of conditions, especially with young-adult onset. The patient was referred to the UK’s 100,000 Genomes Project (100kGP) on suspicion of Gitelman syndrome, but no likely tubulopathy-associated variants have been found and other aspects of her condition are outside the scope of the 100kGP.

Of particular note was that certain steroid medications provoked an atypical reaction resembling the effects of benzodiazepines, leading to a hypothesis of excessive neurosteroid-mediated activation of the GABA_AR. The patient reported that a similar but less acute sense of sedation, distinct from the fatigue, had started some years previously. It had originally occurred only in the latter stages of the menstrual cycle but had progressively spread throughout the cycle until it is now constantly present, although still with some cyclic variation.

The family history included unexplained fatigue in several family members, albeit to a much lesser degree, including a period of quite debilitating fatigue in the mother’s early 20s and a maternal second cousin diagnosed with ME/CFS.

RQ-PCR was conducted for a separate study investigating the AKR1Cs in relation to another condition that the patient also has in a mild form, which she joined after the discovery of the inversion in the 100kGP data. The results from a luteal phase sample showed substantial overexpression of AKR1C2 mRNA, with relative expression approximately 100 times higher than controls, while AKR1C1 and AKR1C3 were comparable with controls (full details seen but not reproduced here as the other study has not yet reported).

Short-read WGS

Whole genome sequence data for the patient and the patient’s mother were provided by the 100kGP. The data included Illumina short reads aligned to GRCh38 with small variants called by Starling (Illumina proprietary workflow), plus structural variants called by Manta [9] and Canvas [10].

Long Read WGS

DNA was extracted from 1 mL of whole blood with the Nanobind CBB Big DNA kit (Circulomics NB-900-001-01). After extraction, the genomic DNA (gDNA) was homogenised by 3-10 passes of needle shearing (26G) and 1h of incubation at 50°C. The gDNA was then split into two aliquots. 4 μg of gDNA were fragmented to ~50kb with the Megaruptor 3 (Diagenode, E07010001 and E07010003) and 10 μg were fragmented to 20-25kb with a gTUBE (Covaris, 520079). All samples were depleted of short DNA fragments (less than 10kb long) with the SRE kit (Circulomics, SS-100-101-01). DNA size distributions were assessed at each step on a FemtoPulse system (Agilent, M5330AA and FP-1002-0275) and all samples were quantified on a Qubit fluorometer (Invitrogen, Q33226).

Two sequencing libraries were generated from 1.5 μg of each gDNA samples with the SQK-LSK109 kit (Oxford Nanopore Technology). 25 fmol of the obtained sequencing libraries were loaded onto a R9.4 flow cell (ONT) and sequenced on a PromethION48 (ONT).

Base-calling was performed with Guppy-4.0.11 using the high accuracy model and a qscore filter of 7. Full details of the protocol can be found at https://research-help.genomicsengland.co.uk/display/GERE/Genomic+Data+from+ONT?preview=/38046759/38047942/v1_protocol_ONT_LSK109.pdf. Sequencing reads were aligned to GRCh38 using minimap2 version 2.17 [11]. QC statistics and plots were generated using Nanoplot version 1.26.0 [12]. Structural variants were called with sniffles version 1.0.11 [13]. The full details of bioinformatics pipeline can be found at https://research-help.genomicsengland.co.uk/display/GERE/Genomic+Data+from+ONT?preview=/38046759/38047944/PromethION%20SV%20calling%20pipeline%20GRCh38.docx.

Bioinformatic analysis

The variant call files provided by the 100kGP were annotated and filtered using seqr [14].

For the ONT data analysis we combined visual inspection in Integrative Genomics Viewer (IGV) [15] with tools including Samtools [16], Minimap2, Flye [17], Racon [18] and Liftoff [19] to define the inversion, phase and assemble the long reads, and to compare the long read data to the short reads.

Serum steroid analysis

The patient provided five samples in total: two from the follicular phase of the menstrual cycle and three from the luteal phase, covering different cycles with symptomatic variation. All samples were taken in the morning. Blood was drawn from the cubital vein, centrifuged at 1790g for 10 minutes, and the serum stored at -20 °C until analysis. The samples were prepared and analysed using our previously described GC-MS/MS method [20] and 82 analytes reported. Concentrations were then compared against a pool of corresponding controls from our records (11 follicular, 19 luteal).

Short-read WGS indicates a large inversion around AKR1C1 and AKR1C2, with a splice-site variant in the 5’UTR of AKR1C2

The Manta SV calls from the 100kGP indicated a heterozygous large inversion of a region from approximately 5kb upstream of AKR1C1 to part way through the extended 5’UTR of the first two transcripts of AKR1C2 (Fig. 1). This variant seemed potentially relevant since AKR1C1 and AKR1C2 are both instrumental in the synthesis of GABAergic neurosteroids. If the inversion alters protein expression in such a way that AKR1C2 is overexpressed and/or AKR1C1 is partially lost then it would seem to fit precisely with the clinical presentation.

In addition to the inversion call we identified a heterozygous splice-donor loss single nucleotide SNV, NG_031852.1:g.15610G>C (rs67245807), in the second intron of the first transcript of AKR1C2, marked with an arrow in Fig. 1. This intron is part of the 5’UTR of AKR1C2 just inside the inverted region.

Both variants were also present in the patient’s mother, suggesting that they were on the same allele. The splice-site SNV was therefore deemed unlikely to be relevant, since the inversion effectively removed the first exon of the extended 5’UTR meaning that transcript 1 would not be produced. This SNV is also relatively common, with an ALFA Allele Frequency of 0.12712 in Europeans (lower in other populations). It is not reported in ClinVar.

Detailed analysis of the inversion was not possible with the short-read data due to the high homology. In addition, similar large inversion calls were present in a number of unrelated individuals in the 100kGP dataset, so that it was unclear whether the inversion was actually real or merely an alignment artefact caused by the homology. Long read sequencing was therefore sought instead.

Characterising the variant from the long read data

Initial overview

We generated two different libraries using DNA extracted from blood, with N50s of 23kb and 37kb, and ran one flowcell for each library. Both flowcells generated approximately 90Gbases of mappable data, providing a total coverage in the region of 60x.

SV calling with Sniffles indicates several insertions and deletions of between 50 and 420 kb in the region of AKR1C1 and AKR1C2 but does not detect an inversion. On visual inspection in IGV, however, signs of a heterozygous large SV spanning the region can clearly be seen, although the exact nature of the SV is not immediately obvious (Fig. 2).

After initial phasing with Samtools Phase it can be seen that the REF allele looks normal while the ALT allele has several features suggesting an inversion, including mirrored insertions and deletions where a segment is deleted from one side of the region and the reverse complement of that segment inserted at the equivalent position on the other side. These are the variants that Sniffles had detected. Due to the high homology, however, there is not the classic inversion pattern of reads split into supplementary alignments on opposite strands. Instead, the reads remain intact but with clusters of high-density SNVs interspersed with sections that appear normal, as well as the mirrored insertions and deletions.

Confirming the inversion

Breakpoints were estimated by the loci at which variants in the ONT sequence started to differ significantly from those in the 100kGP sequence. To confirm the inversion we created a new reference sequence consisting of the GRCh38 consensus sequence with the region between the proposed breakpoints inverted. When mapped to this inverted reference the ALT allele fits cleanly while the REF allele now shows the same pattern of variants as the ALT allele had when aligned to the standard reference, but with the pattern reversed (Additional file 1 - Fig. S1). The span of the inversion was therefore established as chr10:4962657-5004621 - i.e. 758bp before exon 1 of AKR1C1 and 602bp upstream of the first coding exon of AKR1C2.

Homology hinders phasing

The phasing itself presents some challenges. The extreme homology means that several reads contained entirely within the inversion cannot be confidently assigned to either allele at the first attempt, then once aligned to the inverted reference it is clear that some of the reads that had been confidently assigned are actually on the wrong allele. The phasing was therefore refined with two further rounds of Samtools Phase on each allele and the final few aberrant reads were moved manually using basic bash commands. This produces two cleanly phased BAM files, one of which closely matches the standard reference while the other closely matches the inverted reference sequence (Additional file 1- Fig. S1).

Additional insertion near AKR1C2

In addition to the large inversion, the long reads reveal a novel sequence insertion of 310-330bp on the inverted allele at chr10:4961564-4961565, 1092bp upstream of the first breakpoint. It is not possible to determine the precise length due to a homopolymer run at one end, the length of which varies significantly between the reads. This insertion appears to be nsv4544854, which can also be detected visually in the 100kGP data (though it hadn’t been called), but in a slightly different position because repetitive sequence at the end affects short-read alignment.

Assembling the alignments

Attempts to assemble the ALT allele using Canu and Minimap2-Miniasm appeared to be confounded by the homology in the region but good results were obtained using Flye, either with 5 rounds of polishing as part of the Flye pipeline or with one round of polishing by Flye followed by one round of polishing with Racon. Accuracy was assessed by mapping the raw reads back to the assembly for visual inspection in IGV (Additional file 1- Fig. S2). We then followed the same process to assemble the REF allele for comparison.

Liftoff was used to map the RefSeq curated subset of genes and functional element annotations (May 2021) to each of the assembled alleles. As can be seen in Fig. 3, the REF allele shows the expected position and orientation of AKR1C6P, AKR1C1, AKR1C2, AKR1C3 and AKR1C8P, while in the ALT allele AKR1C1 has been swapped with the coding region of AKR1C2.

Comparison with short-read sequence

The inversion revealed by the ONT sequence is nearly 10kb shorter than suggested by the short read data. To investigate this discrepancy we created two simulated reads from the 10kb sections immediately outside of the breakpoints at each end of the ALT assembly. These simulated reads were then aligned back to the GRCh38 reference both in their correct positions and at the opposite sides (Fig. 4a). This enabled a direct comparison of the two sides of the ALT allele, both with each other and with the short reads.

Approximately 5kb out from the real breakpoints, where Manta had called an inversion due to split and RR/LL oriented reads, there is a short section of 31bp where the sequence from the AKR1C1 side appears to have been copied and pasted over the equivalent position on the AKR1C2 side (Fig. 4b). Both of those loci therefore match the reference at the AKR1C1 side, which is what had caused the short reads to misalign. All of the apparently split or wrongly oriented read pairs can actually be aligned to GRCh38 in a normal configuration that matches the ALT assembly.

Comparing the simulated reads immediately outside of the breakpoints also revealed why the short-read sequence didn’t detect the real inversion: in the reference sequence the first approx. 570 bases just outside each end of the inversion are almost identical, with only 3 bases different; in the ALT allele those regions are completely identical, with both sides the same as the reference on the AKR1C1 side. In other words, it appears to be another copy-and-paste. There is then only minimal variation for nearly 300bp after that.

Serum steroid analysis suggests a functional effect with increased AKR1C2 and reduced AKR1C1 activity

One of the patient’s two follicular phase samples was from a normal cycle while the other was from what she described as a “hangover cycle”, where luteal symptoms persist for several days into the follicular phase. This was reflected in the results, with the first having a P concentration comparable with controls at 0.156 nmol/L, while in the “hangover” sample it was much higher at 4.06 nmol/L. Luteal samples also showed significant variation in absolute steroid concentrations but the relative ratios between steroids were largely similar across the luteal samples.

Relative concentrations of selected serum steroids (Table 1) were calculated to overcome concentration differences and to reflect the relative activities of AKR1C1 and AKR1C2. Conjugated steroids were high relative to the corresponding unconjugated steroids across the board in the patient compared to controls. We therefore used the sum of free + conjugated steroid, where conjugates were present, as better representing total production. The progesterone-derived AKR1C2 substrates 5α-dihydroprogesterone (5α-DHP) and 5β-dihydroprogesterone (5β-DHP) are not included because their concentrations were below the limit of quantification.

Abbreviation	Steroid	Significance
Allo	Allopregnanolone	Inhibitory neurosteroid
Preg	Pregnanolone	Inhibitory neurosteroid
Iso	Isopregnanolone	Counter-inhibitory neurosteroid
P	Progesterone	Precursor of 20α-reductase (AKR1C1)
17-OHP	17α-Hydroxyprogesterone	Precursor of 20α-reductase (AKR1C1)
20α-DHP	20α-Dihydroprogesterone	20α-reduced metabolite of P
17,20α-diOHP	17α,20α-Dihydroxy-4-pregnen-3-one	20α-reduced metabolite of 17-OHP
DHT	5α-Dihydrotestosterone	Precursor of 3α-reductase (AKR1C2)
Adiol	3α-Androstanediol	3α-reduced metabolite of DHT

Table 1 Steroids for ratios to assess the functional impact of the inversion. Conjugates were included for all steroids except progesterone and 17α-hydroxyprogesterone, which do not undergo significant conjugation.

Results are shown in Fig. 5. The bars represent medians, with results for the patient individually marked and error bars for the controls showing the interquartile range, to allow for skewed data distribution and non-constant variance.

Neurosteroid homeostasis

The concentration of the main inhibitory neurosteroids, allo and preg, is raised relative to the main counter-inhibitory neurosteroid, iso (Fig.5a).

Precursor-product ratios of AKR1C1 (20α-reduction)

The ratios of 17-OHP/17,20α-diOHP (Fig. 5c) and P/20α-DHP (Fig. 5d) were used to assess the 20α-reduction activity of AKR1C1. 17-OHP/17,20α-diOHP is raised, except in one luteal sample. This sample contained a relatively high cortisol concentration due to higher adrenocortical stimulation, which would have resulted in increased utilisation of its precursor, 17-OHP, leaving less available for 20α-reduction (explained further in Additional file 2). P/20α-DHP is raised in the follicular phase but normal or slightly low in the luteal phase.

To see if the normal-low luteal P/20α-DHP ratio could be due to excess P being diverted down the 5α pathway instead of the 20α pathway we examined the sum of all 5α progesterone metabolites versus the sum of all 5β progesterone metabolites, as an indication of relative 5α-reductase (5α-R) activity (Fig. 5e). The 5α/5β ratio for progesterone metabolites (5α/5β P) is high in the luteal phase.

Precursor-product ratio of AKR1C2 (3α-reduction)

The ratio of DHT/Adiol, an indicator of 3α-reduction by AKR1C2, is low (Fig. 5b).

Long-read sequencing has not yet filtered through to routine clinical practice, partly due to the high cost and analytical challenges. Variants in difficult-to-align regions may not show the same patterns, such as split or supplemental reads, that would be seen with the same variant in a more straightforward part of the genome. The different patterns may cause SVs to be missed by some variant callers, and while algorithmic tools are rapidly improving they can still struggle, potentially introducing errors that would hinder the downstream analysis or affect interpretation of the variant. In this case we found that despite high coverage, phasing was only partly successful due to the homology; several rounds of Samtools Phase and additional manual manipulation were required to define fully phased alleles that could then be assembled. Even with clean phasing, the homology still confounded some commonly used assemblers. Despite the challenges, however, long-read sequencing can already offer real hope for desperate patients by illuminating regions where short reads fail.

In our case the original short-read sequence did indicate an inversion but what it had picked up was, in fact, an alignment artefact due to a small copy-and-paste error some 5kb further out. The true inversion is effectively invisible to short reads since any reads spanning the breakpoints would align directly to the other end, while only those pairs with the longest inserts would reach from inside the inversion to the more divergent region further out and so be correctly aligned to opposite ends. In the patient’s short-read sequence there are two such correctly aligned pairs, in the mother's there are none at all. Of course, that also means there's no way of knowing how common this inversion is, since it would not be detected in any of the short-read-based datasets currently available.

The narrower span of the true inversion alters interpretation of the splice-donor loss SNV, which is now seen to be outside of the inversion and therefore part of the extra 5’UTR attached to AKR1C1. The narrower span also makes it more likely that the inversion is functional, since the breakpoints are much closer to the coding regions of both genes leaving more of the regulatory regions outside.

AKR1C2 has 4 coding transcript variants in the Consensus Coding Sequence (CCDS): transcripts 1 & 2 (CCDS7062) differ only in their 5’UTR and both encode the full AKR1C2 protein of 323 amino acids (isoform 1). Transcript 3 (CCDS44350) encodes isoform 2, a truncated version of only 139 amino acids that is predicted to be inactive. Transcript 4 (CCDS81438) lacks the extended 5’UTR of transcripts 1 & 2 and is also missing one exon from the coding region; it encodes the 297 amino acid isoform 3. The RefSeq curated subset also includes a fifth transcript variant, which encodes the full isoform 1 of transcripts 1 & 2 but without the extended 5’UTR. AKR1C1 has only one transcript in the CCDS.

The high homology on opposite strands gives this inversion rotational near-symmetry and effectively swaps AKR1C1 with the coding region of AKR1C2. This has a number of consequences:

1. The initial non-coding exons of AKR1C2 transcripts 1 & 2 are still in their correct place, meaning that transcription could presumably be initiated in the normal way and so potentially result in a rearranged AKR1C1 mRNA with those additional exons added at the 5’ end, and possibly a retained intron immediately upstream of the first coding exon. It seems likely that this would affect expression and/or translation of AKR1C1.

2. Without those first non-coding exons, AKR1C2 would only have transcript variants 3, 4 and 5. Transcripts 4 and 5 would presumably still produce a functional protein but it seems likely that expression would be altered.

3. The orientation and positions of AKR1C1 and 1C2 are reversed relative to AKR1C6P, AKR1C3, AKR1C8P, and to all of the regulatory regions outside of the breakpoints. These include three transcriptional cis-regulatory regions on the forward strand in the RefSeq curated subset, marked A, B and C in Fig. 3, that have been shown to positively regulate AKR1C1 with high confidence in K562 cells. The middle of these (B) also positively regulates AKR1C2. GeneHancer indicates other regulatory elements with interactions potentially disrupted by the inversion, as does the ORegAnno database (both shown in Fig. 6), and the MoLoTool sequence motif location tool [21] used with the assembled alleles suggests several candidate steroid hormone receptor binding sites close to the breakpoints.

Expression of the different AKR1C proteins and isoforms is tissue-specific, with AKR1C4 largely confined to the liver while 1C1-3 are expressed in multiple tissues with distinct distribution patterns [5, 23]. The main roles of the AKR1Cs include regulation of active and inactive intracellular steroid hormone concentrations, prostaglandin metabolism, synthesis and binding of bile acids, and as Phase I metabolising enzymes including as part of the antioxidant defence system [5]. All of the AKR1C enzymes can accept a wide range of substrates and exhibit 3-, 17- and 20-ketosteroid reductase activity, but with different preferences for substrate, stereochemistry and reduction position [5].

AKR1C2 preferentially acts as a 3α-reductase, and is the primary enzyme for synthesis of the inhibitory GABAergic neurosteroids allo, preg and allotetrahydrodeoxycorticosterone (THDOC) from 5α-DHP, 5β-DHP and 5α-dihydrodeoxycorticosterone respectively [5, 24]. Allo and preg in particular are potent positive modulators of the GABA_AR, making the receptor more responsive to circulating GABA, which in turn suppresses signalling in both the central nervous system and peripheral tissues and so has wide-ranging effects, particularly in the brain [25–27].

AKR1C1 acts predominantly as a 20α-reductase (20α-R) that reduces active steroids to their inactive metabolites, particularly P to 20α-DHP. AKR1C1 is also thought to be the main enzyme producing the 3β metabolites of 5α-reduced steroids [28]; these include iso, which counters the effects of allo and preg at the GABA_AR [26, 27].

As shown in Fig. 7, if AKR1C1 is underactive and AKR1C2 overactive then more progesterone will likely be available for the 5α- and 5β-reduction pathways instead of being deactivated to 20α-DHP, and will then be 3α-reduced to allo and preg rather than 3β-reduced to iso, disrupting neurosteroid homeostasis.

The overlapping activities of the AKR1Cs and the complex pathways of steroid metabolism more widely can make it challenging to assess individual enzymes from peripheral blood. Hepatic AKR1C4 in particular has high catalytic efficiency for both 20α- and 3α- reduction [5] and would likely compensate for either excess or deficiency in the other enzymes. Furthermore, it is known that locally synthesised steroids, e.g. in the brain, may reach functionally significant levels without any corresponding variation in peripheral blood [27]. However, neurosteroids can cross the blood-brain barrier so that fluctuations in circulating concentrations, e.g. due to the menstrual cycle, can be reflected in the brain [29], and serum results may still indicate a general tendency that could potentially be more pronounced in other tissues.

The high ratio of (allo+preg)/iso directly supports the hypothesis of a relative increase in inhibitory GABAergic neurosteroids in our patient, with the higher value in the luteal phase corresponding to increased symptoms in the latter part of the menstrual cycle.

The high ratios of 17-OHP/17,20α-diOHP and follicular P/20α-DHP are consistent with partial loss of AKR1C1. The normal-low luteal P/20α-DHP ratio is surprising but could be explained by compensatory diversion of P down the 5α pathway.

5α-R is ubiquitous while the majority of 5β-reduction occurs in the liver [30], meaning that the 5α/5β ratio provides a measure of where 5-reduction is occurring as well as indicating possible alterations in enzyme capacity. The high 5α/5β ratio for P metabolites in the luteal phase suggests a shift to extrahepatic 5-reduction of P. Taken together, the P/20α-DHP and 5α/5β ratios suggest that intracellular P may be accumulating in extrahepatic tissues due to limited 20α-reduction by AKR1C1 and is then being metabolised by 5α-R instead.

The low ratio of DHT/Adiol is consistent with increased 3α-reduction by AKR1C2, with the caveat that adiol can also be synthesised from allo, via androsterone, bypassing DHT and potentially altering the DHT/Adiol ratio. The overexpression of AKR1C2 mRNA in luteal phase also supports the hypothesis of increased AKR1C2.

If this variant is indeed causing overactivity of AKR1C2 while limiting AKR1C1 then it could be part of the underlying cause not only of the fatigue and sedation, due to the sedative effects of allo and preg, but also of the endocrine, renal and gastrointestinal dysfunction in this patient. Excessive GABAergic inhibition in the hypothalamus could suppress the release of corticotropin releasing hormone (CRH) and antidiuretic hormone (ADH), resulting in cortisol deficiency with diabetes insipidus [31, 32]. AKR1C1 (but not 1C2) is highly expressed in the proximal tubule of the kidney [33, 34], where impaired progesterone clearance could interfere with sodium reabsorption through inhibition of the mineralocorticoid receptor and other electrolyte-regulating proteins, resulting in renal salt wasting [35, 36]. AKR1C1, 1C2 and the GABA_AR are all active in the gastrointestinal tract, where GABA is involved in regulating secretion and motility [33, 34, 37].

More broadly, variants in AKR1C1 and 1C2 could potentially affect any tissue where steroid hormones undergo 20α- or 3,5-reduction, likely in a tissue-specific manner, with implications for a variety of conditions where the aetiology is unknown or a subset of cases remains unexplained. This would be particularly the case for multi-system chronic fatigue disorders like ME/CFS, fibromyalgia and long covid. Gastrointestinal dysfunction and hypocortisolism are common in these conditions, as is dysregulation of the autonomic nervous system [38–40] which may similarly be modulated by neurosteroids [41, 42]. These conditions are also more common in women [43, 44], as might be expected if higher P is a factor. We therefore plan further investigation including RNA-seq, in vitro functional assay and further metabolomic analysis, and hope to extend the study to other patients.

If this hypothesis is correct then it also raises the possibility of treatment. The commonly used and well tolerated bile acid ursodeoxycholic acid (UDCA) is a selective AKR1C2 inhibitor that can penetrate the brain [45, 46], and so may address the root cause of excessive inhibitory neurosteroid synthesis, while bioidentical iso (Sepranolone) and the GABA_AR modulating steroid antagonist GR3027 (Golexanolone) are showing promise in clinical trials for a variety of conditions associated with relative allo excess [27].

It seems extremely fortunate that in this case there was an error that could be detected by the short-read alignment, since this inversion would otherwise have been missed. Perhaps it was all part of the same event, with the DNA break being repaired by gene conversion but with the homology leading to the central section being fully inverted with the copy-and-paste errors further out. With such high homology making errors more likely, however, it is impossible to be sure that it was not sheer coincidence.

1) Long-read sequencing has the power to identify and fully characterise structural variants in difficult-to-align regions, although the challenges that these regions present to algorithmic analysis mean that manual oversight of the bioinformatic pipeline is currently still needed to ensure that variants are accurately identified. This has the potential to provide a molecular diagnosis for many more patients with known genetic conditions, and to reveal the underlying causes of conditions of unknown aetiology. We hope that as costs fall and algorithms improve, long-read methodologies will be widely adopted in clinical practice and research.

2) Genetic variants causing upregulation of AKR1C2 combined with partial loss of AKR1C1 may underlie some cases of unexplained chronic fatigue and possibly other conditions too. If they do then the biomarkers discovered in this study might allow better targeting of a specific patient subgroup for research and treatment. Further work is needed to determine if this is the case.

Ethics approval and consent to participate

This study was conducted in accordance with the Declaration of Helsinki. Written informed consent was obtained from all participants.

Consent for publication

Additional written informed consent was obtained from the patient for publication of this paper.

Availability of data and materials

Details of the large inversion are deposited in the European Variation Archive (EVA) with accession number ERZ21343151.

The datasets generated during the current study are not publicly available to protect the privacy of the patient but are available from the corresponding author on reasonable request.

Competing interests

The authors declare that they have no competing interests.

Funding

Long read sequencing was funded by a direct award to Genomics England from the UK Department of Health.

GC-MS/MS analysis was supported by the Ministry of Health of the Czech Republic - RVO (Institute of Endocrinology - EU, 00023761)

Authors’ contributions

JO conceived the hypothesis, performed the bioinformatic analysis and analysed the steroid results, with contributions from AG, MT and GE for the nanopore sequencing and MH for the GC-MS/MS and steroid analysis. JO drafted the manuscript, which GE critically edited. All authors read and approved the final manuscript.

Acknowledgements

This research was made possible through access to data in the National Genomic Research Library, which is managed by Genomics England Limited (a wholly owned company of the Department of Health and Social Care). The National Genomic Research Library holds data provided by patients and collected by the NHS as part of their care and data collected as part of their participation in research. The National Genomic Research Library is funded by the National Institute for Health Research and NHS England. The Wellcome Trust, Cancer Research UK and the Medical Research Council have also funded research infrastructure.

The authors would like to thank Dr Karel Valeš, Dr Filip Španiel and Dr Petra Fürstová for sharing their follicular control data, collected and assayed with funding from the Czech National Institute of Mental Health (NUDZ).

The authors would also like to thank Dr Norman Taylor for his advice on the manuscript.

100kGP 100,000 Genomes Project

17-OHP 17α-Hydroxyprogesterone

17,20α-diOHP 17α,20α-Dihydroxy-4-pregnen-3-one

20α-DHP 20α-Dihydroprogesterone

ADH Antidiuretic hormone

Adiol 3α-Androstanediol

Allo Allopregnanolone

CCDS Consensus coding sequence

CF Chronic fatigue

ME/CFS Myalgic encephalomyelitis/chronic fatigue syndrome

CRH Corticotropin releasing hormone

DHT 5α-Dihydrotestosterone

GABA_AR Gamma-aminobutyric acid type A receptor

gDNA genomic DNA

GWAS Genome-wide association studies

IGV Integrative Genomics Viewer

Iso Isopregnanolone

ONT Oxford Nanopore Technology

P Progesterone

PacBio Pacific Biosciences

Preg Pregnanolone

SNV Single nucleotide variant

SV Structural variant

THDOC Tetrahydrodeoxycorticosterone

UDCA Ursodeoxycholic acid

WGS Whole genome sequence

1. Dibble JJ, McGrath SJ, Ponting CP (2020) Genetic risk factors of ME/CFS: a critical review. Hum Mol Genet 29:R117–R124

2. Das S, Taylor K, Kozubek J, Sardell J, Gardner S (2022) Genetic risk factors for ME/CFS identified using combinatorial analysis. J Transl Med 20:598

3. Montoya JG, Dowell TG, Mooney AE, Dimmock ME, Chu L (2021) Caring for the Patient with Severe or Very Severe Myalgic Encephalomyelitis/Chronic Fatigue Syndrome. Healthcare 9:1331

4. Liu Z, Roberts R, Mercer TR, Xu J, Sedlazeck FJ, Tong W (2022) Towards accurate and reliable resolution of structural variants for clinical diagnosis. Genome Biol 23:68

5. Penning TM, Wangtrakuldee P, Auchus RJ (2019) Structural and Functional Biology of Aldo-Keto Reductase Steroid-Transforming Enzymes. Endocr Rev 40:447–475

6. Rogawski MA, Loya CM, Reddy K, Zolkowska D, Lossin C (2013) Neuroactive steroids for the treatment of status epilepticus. Epilepsia 54:93–98

7. Detlefsen AJ, Paulukinas RD, Penning TM (2023) Germline Mutations in Steroid Metabolizing Enzymes: A Focus on Steroid Transforming Aldo-Keto Reductases. Int J Mol Sci 24:1873

8. Mastrorosa FK, Miller DE, Eichler EE (2023) Applications of long-read sequencing to Mendelian genetics. Genome Med 15:42

9. Chen X, Schulz-Trieglaff O, Shaw R, Barnes B, Schlesinger F, Källberg M, Cox AJ, Kruglyak S, Saunders CT (2016) Manta: rapid detection of structural variants and indels for germline and cancer sequencing applications. Bioinformatics 32:1220–1222

10. Roller E, Ivakhno S, Lee S, Royce T, Tanner S (2016) Canvas: versatile and scalable detection of copy number variants. Bioinformatics 32:2375–2377

11. Li H (2018) Minimap2: pairwise alignment for nucleotide sequences. Bioinformatics 34:3094–3100

12. De Coster W, Rademakers R (2023) NanoPack2: population-scale evaluation of long-read sequencing data. Bioinformatics 39:btad311

13. Sedlazeck FJ, Rescheneder P, Smolka M, Fang H, Nattestad M, Von Haeseler A, Schatz MC (2018) Accurate detection of complex structural variations using single-molecule sequencing. Nat Methods 15:461–468

14. Pais LS, Snow H, Weisburd B, et al (2022) seqr : A web‐based analysis and collaboration tool for rare disease genomics. Hum Mutat humu.24366

15. Thorvaldsdottir H, Robinson JT, Mesirov JP (2013) Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration. Brief Bioinform 14:178–192

16. Danecek P, Bonfield JK, Liddle J, et al (2021) Twelve years of SAMtools and BCFtools. GigaScience 10:giab008

17. Kolmogorov M, Yuan J, Lin Y, Pevzner PA (2019) Assembly of long, error-prone reads using repeat graphs. Nat Biotechnol 37:540–546

18. Vaser R, Sović I, Nagarajan N, Šikić M (2017) Fast and accurate de novo genome assembly from long uncorrected reads. Genome Res 27:737–746

19. Shumate A, Salzberg SL (2021) Liftoff: accurate mapping of gene annotations. Bioinformatics 37:1639–1643

20. Hill M, Hána V, Velíková M, et al (2019) A method for determination of one hundred endogenous steroids in human serum by gas chromatography-tandem mass spectrometry. Physiol Res 179–207

21. Kulakovskiy IV, Vorontsov IE, Yevshin IS, et al (2018) HOCOMOCO: towards a complete collection of transcription factor binding models for human and mouse via large-scale ChIP-Seq analysis. Nucleic Acids Res 46:D252–D259

22. Kent WJ, Sugnet CW, Furey TS, Roskin KM, Pringle TH, Zahler AM, Haussler AD (2002) The Human Genome Browser at UCSC. Genome Res 12:996–1006

23. Ostinelli G, Vijay J, Vohl M-C, Grundberg E, Tchernof A (2021) AKR1C2 and AKR1C3 expression in adipose tissue: Association with body fat distribution and regulatory variants. Mol Cell Endocrinol 527:111220

24. Penning TM, Chen M, Jin Y (2015) Promiscuity and diversity in 3-ketosteroid reductases. J Steroid Biochem Mol Biol 151:93–101

25. Ong J, Kerr DIB (1990) GABA-receptors in peripheral tissues. Life Sci 46:1489–1501

26. Liang JJ, Rasmusson AM (2018) Overview of the Molecular Steps in Steroidogenesis of the GABAergic Neurosteroids Allopregnanolone and Pregnanolone. Chronic Stress 2:247054701881855

27. Bäckström T, Das R, Bixo M (2022) Positive GABA _A receptor modulating steroids and their antagonists: Implications for clinical treatments. J Neuroendocrinol. https://doi.org/10.1111/jne.13013

28. Steckelbroeck S, Jin Y, Gopishetty S, Oyesanmi B, Penning TM (2004) Human Cytosolic 3α-Hydroxysteroid Dehydrogenases of the Aldo-keto Reductase Superfamily Display Significant 3β-Hydroxysteroid Dehydrogenase Activity. J Biol Chem 279:10784–10795

29. Bixo M, Andersson A, Winblad B, Purdy RH, Bäckström T (1997) Progesterone, 5α-pregnane-3,20-dione and 3α-hydroxy-5α-pregnane-20-one in specific regions of the human female brain in different endocrine states. Brain Res 764:173–178

30. Kancheva R, Hill M, Cibula D, Včeláková H, Kancheva L, Vrbíková J, Fait T, Pařízek A, Stárka L (2007) Relationships of circulating pregnanolone isomers and their polar conjugates to the status of sex, menstrual cycle, and pregnancy. J Endocrinol 195:67–78

31. Mody I, Maguire J (2012) The reciprocal regulation of stress hormones and GABAA receptors. Front Cell Neurosci. https://doi.org/10.3389/fncel.2012.00004

32. Korpal AK, Han SY, Schwenke DO, Brown CH (2018) A switch from GABA inhibition to excitation of vasopressin neurones exacerbates the development of angiotensin II-dependent hypertension. J Neuroendocrinol 30:e12564

33. Tissue expression of AKR1C1 - Summary - The Human Protein Atlas. https://www.proteinatlas.org/ENSG00000187134-AKR1C1/tissue. Accessed 11 Feb 2023

34. Tissue expression of AKR1C2 - Summary - The Human Protein Atlas. https://www.proteinatlas.org/ENSG00000151632-AKR1C2/tissue. Accessed 11 Feb 2023

35. Quinkler M, Bumke-Vogt C, Meyer B, Bähr V, Oelkers W, Diederich S (2003) The Human Kidney Is a Progesterone-Metabolizing and Androgen-Producing Organ. J Clin Endocrinol Metab 88:2803–2809

36. Salyer SA, Parks J, Barati MT, Lederer ED, Clark BJ, Klein JD, Khundmiri SJ (2013) Aldosterone regulates Na+, K+ ATPase activity in human renal proximal tubule cells through mineralocorticoid receptor. Biochim Biophys Acta BBA - Mol Cell Res 1833:2143–2152

37. Auteri M, Zizzo M, Serio R (2015) The GABAergic System and the Gastrointestinal Physiopathology. Curr Pharm Des 21:4996–5016

38. Natelson BH (2019) Myalgic Encephalomyelitis/Chronic Fatigue Syndrome and Fibromyalgia: Definitions, Similarities and Differences. Clin Ther 41:612–618

39. Lin Y-J, Ko Y-C, Chow L-H, Hsiao F-J, Liu H-Y, Wang P-N, Chen W-T (2021) Salivary cortisol is associated with cognitive changes in patients with fibromyalgia. Sci Rep 11:1311

40. Komaroff AL, Lipkin WI (2023) ME/CFS and Long COVID share similar symptoms and biological abnormalities: road map to the literature. Front. Med. 10:

41. GABAA receptors modulate sympathetic vasomotor outflow and the pressor response to skeletal muscle metaboreflex activation in humans. https://doi.org/10.1113/JP277929

42. Hamidovic A, Davis J, Soumare F, Naveed A, Ghani Y, Semiz S, Khalil D, Wardle M (2023) Allopregnanolone Is Associated with a Stress-Induced Reduction of Heart Rate Variability in Premenstrual Dysphoric Disorder. J Clin Med 12:1553

43. Pollack B, von Saltza E, McCorkell L, Santos L, Hultman A, Cohen AK, Soares L (2023) Female reproductive health impacts of Long COVID and associated illnesses including ME/CFS, POTS, and connective tissue disorders: a literature review. Front Rehabil Sci 4:1122673

44. Arout CA, Sofuoglu M, Bastian LA, Rosenheck RA (2018) Gender Differences in the Prevalence of Fibromyalgia and in Concomitant Medical and Psychiatric Disorders: A National Veterans Health Administration Study. J Womens Health 27:1035–1044

45. Hara A, Matsuura K, Tamada Y, Sato K, Miyabe Y, Deyashiki Y, Ishida N (1996) Relationship of human liver dihydrodiol dehydrogenases to hepatic bile-acid-binding protein and an oxidoreductase of human colon cells. Biochem J 313:373–376

46. Parry GJ, Rodrigues CMP, Aranha MM, Hilbert SJ, Davey C, Kelkar P, Low WC, Steer CJ (2010) Safety, Tolerability, and Cerebrospinal Fluid Penetration of Ursodeoxycholic Acid in Patients With Amyotrophic Lateral Sclerosis. Clin Neuropharmacol 33:17–21

No competing interests reported.

Download PDF

Journal Publication

published 16 Nov, 2023

Read the published version in Journal of Translational Medicine →

Version 2

posted

You are reading this older preprint version

Read the latest preprint version →

Long read sequencing characterises a novel structural variant, revealing underactive AKR1C1 with overactive AKR1C2 as a possible cause of unexplained severe fatigue

Status:

Journal Publication

Version 2

Abstract

Figures

Background

Materials And Methods

Clinical Details

Short-read WGS

Long Read WGS

Bioinformatic analysis

Serum steroid analysis

Results

Short-read WGS indicates a large inversion around AKR1C1 and AKR1C2, with a splice-site variant in the 5’UTR of AKR1C2

Characterising the variant from the long read data

Initial overview

Confirming the inversion

Homology hinders phasing

Additional insertion near AKR1C2

Assembling the alignments

Comparison with short-read sequence

Serum steroid analysis suggests a functional effect with increased AKR1C2 and reduced AKR1C1 activity

Neurosteroid homeostasis

Precursor-product ratios of AKR1C1 (20α-reduction)

Precursor-product ratio of AKR1C2 (3α-reduction)

Discussion

Conclusions

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Competing interests

Funding

Authors’ contributions

Acknowledgements

Abbreviations

References

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 2