While the quantum metrological advantages of performing non‐Gaussian operations on two‐mode squeezed vacuum (TMSV) states have been extensively explored, similar studies in the context of two‐mode squeezed thermal (TMST) states are severely lacking. This paper explores the potential advantages of performing non‐Gaussian operations on TMST state for phase estimation using parity detection‐based Mach–Zehnder interferometry and compares it with the TMSV case. To this end, a realistic photon subtraction, addition, and catalysis model is considered. A unified Wigner function of the photon subtracted, photon added, and photon catalyzed TMST state is derived, which is used to obtain the expression for the phase sensitivity. The results show that performing non‐Gaussian operations on TMST states can enhance the phase sensitivity for significant squeezing and transmissivity parameter ranges. Because of the probabilistic nature of these operations, it is of utmost importance to consider their success probability. When the success probability is considered, the photon catalysis operation performed using a high transmissivity beam splitter is the optimal non‐Gaussian operation. This contrasts with the TMSV case, where photon addition is observed as the most optimal. Further, the derived Wigner function of the non‐Gaussian TMST states will be useful for state characterization and various quantum protocols.