Molecular dynamics simulations of cytochrome P450 aromatases reveal structural variances across the cat family

Molecular dynamics (MD) simulation method

Fully-atomistic molecular dynamics simulations were conducted using Nanoscale Molecular Dynamics software (NAMD 2.14) (Phillips et al. 2005) with the CHARMM36 force field (Huang and MacKerell 2013), and the systems were solvated using the CHARMM-modified TIP3P water model. Visual Molecular Dynamics Software (VMD) was used for visualisation and analysis (Humphrey et al. 1996). Energy minimisation was performed using the conjugate gradient algorithm. Constraints applied to the protein backbone were gradually removed over a span of 2.04 ns. The temperature of the system was then increased gradually from 0 K to 310 K over 0.62 ns. This was followed by NPT equilibration for 1 ns, and then further pressure equilibration under the NPT ensemble at 310 K and 1 atm for 10 ns (NPT refers to the isothermal-isobaric ensemble, where the number of particles (N), pressure (P), and temperature (T) are held constant during the simulation). For membrane-embedded systems, standard CHARMM-GUI (Jo et al. 2008) protocols were followed for system construction and equilibration. After equilibration, production runs were performed without restraints: 100 ns for monomeric systems, 150 ns for dimeric systems and 350 ns for membrane systems. Further details on the methodology are provided in Supplementary Data S1.

Amino acid sequences and alignment

The seven cat family aromatase sequences in the study include: Homotherium latidens (homotherium) (Barnett et al. 2020), Panthera leo (lion) XP_042797965.1, Panthera tigris (tiger) XP_042844798.1, Prionailurus bengalensis (leopard cat) XP_043411365.1, Puma concolor (puma) XP_025775219.1, Acinonyx jubatus (cheetah) XP_053079996.1, and Felis catus (domestic cat) XP_044915259.1, which were compared to that of human aromatase. The amino acid sequence in the case of Homotherium latidens was computationally derived from the DNA sequence provided by Prof. Michael Westbury from the University of Copenhagen (Barnett et al. 2020; Westbury et al. 2021). Felis catus aromatase sequence (XM_006932525.3) was used as the reference genome to extract the approximate region of Homotherium latidens CYP19A1 gene and manually mapped for exon identification. The exons were spliced and translated using the Expert Protein Analysis System (ExPASy) translate tool (Artimo et al. 2012), which provided the final amino acid sequence of Homotherium latidens aromatase. The remaining six cat family aromatase sequences (listed in Table 1) were obtained from the NCBI database (Sayers et al. 2022).

Clustal Omega in default mode (Sievers et al. 2011) was used to perform sequence alignments on these seven cat family aromatase species and human aromatase (Protein Data Bank (PDB) ID: 3EQM) (Ghosh et al. 2010) and was visualised using the BioEdit Software (Hall 1999), as shown in Fig. 1. The aromatase sequences within the cat family exhibited approximately 99% conservation, whereas sequence identity decreased to around 86% if compared to human aromatase (Supplementary Information; Table S1). Position-specific conservation percentages were calculated, revealing high conservation (> 95%) of key functional regions, including the haem-binding site, cytochrome P450 oxidoreductase (CPR) binding residues, and six substrate recognition sites (SRSs) (Supplementary Information; Table S2 and Fig. S3).

Fig. 1.

Multiple sequence alignment (Edgar and Batzoglou 2006) of aromatase sequences from Homo sapiens (PDB: 3EQM), Homotherium latidens (Westbury et al. 2021), Panthera leo (NCBI: XP_042797965.1), Panthera tigris (NCBI: XP_042844798.1), Prionailurus bengalensis (NCBI: XP_043411365.1), Puma concolor (NCBI: XP_025775219.1), Acinonyx jubatus (NCBI: XP_053079996.1), Felis catus (NCBI: XP_044915259.1), aligned using Clustal Omega (Sievers et al. 2011) and compared in Bio-edit software (Hall 1999). The residues 1–54 belong to the transmembrane (TM) region, α helices are indicated and labelled in blue boxes from A–L, four β strands are labelled in red boxes, while K′–L loop and haem-binding regions are indicated by green boxes. The labelling of secondary structure elements is according to the X-ray crystal structure study by Ghosh et al. (2010) and visually represented in Fig. S2.

Homology modelling

