Low abundance of NDUFV2 and NDUFS4 subunits of the hydrophilic complex I domain and VDAC1 predicts mammalian longevity.

Mitochondrial reactive oxygen species (ROS) production, specifically at complex I (Cx I), has been widely suggested to be one of the determinants of species longevity. The present study follows a comparative approach to analyse complex I in heart tissue from 8 mammalian species with a longevity ranging from 3.5 to 46 years. Gene expression and protein content of selected Cx I subunits were analysed using droplet digital PCR (ddPCR) and western blot, respectively. Our results demonstrate: 1) the existence of species-specific differences in gene expression and protein content of Cx I in relation to longevity; 2) the achievement of a longevity phenotype is associated with low protein abundance of subunits NDUFV2 and NDUFS4 from the matrix hydrophilic domain of Cx I; and 3) long-lived mammals show also lower levels of VDAC (voltage-dependent anion channel) amount. These differences could be associated with the lower mitochondrial ROS production and slower aging rate of long-lived animals and, unexpectedly, with a low content of the mitochondrial permeability transition pore in these species.

Introduction

Complex I (Cx I) (NADH-ubiquinone oxidoreductase; EC 1.6.5.3) is an electron entry point in the mitochondrial respiratory electron transport chain (ETC). Cx I catalyses NADH oxidation reducing ubiquinone to ubiquinol, importantly contributing to the proton motive force used to synthesize ATP by the oxidative phosphorylation [1]. Cx I also produce reactive oxygen species (ROS), initially superoxide radicals, which can damage all cellular components. Although at least 11 sites producing ROS have been identified, Cx I and complex III (Cx III) are conventionally recognized as the major sources of ROS at the ETC [2]. Mitochondrial ROS production (mitROSp) has been considered one of the main effectors responsible for aging and longevity [3,4].

Low rates of mitROSp have been described in many long‐lived mammalian and bird species [[3], [4], [5]]. These studies generally demonstrated the existence of a negative correlation between mitROSp and longevity. Among the two main ROS generating ETC complexes, the low ROS production of various long-lived species has been localized at Cx I [[6], [7], [8]]. Interestingly, different pro-longevity nutritional and pharmacological interventions like dietary (DR) and methionine restriction, and rapamycin treatment have been also associated with decreased mitROSp at Cx I [3,9]. Although is a matter of controversy, within complex I longevity-related ROS production has been observed in the hydrophilic domain [10], which could be due either to the flavin semiquinone, to some of the eight different FeS clusters of this domain [[11], [12], [13]] or both to flavin and FeS centers.

The underlying mechanism responsible, at least in part, for the low mitROSp of some long‐lived species has been attributed to a smaller degree of electronic reduction of Cx I under basal conditions [4,8]. Additional studies demonstrated that a lower Cx I content could also explain the low ROS production of those long-lived species [14,15]. Pro-longevity nutritional interventions also induced decreased Cx I content concomitantly with a lower mitROSp [10,[16], [17], [18], [19]].

Mammalian Cx I is the largest component of the ETC built of 45 different subunits in mammals [1]. The conserved L-shaped core, formed by 14 subunits sufficient for catalysis, is surrounded by 31 accessory subunits forming an interlinked shell with unclear function [20,21]. Among the 14 core subunits, the 7 mitochondrial-encoded ND subunits are present in the hydrophobic membrane domain, and the other 7 nuclear-encoded (NDUF) subunits are present in the hydrophilic matrix domain [22]. The 31 supernumerary NDUF accessory subunits are also nuclear coded [20]. However, it is totally unknown if some particular Cx I subunits, especially some NDUF subunits of the Cx I hydrophilic domain, could be involved in the determination of the longevity-related low complex I ROS production of long-lived animal species.

