Bosonic Gaussian thermal states form a fundamental class of states in quantum information science. This paper explores the information geometry of these states, focusing on characterizing the distance between two nearby states and the geometry induced by a parameterization in terms of their mean vectors and Hamiltonian matrices. In particular, for the family of bosonic Gaussian thermal states, we derive expressions for their Fisher-Bures, Kubo-Mori, and $α$-$z$ information matrices with respect to their mean vectors and Hamiltonian matrices. An important application of our formulas consists of fundamental limits on how well one can estimate these parameters. We additionally establish formulas for the derivatives and the symmetric logarithmic derivatives of bosonic Gaussian thermal states. The former could have applications in gradient descent algorithms for quantum machine learning when using bosonic Gaussian thermal states as an ansatz, and the latter in formulating optimal strategies for single parameter estimation of bosonic Gaussian thermal states. Finally, the expressions for the aforementioned information matrices could have additional applications in natural gradient descent algorithms when using bosonic Gaussian thermal states as an ansatz.