Sequence alignments revealed that the aromatases of lion, tiger, and leopard cat were identical (Table S1); therefore, the tiger sequence was selected as representative structure for subsequent analysis (as mentioned in Table 1). The 3D structures of the remaining five cat family aromatases were constructed using the homology modelling tool SWISS-MODEL (Waterhouse et al. 2018), with the X-ray crystal structure of human aromatase (PDB ID: 3EQM (Ghosh et al. 2010), Fig. S2) serving as the template. The quality of the homology-modelled structures was evaluated using the Quality Assessment module of the SWISS-MODEL server (Waterhouse et al. 2018) (Table S3). Metrics such as QMEAN Z-scores (Benkert et al. 2008), low molprobity values (~2.12) (Chen et al. 2010) and > 94% of residues in favoured regions of the Ramachandran plot indicated that the modelled structures were of high quality. However, further refinement, particularly around proline and glycine residues, was anticipated to be resolved during the relaxation process of molecular dynamics simulations.

Molecular dynamics simulation systems

Aromatase monomers for five cat family members (Table 1) and human were constructed both with and without haem. A previously reported (Martin et al. 2015) dimeric structure of human aromatase was used as a template for the construction of Homotherium and puma aromatase dimers, representing the extinct and extant species of cat family. All simulations were conducted in the absence of androstenedione (ASD) and in the resting state, as this study focuses solely on observing dimeric interactions. Interactions in the presence of the substrate are beyond the scope of this work and are currently being investigated in our lab. Further details on the construction of monomeric and dimeric systems, both in solution and in the membrane, can be found in the Supplementary Data S1.

Sequence alignment

The amino acid sequence alignments, shown in Fig. 1, for the human and cat family aromatases produced a sequence identity of > 99% within the cat family and ~86% identity between cats and the human sequence (Table S1). This high degree of conservation within the cat family was expected for an essential reproductive enzyme, further reinforcing that the structure and function of aromatases are preserved across all feline species included in this study. Specifically, positively charged residues on the enzyme’s surface, such as ¹⁰⁸K, ¹⁴²K, ¹⁴⁵R, ¹⁵⁰K, ³⁵²K, ³⁵⁴K, ⁴²⁰K, and ⁴²⁵R (Fatima et al. 2020) are completely conserved in both humans and within the cat family members investigated. In the human aromatase, these residues are implicated in binding to the redox partner, cytochrome P450 oxidoreductase (CPR), promoting electron transfer from NADPH to the haem centre (Table S2) (Ghosh et al. 2010; Gilardi and Di Nardo 2017; Di Nardo et al. 2021). Active site residues, including ³⁰⁶A, ³⁰⁹D, ³¹⁰T, ⁴⁷⁸S, and ⁴⁸⁰H (Fig. 1) (Ghosh 2023), showed 100% conservation, indicating that the active site is similar for the cat family members and human aromatase (Ghosh et al. 2010; Di Nardo et al. 2021). Additionally, with over 95% conservation observed in the substrate-binding sites (Table S2), it is expected that the substrate, androstenedione (ASD), would participate in the catalytic cycle of cat family aromatase in a manner similar to its role in humans. This high level of conservation supports analogous CPR-aromatase and substrate-aromatase interactions in the cat family and in humans. A non-conserved non-polar to polar alteration of ⁴³⁶G (human) → ⁴³⁶S (cat family) was observed in the haem-binding domain, which being adjacent to the iron coordinating cysteine (CYS) residue, may slightly affect the binding affinity for the haem (Fatima et al. 2020). Within the cat family, sequence comparisons revealed that the extinct genus Homotherium shares over 99% sequence similarity with other feline members. Differences were identified between Homotherium and extant species, such as the puma, that includes ¹⁰⁴F (Homotherium) → ¹⁰⁴V (puma), ¹⁶⁸V (Homotherium) → ¹⁶⁸I (puma), and ³³⁹I (Homotherium) → ³³⁹V (puma). Key structural motifs, including the four-helix bundle (helices D, E, I, and L) and substrate recognition sites (SRSs), have been preserved, highlighting strong selective pressures to maintain enzymatic specificity and efficiency (Tiwary and Li 2009; Martin et al. 2015; Di Nardo et al. 2021; Zhang et al. 2022). These findings reflect the remarkable stability of aromatase sequences over 15 million years, while leaving room for ongoing evolutionary refinement that would likely require even longer evolutionary timescales.