On the other hand, the mitochondrial permeability transition pore (mPTP) is a protein complex that permits diffusion of molecules of molecular mass up to 1500 Da across the mitochondrial inner membrane [23,24]. Many substrates can pass from mitochondrial matrix to cytosol through this pore. The mPTP opening is promoted through oxidation of lipids and proteins of mitochondrial membrane by ROS [23] and is related to pathological stages. However, some studies suggest that mPTP opening by mitROS can also drive the progression of aging [23,25]. It has been hypothesized that mPTP comprises inner and outer mitochondrial membrane proteins such as the ATP synthase, key in the channel formation, and a diversity of regulatory components including the voltage‐dependent anion channel (VDAC), the pro-death Bcl-2 family member proteins Bax and Bak, the translocator protein (TSPO), the hexokinase (HK), the adenine nucleotide translocase (ANT), the phosphate carrier (PiC), and the cyclophilin-D [[26], [27], [28]], all of which have been studied as candidates for exact molecular mPTP identity which remains uncertain. Recent studies have shown that decreased VDAC promotes longevity by decreasing mitochondrial permeability in C. elegans [29].

The purpose of this study was to investigate the abundance of Cx I hydrophilic domain selected subunits and VDAC as possibly related to the low mitROSp of long‐lived mammalian species. We used droplet digital PCR (ddPCR) and western blot methods to define the steady-state levels of gene expression and protein content of mitochondrial electron transport chain Cx I subunits and VDAC in heart tissue of eight mammalian species showing more than one order of magnitude of difference in longevity, from 3.5 years in mice to 46 years in horses. Cx I hydrophilic subunits were selected based on: i) previous partial proteomic studies of the ETC revealing the subunits more acutely decreasing in dietary restriction or varying in aged or longevity mutant mice [10]; ii) their core (NDUFV2 and NDUFS3) or supernumerary (NDUFS4 and NDUFS5) character; iii) their especial situation in the electron path of the hydrophilic domain (NDUFV2); or iv) their lack of variation in those longevity models to be potentially used as a reference (NDUFA9). Abundance of the VDAC component of the mPTP was also measured due to its possible relationship with mitROSp and aging.

Results

Multivariate statistics reveals a species-specific Cx I profile

Multivariate statistics was applied to determine whether specific Cx I subunits gene expression (ndufv2, ndufs3, ndufs4 and ndufs5) and protein content (NDUFV2, NDUFS3, NDUFS4 and NDUFA9) differ among mammals. Non-supervised principal component analysis (PCA) revealed the existence of a species-specific gene and protein profile of the Cx I, capable to explain up to 65.4% samples variability (Fig. 1A). A hierarchical clustering of the samples represented by a heat map revealed specific Cx I patterns for rodents (mouse, rat and guinea pig) and non-rodents (rabbit, pig, cow and horse) (Fig. 1B and C). These results were confirmed by performing a supervised analysis, such as partial least squares discriminant analysis (PLS-DA) (Fig. 1D). Cross-validation values of PLS-DA model scored a maximum Q2 of 0.3 and R2 of 0.4 when using only one component (Fig. 1E). Permutation tests (1000 repeats) yielded a low p = 0.02, indicating that none of the distributions formed by the permuted data was better than the observed statistic based on the original data (Fig. 1F). The inter-group discriminating power of the different measured features is evaluated by a variable importance projection (VIP) score > 1.5, which was obtained only for NDUFV2 (Fig. 1G). These results were supported by applying a Random Forest (RF) classification algorithm which revealed a species overall classification error of 0.26, being ndufv2 the variable with the highest contribution to classification accuracy (Fig. 1H and I).

Fig. 1

Multivariate statistics reveals a species-specific gene expression and protein content of Cx I. A) Principal component analyses (PCA) representation of gene expression and protein content of CI subunits. X: Principal component 1 (PC1); Y: Principal component 2 (PC2); Z: Principal component 3 (PC3). B) Hierarchical clustering of individual animal samples according to gene expression and protein content of CI. C) Hierarchical clustering of animal species according to average sample values of gene expression and protein content of CI. D) Partial least squares discriminant analysis (PLS-DA) representation of gene expression and protein content of CI. X: Component 1 (C1); Y: Component 2 (C2); Z: Component 3 (C3). E) Cross validation (CV) analyses (10-fold CV method) of the PLS-DA model. F) Permutation test (1000 repeats) using separation distance. G) Variable importance projection (VIP) scores indicating the elements which contribute the most to define the first component of a PLS-DA. H) Random Forest (RF) classification algorithm. I) VIP scores for RF.

The species-specific Cx I profile is associated with animal longevity

