Bayesian inference on large-scale spatio-temporal models is limited by its computational feasibility, a trend that is further exacerbated by the continuous increase in data availability and model refinements. To address this issue, we present a double-layer distributed-memory parallelization strategy for the popular integrated nested Laplace approximations (INLA) method. First, we perform in parallel the different objective function evaluations that happen during the hyperparameter minimization process. Each of these function evaluations requires the assembly and decomposition of large, often structured, sparse precision matrices. A second layer of parallelism is thus found in the decomposition of these sparse matrices. They are handled by a GPU-accelerated, distributed memory, direct solver called Serinv, which we have integrated into our framework. The latter, named pyINLA, is a novel Python package that aims at providing modern coding practices and features for spatio-temporal modeling within the INLA method. PyINLA outperforms the scalability of state-of-the-art packages, allowing the entire framework to run on distributed memory, GPU-accelerated, systems. We showcase the computational performance of our framework on the CSCS Alps supercomputer, scaling further a large spatio-temporal air temperature prediction model made of a spatial mesh composed of 2865 edges over the course of 365 days.