Sequence comparisons also revealed that the 14% difference between the cat family and human aromatase sequences stems from 69 differing residues (Table 2). Of these, 45 amino acid residues represent conserved chemical properties (Table 2a), thus are likely to be important for maintaining aromatase structure, activity and/or function. The remaining 24 (Table 2b, c) are non-conserved changes. Most of these non-conserved residues were amino acid residues that were charged in humans which were replaced by either hydrophobic or polar residues in cats (Table 2c). To examine the structural orientation of these charged residues in humans that are not conserved in cats, the 3D structures of Homotherium, cheetah, tiger, and domestic cat were generated using the X-ray structure of the human aromatase (Ghosh et al. 2010) as a template for homology modelling, as described in the Materials and methods section. Further, molecular dynamics simulations were carried out to achieve equilibrated structures.

MD simulations of aromatase monomer

Stability of monomeric aromatases

Fig. 2a, b plot the root mean square deviations (RMSD) of the protein backbone relative to the protein’s initial conformation, for all six aromatase species over a 100 ns simulation period. These simulations were conducted both without and with haem cofactor, respectively, to evaluate the stability of the equilibrated structures. Notably, the N-terminal region (residues 1–44) was excluded from these simulations, as the crystal structure of human aromatase lacks the N-terminal segment and includes only the globular domain (Ghosh et al. 2010). The 3D models incorporating haem were constructed with coordination bonds, where iron (Fe) was coordinated to both the axial water molecule and the cysteine residue (⁴³⁷C) to represent the resting state of haem (Fig. S4). All aromatase structures without haem exhibited higher RMSD values (~3 ångström (Å), Fig. 2a) compared to those with haem (RMSD ~2.5 Å, Fig. 2b). Additionally, RMSF (root mean square fluctuations) values (Fig. 2c, d) showed higher fluctuations in the β4–5 region, K′–L loop residues (400–460), and haem-binding residues in absence of haem. Hence, the RMSD and RMSF results suggest that haem plays a key role in stabilising the overall 3D structure. The lower RMSD values in the haem-bound systems are likely due to stabilising interactions between aromatase residues ¹¹⁵R, ¹⁴¹W, ³⁷⁵R, and ⁴³⁸A with the haem molecule (as shown in Fig. S7).

Fig. 2.

Panels (a, b) and (c, d) show the RMSD and RMSF plots, respectively, for six species: human (green), puma (magenta), Homotherium latidens (orange), cat (blue), cheetah (red), and tiger (yellow). Left panels correspond to simulations without haem, and right panels include haem. For each species, the bold lines represent the average across three replicate simulations, with lighter shades denoting individual replicates. These simulations were conducted in water as the solvent and ionised with 0.15 M NaCl forming a simulation box of 145 Å × 145 Å × 145 Å. The residues highlighted show a large difference in fluctuations in the absence/presence of haem, indicated as β4–5, K′–L loop and haem-binding regions.

Surface electrostatic potentials of monomer aromatases

Analysis of the equilibrated haem-containing aromatase structures revealed that the non-conserved charged residues in the cat family aromatases were located on the G-helix, B-helix, C-helix, and β-loop 6–7 regions (refer to Fig. S5). Of all the non-conserved residues, approximately 33% were situated in a region that was earlier proposed as a potential dimeric interface of the human aromatase structure by Fatima et al. (Martin et al. 2015; Fatima et al. 2020) (Fig. 3a–c). It was hypothesised that homodimerisation could potentially enhance the catalytic efficiency of androstenedione conversion to oestrone (Martin et al. 2015; Holien et al. 2017; Fatima et al. 2020). Their analysis suggested that dimerisation occurs through favourable electrostatic interactions (a distinct pattern of alternating positive and negative charges) at the proposed interface between the two monomers (Fatima et al. 2020). Therefore, the distribution of charged residues on the proposed dimerisation interface and the resulting surface electrostatic potential were analysed for human and cat family aromatases. Fig. 3a–c show that all non-conserved charged residues, namely ¹³⁹E_h → T_c, ¹⁴³T_h → A_c, ¹⁴⁴T_h → V_c, ³⁷⁹E_h → Q_c, ⁹⁰S_h → C_c, ¹³⁰K_h → N_c, ²⁶⁰A_h → E_c, and ²⁶⁹E_h → A_c, were located on the proposed dimeric interface of these enzymes (subscripts h and c denote human and cat family species, respectively, and superscripts indicate the positions of non-conserved amino acid residues). For ease of comparison, only the interfaces of human, puma, and Homotherium aromatases are shown in Fig. 3, as most cat family species share similar residues. Puma and Homotherium represent extant and extinct species, respectively, with human included as a reference. Electrostatic surfaces for the rest of the cat family species are provided in Fig. S8. The absence of ¹³⁰K_h and presence of ²⁶⁰E_c reduces the overall electrostatic potential at the centre of the proposed dimerisation interface in cat family aromatases, which is positively charged in human aromatase (Fig. 3d–f). Thus, the absence of charged residues on the surface of cat family aromatases ³⁷⁹E, ¹³⁹E, ¹³⁰K and ²⁶⁹E revealed subtle differences between the surface electrostatic properties of cat family and human aromatase.