Gene expression and protein content of Cx I were also correlated with mammalian longevity. Specifically, long-lived animals showed decreased gene expression of ndufv2 and ndufs3 (Fig. 2A). Moreover, protein content of NDUFV2 and NDUFS4 were also decreased in long-lived mammals (Fig. 2B). Yet no correlation between protein content and gene expression was found (Supplementary Fig. 1).

Fig. 2

Cx I gene expression (A) and protein content (B) are linearly correlated with mammalian longevity. Pearson correlation was performed between maximum longevity (ML) and gene expression or protein content. Linear regression was applied when significant relationships were found. R2(ndufv2) = 0.48; R2(ndufs3) = 0.20; R2(NDUFV2) = 0.10; R2(NDUFS4) = 0.15. NDUFA5 protein results aren't reported due to technical issues. Expression of ndufs9 is not reported since is used as a housekeeping to normalize the data. Minimum signification level was set at p<0.05. Gene expression and protein content were log-transformed to accomplish the assumptions of normality.

The VDAC content is associated with animal longevity

VDAC is a component of the mitochondrial permeability transition pore (mPTP) and is commonly used as a loading control for mitochondrial proteins. In order to check its suitability when comparing different animals, we've compared total VDAC content by performing an ANOVA analysis. Unexpectedly, we've found that total VDAC content changes across the different animal species and is negatively correlated with animal longevity (Fig. 3A). Furthermore, VDAC total content was positively correlated with NDUFV2 and NDUFS4 protein content (Fig. 3B and C, Supplementary Fig. 1).

Fig. 3

VDAC-1 changes across animal longevity and its association with Cx I subunits and animal longevity. A) Pearson correlation between VDAC-1 protein content and animal longevity. R2(VDAC) = 0.18. B) PGLS regression between VDAC-1 protein content and animal longevity according to phylogenetic tree reported in Fig. 3A. C) Pearson correlation between VDAC-1 and complex 1 subunits. Linear regression (LR) model was performed when significant relationships were found. R2(VDAC vs. NDUFV2) = 0.13; R2(VDAC vs. NDUFS4) = 0.19. Minimum signification level was set at p<0.05. Protein content was log-transformed to accomplish the assumptions of normality. D) Longevity model of longevity modulation via Cx I, ROS and VDAC.

ndufv2 and VDAC correlate with longevity after controlling for phylogenetic relationships

Animal species are evolutionarily related, and closely related species often have similar traits due to inheritance from a common ancestor. Most of statistical analysis, such as linear regression, assume the independence of the data, which might not be accomplished from data obtained from these close-related species. In order to find associations between longevity and gene expression and protein content of selected Cx I subunits and VDAC protein content, we have applied a phylogenetic comparative method, such as phylogenetically generalised least squares (PGLS) regression. A phylogenetic tree allowing to evolutionary relate the species in our study was inferred and is constructed in Supplementary Fig. 2A.

First of all, under the assumption of a Brownian motion model of evolution (a branching, random walk of trait values from an ancestral value at the root to the tips of the tree [30]), we have estimated the Pagels λ. It allows to measure phylogenetic signal and indicates the relative extent to which a traits’ correlation among close relatives match a Brownian motion model of trait evolution. Pagels λ range from 0 to 1, where λ = 1 indicate that trait similarities between species are influenced by phylogenetic relationships; λ = 0 indicate that trait similarities between species are independent of phylogenetic relationships; and 0<λ<1 indicate different levels of phylogenetic signal. According to the estimated λ value, we have classified the measured traits according to its association degree with phylogenetic relationships: i) independent (λ = 0): ndufs5 and NDUFS3; ii) low dependence: ndufs3 and NDUFA9; iii) mild dependence: ndufv2 (λ = 0.65), NDUFV2 (λ = 0.81) and VDAC (λ = 0.53); and strong dependence (λ = 1): ndufs4 and NDUFS4. Finally, we have applied a PGLS regression, which revealed that only gene expression of ndufv2 (p = 0.041, r = −0.80) (Supplementary Fig. 2B) and protein content of VDAC (p = 0.015, r = −0.79) (Supplementary Fig. 2C) were positively correlated with animal longevity after controlling for phylogenetic relationships.

Discussion

