We have extended the canonical tree tensor network (TTN) method, which was initially introduced to simulate the zero-temperature properties of quantum lattice models on the Bethe lattice, to finite temperature simulations. By representing the thermal density matrix with a canonicalized tree tensor product operator, we optimize the TTN and compute accurately the thermodynamic quantities, including the internal energy, specific heat, spontaneous magnetization, etc., at various temperatures. By varying the anisotropic coupling constant $Δ$, we obtain the phase diagram of the spin-1/2 Heisenberg XXZ model on the Bethe lattice, where three magnetic ordered phases, namely the ferromagnetic, XY, and antiferromagnetic phases, are found in low temperatures and separated from the high-T paramagnetic phase by a continuous thermal phase transition at Tc. The XY phase is separated from the other two phases by two first-order phase transition lines at the symmetric coupling points $Δ$=$±$1. We have also carried out a linear spin wave calculation on the Bethe lattice, showing that the low-energy magnetic excitations are always gapped, and find the obtained magnon gaps in very good agreement with those estimated from the TTN simulations. Despite the gapped excitation spectrum, Goldstone-like transverse fluctuation modes, as a manifestation of spontaneous continuous symmetry breaking, are observed in the ordered magnetic phases with |$Δ$|$łeq$1. One remarkable feature there is that the transverse correlation length reaches $ξ$c=1/ln(z-1) for T$łeq$Tc, the maximal value allowed on a z-coordinated Bethe lattice.