Fig. 3.

Panels (a), (b), and (c) illustrate the structural models of aromatases from human (green), puma (magenta), and Homotherium (yellow) highlighting the differences in the residues located at the dimerisation interface. The van der Waals (vdW) residues are colour-coded as: red for negative amino acids, blue for positive, ice-blue for polar, and grey for hydrophobic residues. Panels (d), (e), and (f) present the averaged electrostatic surfaces of human, puma and Homotherium aromatases, particularly focusing on the dimeric interface (Martin et al. 2015; Fatima et al. 2020), generated using Visual Molecular Dynamics (VMD). Additional details on the electrostatic surfaces of other cat family species are provided in Supplementary Information, Fig. S8.

Although distinct stripes of opposite charges were observed in human aromatase (Fig. 3d), they appeared more diffused in cat family counterparts (Fig. 3e, f). This suggests reduced stability of cat family aromatase homodimers due to a potential weakening of monomer-monomer interactions at the proposed dimerisation interface. To investigate this further, homodimer structures of human, puma, and Homotherium aromatases were constructed (refer to Table 1), using the known human aromatase homodimer as a reference structure. Additional details on the construction of aromatase homodimers in solution and within membranes can be found in sections S1.2, S1.3 of the Supplementary Information.

MD simulations of aromatase dimers

Fig. 4 compares the relative stability of human, puma and Homotherium homodimers in solution (Fig. 4a–c) and in membrane (Fig. 4d–f), respectively. The homodimer simulations in solution excluded the N-terminal regions (residues 1–44) to focus on the globular domains of the dimers in solution (Fig. S6b). The N-terminal region corresponds to the transmembrane domain of the aromatase when embedded in a membrane. Fig. 4a–c shows that the human dimer in solution exhibited lower RMSD values (~4 Å) compared to the cat family dimers in solution (~5 Å), suggesting stronger interfacial interactions in the case of the human homodimer. This is consistent with the proposed weakening of electrostatic interactions across the proposed dimerisation interface in the cat family aromatases (Fig. 3). Further, to mimic their natural state, membrane-embedded dimeric systems were also constructed for the human, puma, and Homotherium homodimers (Fig. S6c). Human aromatase homodimers exhibited significantly lower RMSD values (Fig. 4d) compared to both cat family homodimers (Fig. 4e, f), indicating greater stabilisation of the human homodimer in the membrane. Notably, the difference in stability between the human and cat family homodimers is more pronounced in membrane than in solution. The RMSF values (Fig. S10) also showed higher fluctuations in the I, J, and K helices and beta strands 1–2 and 4–5 of Homotherium and puma compared with the corresponding monomeric fluctuations (Fig. 2b). This suggests that dimer formation influences residue-residue interactions at the interface, which in turn affects fluctuations of the connected J and K helices on the outer end of interfaces. Since the I and K helices form the substrate-binding sites (SBS 4 and 5; Table S2) (Di Nardo et al. 2021), dimer formation could influence substrate interaction and activity towards aromatase. Further to understand the dimeric interactions, hydrogen bonding and salt bridge formation between interacting residues at the dimerisation interface were monitored throughout the simulation (Fig. S11–S14).

Fig. 4.

Panels (a–f) represent RMSD profiles of dimeric aromatase systems for human (light and dark green), puma (pink and magenta), and Homotherium latidens (orange and brown) simulated over 150 ns in solution and for 350 ns within a membrane. Bold lines indicate average RMSD values, whereas dotted lines represent individual chain RMSD values (a and b) for the trajectory sets T1 and T2.

Aromatase dimeric interactions