Longevity-associated Cx I profile

In this investigation we have found that the mitochondrial Cx I profile is associated with mammalian longevity. Our model revealed that Cx I accounts for 40% of inter-species variation, being ndufv2 gene expression and protein content, followed by NDUFS4 protein, the highest longevity predictor. To the best of our knowledge, ours is the first study elucidating the genetic changes in particular Cx I subunits across mammals of different orders revealing new insights on animal longevity.

Low Cx I content in long-lived mammalian species

Mitochondrial-related genes have been positively selected in long-living invertebrates [31] and mammals, including rodents [32], bats [33] and primates [34]. Our results showed a decreased matrix Cx I subunits gene expression in the same direction that protein content. These findings agree with previous studies that show a positive selection during evolution of gene expression of Cx I subunits, as well as genes that regulate its assembly [35]. Our study design revealed decreased gene expression and/or amount of particular Cx I subunits, suggesting that longevity is associated with a genetic modulation targeted to maintain a lower content of those hydrophilic Cx I matrix domain subunits.

The core Cx I structure is highly conserved across animal species [36] and maintains a fixed stoichiometric relationship of 1:1 in their forming subunits [37]. Since we have found different protein content of the core subunits NDUFV2 and NDUFS4 in different animal species, we suggest that the lowest amount of core subunits determines the total amount of fully assembled Cx I. Accordingly, several studies extrapolated changes in specific Cx I subunits to total Cx I content variations [14,18,38,39]. Furthermore, we are measuring the total content of these proteins, but not its intracellular location. Therefore, we cannot discard that these stoichiometric differences are due to the presence of free subunits at the mitochondrial matrix. Supporting the latter, it has been described different protein turnover rates for each Cx I subunit (up to 4.6-fold difference), which is higher in matrix subunits compared to those located in the membrane [40].

Although ROS can be produced at many different mitochondrial sites, studies including those with endogenous reporters in the absence of inhibitors have localized the main ROS generator sites at the hydrophilic Cx I domain and at Cx III [2,41]. Among these two, the low ROS production of long-lived species and dietary restricted animals has been localized at Cx I [3,[6], [7], [8], [9]], particularly to its matrix hydrophilic domain [10] which catalyses oxidation-reduction reactions, while the hydrophobic membrane domain catalyses proton pumping (Fig. 4). ROS production at the hydrophilic Cx I domain can occur either at the flavin (IF site) or at some of its FeS clusters [[11], [12], [13]], which can also further contribute to secondary ROS formation due to Fe liberation under oxidative stress [1]. Adscription of ROS production to the IF site has been based on experiments increasing or decreasing flavin reduction concomitant with increases and decreases in ROS production respectively. However, these experiments can't rule out Cx I FeS clusters as potential ROS generators because they are located downstream of the flavin, and then increases and decreases in FeS reduction will also occur in those experiments. Studying the amounts of particular Cx I subunits of the hydrophilic Cx I domain can further localize the sites, subunits and genes responsible for the longevity related decrease in mitROSp.

Fig. 4

Low NDUFV2 and NDUFS4 in long-lived species can decrease mitROSp from the hydrophilic domain of Cx I. The figure shows the pathway of electron transfer. The human nomenclature for Cx I subunits (NDUF) is used. Electron transfer is indicated by solid arrows from FMN to Fe S clusters N1a, N3, N1b, N4, N5, N6a, N6b and N2 which reduces the quinone. The conserved off-pathway bifurcation from the flavin to FeS cluster N1a (in NDUFV2) is shown as alternative to main electron transfer to N3. Edge-to-edge distances between cofactors indicated in Å correspond to those described for the mammalian ovine enzyme studied by cryo-Electron Microscopy [20]. The distance between iron-sulphur clusters N1a and N3 (T. thermophilus, Sazanov 2015) is too long for physiological electron transfer to occur either edge-to-edge (19,4 Å) or centre-to-centre (22.3 Å, shown inside brackets). Therefore, the electron transferred from FMN to 2Fe2s cluster N1a could then reduce oxygen to superoxide radical at NDUFV2. A low amount of NDUFV2 in long-lived species would divert less electrons from the main pathway and would thus lead to a lower rate of complex I ROS production. A lower amount of the accessory subunit NDUFS4 in long-lived species can also downregulate ROS generation at N1a through its known interaction with NDUFV2.

