A set of unified relativistic mean-field equations of state for hyperonic compact stars recently built in [M. Fortin, Ad. R. Raduta, S. Avancini, and C. Providencia, Phys. Rev. D {bf 101}, 034017 (2020)] is used to study the thermal evolution of non-magnetized and non-rotating spherically-symmetric isolated and accreting neutron stars under different hypothesis concerning proton $S$-wave superfluidity. These equations of state have been obtained in the following way: the slope of the symmetry energy is in agreement with experimental data; the coupling constants of $Lambda$ and $Xi$-hyperons are determined from experimental hypernuclear data; uncertainties in the nucleon-$Sigma$ interaction potential are accounted for; current constraints on the lower bound of the maximum neutron star mass are satisfied. Within the considered set of equations of state, the presence of hyperons is essential for the description of the cooling/heating curves. One of the conclusions we reach is that the criterion of best agreement with observational data leads to different equations of states and proton $S$-wave superfluidity gaps when applied separately for isolated neutron stars and accreting neutron stars in quiescence. This means that at least in one situation the traditional simulation framework that we employ is not complete and/or the equations of state are inappropriate. Another result is that, considering equations of state which do not allow for nucleonic dUrca or allow for it only in very massive NS, the low luminosity of SAX J1808 requires a repulsive $Sigma$-hyperon potential in symmetric nuclear matter in the range $U_Sigma^{(N)}approx 10-30$ MeV. This range of values for $U_Sigma^{(N)} $ is also supported by the criterion of best agreement with all available data from INS and XRT.