Chalcogenide perovskites are being actively considered for photovoltaic, optoelectronic, and thermoelectric applications due to their high carrier mobility, strong light absorption, long-term stability, and environment-friendliness. For all these applications, thermal properties play a key role in determining the performance and lifetime of perovskite systems. In this work, we have developed a machine-learning Gaussian approximation potential to study the structural and thermal transport properties of chalcogenide perovskite CaZrS3. We show that the GAP achieves a DFT-level accuracy in describing both cubic and orthorhombic CaZrS3, with 2-4 orders of magnitude reduced computational cost. Specifically, we applied the GAP to predict the lattice thermal conductivities (κL) and phonon properties of orthorhombic CaZrS3 from 200 to 900 K by considering four-phonon processes. Compared to its counterpart CaZrSe3, the CaZrS3 exhibits comparably low but relatively more anisotropic κL mainly due to its strong anharmonicity and anisotropic group velocities. Specifically, its thermal conductivities along the a- and c-axis are close and notably lower than that along the b-axis. Optical phonons contribute as high as nearly half of the total thermal conductivity throughout the entire temperature range. Particularly, we observe non-negligible suppression of κL by four-phonon scattering, which is 10% at 300 K and increases to 23% at 900 K. The four-phonon scattering is dominated by the redistribution process, which has a large phase space comparable to that of three-phonon processes within the most frequency range. These results provide a thorough understanding of the phonon transport in orthorhombic CaZrS3 and will be helpful for tailoring the thermal properties and thus the performance of perovskites for potential applications.