Low abundance of NDUFV2 and NDUFS4 subunits is associated with mammalian longevity

The functional role of the Cx I accessory subunits is fragmentary and not clear [21,22] although it has been proposed that they have regulatory roles participating in Cx I stabilization, assembly, subunit interactions or limitation of ROS production [21].

Among the CxI subunits selected, NDUFVS3 and NDUFS5 did not show correlation with species longevity in spite of their decrease in DR and in the long-lived ICRFa mouse strain (moderate decreases), and the increase in NDUFS3 with age in C57Bl/6 mice [10]. The lack of correlation with longevity observed for NDUFS5 is not strange however, since recent description of the near atomic structure of bovine and ovine Cx I at around 4 Å resolution by cryoEM localized NDUFS5 at the bottom of the hydrophobic domain near the intermembrane space [20,21], far away from the Cx I hydrophilic domain sites of ROS production. Concerning NDUFS3, it is an accessory subunit without cofactor that does not take part directly on the main electron transfer but rather provides the platform required for complex assembly [22]. The lack of correlation of the NDUFA9 subunit with species longevity was expected since it did not show correlation in other longevity models [10], it is a subunit essential for activity that binds NADPH to indirectly regulate the terminal N2 cluster [22], and is present at the bottom of the hydrophilic domain near the hinge connecting the two Cx I domains. This protein also has a structural function because it stabilizes the junction between the membrane and matrix domains of Cx I [42].

The remaining two subunits studied, the NDUFV2 and NDUFS4 proteins, showed negative correlation with species longevity. In principle, regulation of ROS production should be left for accessory subunits because the central core is needed for the most fundamental function, electron transfer coupled to H+ pumping and ATP production. However, NDUFV2, which was lower in long-lived species both at mRNA and protein levels, is a special case because it does not take part in the main electron path of the hydrophilic domain of Cx I. Instead, the highly conserved off-pathway N1a FeS cofactor of NDUFV2 can take one electron from FMN but cannot give it to N3 because the distance from N1a (in subunit NDUFV2) to N3 (19.4 to 22.3 Å, Fig. 4) is too long compared to the maximum distance for physiological electron transfer between redox centers, 14 Å [43]. It has been proposed that N1a can function to prevent excessive ROS production [44,45]. The very low one electron potential of N1a (around −370mV) compared to all the other isopotential FeS clusters (N3 to N6b at around −250 mV) makes it highly susceptible to one electron reduction of oxygen to superoxide radical (Fig. 4). The low NUDFV2 (and N1a) amount present in long-lived species agrees with the hypothesis that it can serve to limit ROS production, since the smaller the NDUFV2 amount, the less the number of electrons transferred from FMN to N1a, and thus less would reduce O2 to O2.-, while more electrons would follow the central path from the flavin to N3 and then the rest of the FeS until N2 (Fig. 4). This could help to explain the smaller amount of ROS produced per unit electron flow in long-lived species. This mechanism is simple and avoids the need for electrons to travel bidirectionally between the flavin and N1a [44]. Therefore, we postulate that the amount of NDUFV2 subunit could serve to modulate ROS production in each species in accordance with its longevity and would help to explain the functional significance of the off-path situation of cluster N1a. Interestingly, it has been found that phosphorylation of NDUFV2 tyrosine 118 decreases ROS production leading to cardioprotection after ischemia-reperfusion [46].

The supernumerary non cofactor containing NDUFS4 was the other subunit showing less protein level amount in long-lived species, although its mRNA levels did not show significant differences. Such low level agrees with its strongly decreased levels in long-lived DR mice [10]. The accessory NDUFS4 subunit is located where the NADH dehydrogenase domain meets the rest of the hydrophilic arm [21] and is not likely to play a catalytic role but it is important for assembly of the complex [47]. Various NDUFS4 mutations lead to total loss of Cx I activity [48].

