It is well-known that the upper ocean heat content (OHC) variability in the tropical Pacific contains valuable information about dynamics of El Niño–Southern Oscillation (ENSO). Here we combine sea surface temperature (SST) and OHC indices derived from the gridded datasets to construct a phase space for data-driven ENSO models. Using a Bayesian optimization method, we construct linear as well as nonlinear models for these indices. We find that the joint SST-OHC optimal models yield significant benefits in predicting both the SST and OHC as compared with the separate SST or OHC models. It is shown that these models substantially reduces seasonal predictability barriers in each variable – the spring barrier in the SST index and the winter barrier in the OHC index. We also reveal the significant nonlinear relationships between the ENSO variables manifesting on interannual scales, which opens prospects for improving yearly ENSO forecasting.