In solution, the human aromatase dimer initially formed numerous intermolecular hydrogen bonds and salt bridges (Figs S11, 12). However, towards the end of the simulations, only a few residue pairs interacted strongly: (a) ARG293_A⋯GLY289_B and oppositely charged residue pairs (b) GLU139_A⋯LYS130_B, (c) LYS130_A⋯GLU139_B, and (d) LYS390_A⋯GLU92_B (Fig. S11a, c). Here, subscripts A and B denote chain A and chain B of the homodimers. These residue pairs maintained ~3 Å distances, stabilising the dimer through hydrogen bonds and electrostatic interactions. In the puma dimer, LYS293_A⋯GLY289_B and LYS390_A⋯GLU92_B (Fig. S11a, b) formed hydrogen bonding and salt-bridge interactions, whereas THR129_A⋯GLU139_B (Fig. S11d) failed to interact due to repulsion from the negative centre, preventing the formation of stable salt bridges throughout the simulation. In Homotherium, the oppositely charged residue pair GLU92_A⋯LYS390_B (Fig. S11a, b, e) formed salt bridges and hydrogen bonds. Non-conserved mutations (¹³⁰K_h → N_c and ²⁶⁰A_h → E_c in the cat family) resulted in a negative potential at the centre of the dimerisation interface, initially causing repulsion between residue pairs. Over time, the monomers adjusted, resulting in strong interactions between residue pairs, THR139_A⋯GLU129_B and LYS293_A⋯GLY289_B (H-bonding, Fig. S11a, b, e) for the remainder of the simulation. However, in some cases, snapshots (Fig. S9) revealed that the Homotherium homodimer in solution dissociated toward the end of the 150 ns simulation, increasing inter-dimer distance and reducing interactions. This suggests that a membrane could serve as a natural constraint at the N-terminal region, stabilising dimeric interactions. Therefore, these dimeric systems were embedded within a lipid bilayer to assess its effect on dimer stability.

When embedded in the membrane, the human dimer exhibited greater stability, characterised by stronger hydrogen bonding and salt-bridge interactions (Fig. 5, Figs S13, S14). Specifically, residue pairs (1) ARG293_A⋯GLY289_B showed strong H-bonding interactions, while the rest of the oppositely charged residue pairs (2) GLU139_A⋯LYS130_B, (3) LYS389_A⋯ASP232_B, (4) LYS390_A⋯ASP232_B, and (5) LYS236_A⋯GLU379_B maintained ~3 Å distance between them, forming a constant hydrogen bond and salt bridges along with electrostatic interactions between them observed in both the trajectory sets T1 and T2 (Fig. 5a). In case of puma and Homotherium, the interacting residues varied between the two trajectory sets (Fig. 5b–c, Figs S13, S14) due to weaker interactions at the dimeric interface caused by the absence of key charged residues (³⁷⁹E, ¹³⁹E, ¹³⁰K, and ²⁶⁹E) in the cat family. Fig. 6 illustrates that human aromatase dimers form a significantly higher number of stabilising interactions, hydrogen bonds, salt bridges, and electrostatic forces compared to the cat family aromatases, both in solution and membrane environments. The diminished stability and fewer interacting pairs observed in puma and Homotherium suggest evolutionary differences in dimerisation stability and interaction patterns. These differences likely stem from the absence of key charged residues (³⁷⁹E, ¹³⁹E, ¹³⁰K, and ²⁶⁹E) on the surface of cat family aromatases, which, in humans, contribute to critical hydrogen bonding and electrostatic interactions. The reduced propensity for stable dimer formation in cat family aromatases highlight a potential evolutionary divergence, where structural adaptations may have influenced functional dynamics across species.

Fig. 5.

Panels (a–c) represent hydrogen bonding distances tracked throughout 350 ns simulation for dimer structures of human (green), puma (magenta), and Homotherium latidens (yellow) in membrane environments, with 3.5 Å as H-bonding cut-off distance. Structural representations in panels (d–f) highlights key amino acid residues hydrogen bonded in T1 trajectory simulations (detailed information on hydrogen bonding and salt-bridge interactions on trajectories T1 and T2 is provided in the Fig. S13–S15).

Fig. 6.

Panels (a) and (b), show the average number of hydrogen bonding (blue) and salt bridges (red) forming amino acid residue pairs between the human, puma and Homotherium homodimers when in solution (a) and when embedded in the membrane (b), tracked over 150 ns and 350 ns simulation time respectively.

Stability of haem

The stability of haem in aromatase is likely to have a direct impact on the enzyme’s ability to carry out its catalytic function (Gilardi and Di Nardo 2017). Haem, when bound to aromatase, undergoes a series of redox reactions that are essential for the conversion of androgens to oestrogens. Hence, when it becomes unstable, it affects the redox reactions, resulting in a decrease in enzymatic activity, and consequently, a decrease in oestrogen production (Hong et al. 2007). Comparative analysis (Fig. 7) showed that haem interactions were significantly stronger in dimeric, membrane-embedded systems than in monomers or dimers in solution, with a 10-fold increase in interaction occupancy. Cases where occupancy exceeded 100% reflect multiple hydrogen bonds involving the same haem atom. These findings suggest that the membrane environment enhances aromatase stability and function.