It is known that mammalian NDUFSV1, V3, S1, S4 and S6 interact with NDUFV2 [20,21]. Among all those, only NDUFS4 shows a strong decrease in DR [10] the rest being invariant or showing only modest decreases. Therefore, although NDUFS4 does not have an electron transport cofactor, it could modulate ROS production through its interaction with NDUFV2, which would be affected by the lower NDUFS4 levels of long-lived species. In agreement with this possibility, it is known that cAMP associated serine phosphorylation of NDUFS4 by PKA lowers ROS production and increases Cx I activity [49], and heart-specific NDUFS4 knockout mice have normal heart energetics and function together with less ROS generation and better functional recovery after ischemia-reperfusion while reintroduction of NDUFS4 abolishes this protection [50].

Although our finding of smaller amounts of NDUFV2 and NDUFS4 can contribute to decrease Cx I ROS production in long-lived species, association of complex I within supercomplexes could perhaps also contribute to ROS regulation. Possible functional roles of respiratory supercomplexes remain controversial, proposals including roles in stability of individual complexes, CoQ substrate channeling between Cx I and Cx III2, prevention of protein aggregation, and reduction in ROS production. Various authors have suggested that assembly of Cx I into supercomplexes prevents excessive ROS production during oxidation of NADH-linked substrates because more efficient CoQ channeling helps in maintaining the chain in a more oxidized state [51,52]. However, despite progress in the determination of the structures of respiratory supercomplexes their physiological role is not yet understood [22] and more research is needed to clarify if they could also contribute to modulate mitROSp in relation to animal longevity.

There is an association between mitROSp and VDAC content in animal longevity

ROS production promotes the opening of the mitochondrial permeability transition pore mPTP [25] which is related to mitochondrial swelling and fragmentation, ischemia-reperfusion injury, and cell death [1,50,53]. Specifically, ROS generation and reverse electron transport through Cx I increase the sensitivity of mPTP opening [54]. Previous studies have shown that dietary restriction decreases both mPTP [55,56] and mitROSp at Cx I [3,9]. Since mPTP opening is promoted through oxidation of mitochondrial membrane lipids and proteins by ROS [23], the decrease in mitochondrial permeability in DR could be due to the decrease in mitROSp that occurs in this dietary manipulation.

VDAC is a component of mPTP [[26], [27], [28],57]. Some studies indicate that VDAC content regulates the opening of mPTP in response to modifications in mitROSp [23,25]. Since decreases in Cx I content correlate with reduction in mitROSp [10,[16], [17], [18], [19]], Cx I content could regulate the mPTP opening mediated by mitROSp modifications.

On the other hand, our results show that a lower VDAC content is associated with mammalian longevity. Furthermore, VDAC is positively correlated with NDUFV2 and NDUFS4 Cx I subunits which are present at higher levels in short-lived animals. Thus, we hypothesize that a lower level of NDUFV2 and NDUFS4 would lead to lower complex I ROSp, lower ROS-mediated VDAC modification and lower mPTP in long-lived animal species (Fig. 3). This agrees with the recent description that high VDAC levels lead to mPTP and to a shortened lifespan [29], as well as with the decrease in mPTP observed in dietary restriction studies [55,56].

Dealing with Cx I and VDAC variations in evolutionary-related species

Comparative studies across species with different lifespan are a powerful source of information to identify mechanisms linked to extended longevity [58,59]. However, those studies come up with several technical limitations. Specifically, in terms of protein recognition, the presence of SNPs inducing little sequence or amino acid changes might alter protein structure, even in highly conserved structure, possibly affecting antibody recognition. Since we lack of methods to improve inter-species antibody reactivity, we have increased the number of repeated measurements up to 7, in order to capture higher inter-individual variability in each specie, and measured gene expression using a very sensible technique in order to confirm those results. However, we didn't find significative correlations between gene expression and protein content. Although we cannot omit technical issues, the lack of correlation doesn't necessarily invalidate the presented results. Accordingly, it had been reported that mRNA and protein content doesn't need to be necessarily correlated [60], even in mammals [61]. Furthermore, a comparative study has demonstrated an increasing protein divergence with higher evolutionary distance, since expression divergence seems to be more stable [62], which agrees with the obtained data.

