The quantification of hemoglobin oxygen saturation (sO₂) with multispectral optoacoustic (OA) (photoacoustic) tomography (MSOT) is a complex spectral unmixing problem, since the OA spectra of hemoglobin are modified with tissue depth due to depth (location) and wavelength dependencies of optical fluence in tissue. In a recent work, a method termed eigenspectra MSOT (eMSOT) was proposed for addressing the dependence of spectra on fluence and quantifying blood sO₂ in deep tissue. While eMSOT offers enhanced sO₂ quantification accuracy over conventional unmixing methods, its performance may be compromised by noise and image reconstruction artifacts. In this paper, we propose a novel Bayesian method to improve eMSOT performance in noisy environments. We introduce a spectral reliability map, i.e., a method that can estimate the level of noise superimposed onto the recorded OA spectra. Using this noise estimate, we formulate eMSOT as a Bayesian inverse problem where the inversion constraints are based on probabilistic graphical models. Results based on numerical simulations indicate that the proposed method offers improved accuracy and robustness under high noise levels due the adaptive nature of the Bayesian method.