Fig. 7.

Percentage occupancy of hydrogen bonding interactions between aromatase residues and haem, tracked over the 150 ns of simulation time, for monomeric and dimeric structures in solution, and 350 ns in membrane-embedded dimeric systems (stars represent data points, and bars the standard deviation).

Stability of aromatase as a dimer

Cytochrome P450 enzymes have frequently demonstrated dimer formation in crystals through the interaction of the helices F and G regions (Ghosh et al. 2010; Reed and Backes 2017). In 2015, evolutionary studies suggested that dimerisation of human aromatase increased the enzymatic activity and efficiency (Martin et al. 2015; Fatima et al. 2020). Our simulations show that human aromatase dimers are significantly more stable in membranes, whereas cat family dimers exhibit lower stability. To simplify understanding, the key interactions observed at the dimeric interface are categorised as top, middle and bottom interactions (refer Fig. 5d–f). In human aromatase, H-bonding interaction between the residues ARG293(+ve) and GLY289 occurred at the top region of dimeric interfaces. In the cat family, ARG(+ve) is substituted by LYS293(+ve), which is a conserved substitution and maintains similar interaction patterns. In the middle region, electrostatic interactions, such as those between residues LYS130(+ve) and GLU139(−ve) in human aromatase, contribute to dimer stability. However, in cat family aromatases, LYS130 is replaced by ASN130 (polar) and GLU139(−ve) is replaced by THR139, a polar residue, which weakens the electrostatic interactions. At the lower dimer interface, interactions involving (+ve)LYS390⋯ASP232(−ve) and (+ve)LYS390⋯GLU92(−ve) were stabilised in human aromatase due to the constraints provided by the membrane and helped in the formation of strong electrostatics, H-bonding and salt-bridge interactions. On the other hand, these interactions were weaker in the case of cat family aromatase due to the absence of charge residues on the interface (Fig. 3). Additionally, haem structure stability was higher in homodimeric forms embedded in membrane compared to monomeric forms, as evidenced from Fig. 7. This comparative difference in stability likely contributes to higher catalytic efficiency observed in human aromatase, compared with some species. Overall, human homodimer aromatase appeared to be a stable aromatase system when embedded in the membrane that increased haem stability, as shown in Fig. 7, while also exhibiting strong dimeric interactions (Figs 5a and 6b). Puma aromatase, on the other hand, achieves an intermediate level of dimeric interactions (Figs 5b and 6b) but loses stability as an aromatase system, as observed with lower haem stability (Fig. 7). In contrast, Homotherium aromatase has the fewest dimeric interactions (Figs 5c and 6b), forming the least stable homodimer. These suggests human forms a stable homodimer, puma forms a homodimer which is less stable compared to human and Homotherium does not favour dimerisation interactions compared to rest of the species. To assess the stability and re-association of the dimeric interface, a 100 ns test simulation was performed starting with the two monomers placed 10 Å apart (Fig. S16). Over the course of the simulation, the monomers gradually moved closer, resulting in the decrease of the centre-of-mass distance and formation of interfacial contacts (Fig. S16). Thus, leading to an increase in the number of hydrogen bonds and salt bridges formed at the dimeric interface (Fig. S17). A comparison with the 350 ns membrane-embedded dimer simulation shows that the key interacting residues are largely restored, indicating reformation of the native interface (Table S6).

While comparative biochemical data on purified cat aromatase is limited, one study examining aromatase active in cat testes (Meyer et al. 2018) reported activity levels of approximately 0.45 pmol/min.mg protein, which is markedly lower than values commonly reported for human placental aromatase, which range from 3–4 pmol/min.mg protein (Okubo et al. 1996). While the assays were conducted in different tissue contexts and conditions, this difference provides in vitro experimental support for our computational observation that reduced dimer stability in the cat family may contribute to lower catalytic efficiency. Our simulations show that the human homodimer forms more stable interfacial contacts compared to both puma and Homotherium dimers, particularly in membrane environments. Given our earlier in vitro studies implicating oligomerisation in steroidogenic cytochrome P450 enzymes (Praporski et al. 2009; Ghosh et al. 2011), we are confident that homodimerisation influences catalytic behaviour, however, we acknowledge that further biochemical validation is desirable.