Even after dealing with technical issues, inter-species studies have to face that evolutionary relationships don't allow for independence of the data [63]. Therefore, it's important to elucidate whether a specific trait correlates with longevity differences, or alternatively, these differences arise because of the data similarity. To overcome this limitation, several longevity cross-species studies in birds [64], bats [65], primates [34], mammals [66,67] and vertebrates [68] have applied phylogenetic comparative methods, such as PGLS. In our data, PGLS revealed that ndufv2 gene expression and VDAC protein content changes associated with longevity weren't due interspecies relationships. Although robust, the main limitation of the methods for phylogenetic analysis is that they rely on the construction of a phylogenetic tree, which is unknown. Phylogenies are estimated, mainly by aligning homologue gene sequences and using models of mutation to infer most-likely evolutionary histories. However, depending on the genes chosen, how the sequences are aligned and which method is used to infer evolutionary histories we get different phylogenies. And errors in phylogenetic inferences propagate to errors in phylogenetic analyses [30]. These limitations arise the need of the development of more robust bioinformatic methodologies for phylogenetic tree inferences.

Conclusions

Our study indicates that long-lived mammalian species have low levels of NDUFV2 and NDUFS4 subunits of the hydrophilic domain of complex I which could be responsible in part for their lower mitROSp and superior longevity. The lower VDAC content of these animals could also contribute to their high longevity due to lower ROS-mediated opening of the mitochondrial permeability transition pore.

Methods

Animals

Mammalian species included in the study were male adult specimens with an age representing their 15–30% of their longevity. The recorded values for longevity (in years) were: mouse (Mus musculus), 3.5; rat (Rattus norvegicus), 4.5; gerbil (Meriones unguiculatus), 6.3 guinea pig (Cavia porcellus), 8; rabbit (Oryctolagus cuniculus), 13; pig (Sus scrofa), 27; cow (Bos taurus), 30; and horse (Equus caballus), 46. Rodents and rabbits were obtained from rodent husbandries and sacrificed by decapitation, whereas sheep, pigs, cow, and horses were obtained from abattoirs. The animal care protocols were approved by the Animal Experimentation Ethics committee of the University of Lleida. Heart ventricles from 5-7 animals were removed and immediately frozen in liquid nitrogen and transferred to −80 °C until analyses.

Sample homogenization and quantification

Heart tissue (≈50 mg of whole tissue) was homogenized in a buffer containing 180 mM KCl, 5 mM MOPS, 2 mM EDTA, 1 mM DTPAC adjusted to pH = 7.4. Prior homogenization, 1 μM BHT and a mix of proteases inhibitors (GE80-6501-23, Sigma, Madrid, Spain) and phosphatase inhibitors (1 mM Na3VO4, 1 M NaF) were added. After a brief centrifugation (1000 rpm for 3 min at 4 °C), supernatants protein concentration was measured using the Bradford method (500-0006, Bio-Rad Laboratories, Barcelona, Spain).

Protein content determination

The amount of the specific Cx I subunits NDUFV2, NDUFS3, NDUFS4, NDUFS5 and NDUFA9, as well as VDAC-1/Porin, were estimated using western blot analyses as previously described [39]. Anticipating potential technical limitations mainly due to antibody reactivity, and before performing the comparative approach, we have evaluated the similarity of the protein sequences using a BLAST algorithm (Supplementary Table 1). Although the exact target sequence of the commercial antibodies is unknown, we have found a high similarity between sequences, and have chosen the most appropriate antibodies available according to their proved or predicted reactivity against the species of our study.

Briefly, heart homogenates were mixed with a buffer containing 62.5 mM Tris-HCl pH 6.8, 2% SDS, 10% glycerol, 20% β-mercaptoethanol and 0.02% bromophenol blue, and heated for 3 min at 95 °C. Then, proteins were subjected to one-dimensional electrophoresis on SDS and transferred to PVDF membranes. Membranes were maintained in blocking solution containing Tris 2 M, NaCl 2.5 M, 5% BSA and 0.01% Tween for 1 h at room temperature. Immunodetection was performed using antibodies against NDUFV2 (SAB2107279, Sigma, Madrid, Spain), NDUFS3 (459130, ThermoFisher, Barcelona, Spain), NDUFS4 (ab96549, Abcam, Cambridge, UK), NDUFS5 (ab188510, Abcam, Cambridge, UK), NDUFA9 (459100, ThermoFisher, Barcelona, Spain), and VDAC-1/Porin (ab15898, Abcam, Cambridge, UK).

