We study a class of Stochastic Differential Equations (SDEs) with jumps modeling multistage Michaelis--Menten enzyme kinetics, in which a substrate is sequentially transformed into a product via a cascade of intermediate complexes. These networks are typically high-dimensional and exhibit multiscale behavior with a strong coupling between different components, posing substantial analytical and computational challenges. In particular, the problem of statistical inference of reaction rates is significantly difficult and becomes even more intricate when direct observations of system states are unavailable and only a random sample of product formation times is observed. We address this problem in two stages. First, in a suitable scaling regime consistent with the Quasi-Steady State Approximation (QSSA), we rigorously establish a stochastic averaging principle yielding a reduced model for the product-substrate dynamics. Guided by the reduced-order dynamics, we next construct a novel Interacting Particle System (IPS) that approximates the product-substrate process at the particle level. This IPS plays a pivotal role in the inference methodology, and we prove several non-asymptotic bounds and limiting results for this system. These results facilitate the construction of an estimator based on a product-form approximate likelihood requiring only a random sample of product formation times. This approach does not need access to the system states, and we rigorously prove consistency of the estimator.