Secondary antibodies were anti-mouse (GENA931, Sigma, Madrid, Spain) and anti-rabbit (31460, ThermoFisher, Barcelona, Spain). Bands were visualized using an enhanced chemiluminescence HRP substrate (Millipore, MA, USA). Signal quantification and recording was performed with ChemiDoc equipment (Bio-Rad Laboratories, Barcelona, Spain). The protein amount was calculated from the ratio of their densitometry values referred to its respective Coomassie staining (1610436, Bio-Rad Laboratories, Barcelona, Spain) (Supplementary Fig. 3). Due to technical issues (antibody recognition), NDUFA5 amount is not reported. Lack of protein recognition might be due to technical issues, not to missing proteins, such as for NDUFS3 in guinea pig.

Primers design

Gene cDNA sequences coding for the Cx I subunits NDUFV2 (ndufv2), NDUFS3 (ndufs3), NDUFS4 (ndufs4), NDUFS5 (ndufs5) and NDUFA9 (ndufa9) were obtained from Ensmbl (http://www.ensmbl.org). Due to gene cDNA sequences limitations for gerbil, that specie wasn't included in the gene expression analyses. Degenerate primers were designed to amplify conserved regions among mammalian sequences using the software PriFi [69], and are listed in Supplementary Table 2. Primers were purchased from Isogen (LifeSciences, Utretch, Netherlands).

Gene expression: Droplet digital PCR

Prior DNA amplification, RNA from 15 mg whole heart tissue was extracted using RNeasy Fibrous Tissue Mini Kit (Qiagen, Hilden, Germany), and retro-transcribed to cDNA using the High-Capacity cDNA Reverse Transcription kit (Applied Biosystems, CA, USA).

For DNA absolute quantification (copies/μl), reaction mixture contained 1x of EvaGreen ddPCR Supermix, 200 nM primers and 0.01–16 ng of template cDNA. For droplet generation, 20 μL of reaction mixture and 70 μL of Droplet Generation Oil for EvaGreen were loaded in the Droplet Generation Cartridge, which was placed into the Droplet generator. From each PCR reaction mixture, 20 μL were transferred to a 96-well PCR plate, which was sealed with a foil heat using PX1 PCR plate sealer. Amplification was performed in a C1000 Touch Thermal Cycler following an initial DNA Polymerase activation (95 °C, 5 min), and 40 cycles consisting of a DNA denaturation (95 °C, 30 s), primer annealing (58 °C, 1 min) and extension (60 °C, 1 min). A final dye-stabilization step was included (4 °C 5 min, and 90 °C 5 min). Droplets were read with a QX200 Droplet Reader and analysed using QuantaSoft software (Bio-Rad). The results from more than 12,000 droplets were accepted and normalised to an appropriate housekeeping gene (ndufa9) as suggested by previous work [10]. Values are reported as cDNA gene units per cDNA housekeeping units. All equipment and reagents were purchased from Bio-Rad (Bio-Rad Laboratories, Barcelona, Spain).

Statistics

Multivariate statistics was performed using Metaboanalyst software [70]. Pearson correlation and Pearson partial correlation were performed using RStudio (v1.1.453). Linear models and phylogenetic generalised least squares regression (PGLS) were constructed using RStudio (v1.1.453). The phylogenetic tree was constructed using taxa names as described previously [71]. Functions used were included in the package caper. Data was log-transformed and mean-centred prior statistical analyses in order to accomplish the assumptions of linearity.

Funding

This work was supported by the Spanish Ministry of Economy and Competitiveness, Institute of Health Carlos III (grant number PI14/00328), the Spanish Ministry of Science, Innovation and Universities (RTI2018-099200-B-I00), and the Generalitat of Catalonia, Agency for Management of University and Research Grants (2017SGR696) and (SLT002/16/00250) to R.P, and PR [19] BIO MET 0155 to GB.This study has been co-financed by funds from the (“A way to build Europe”). IRBLleida is a CERCA Programme/Generalitat of Catalonia.

Declaration of competing interest

The authors declare no conflict of interest.