An extremely useful evolution equation that allows systematically calculating the two-time correlation functions (CF's) of system operators for non-Markovian open (dissipative) quantum systems is derived. The derivation is based on perturbative quantum master equation approach, so nonMarkovian open quantum system models that are not exactly solvable can use our derived evolution equation to easily obtain their two-time CF's of system operators, valid to second order in the system-environment interaction. Since the form and nature of the Hamiltonian are not specified in our derived evolution equation, our evolution equation is applicable for bosonic and/or fermionic environments and can be applied to a wide range of system-environment models with any factorized (separable) system-environment initial states (pure or mixed). When applied to a general model of a system coupled to a finite-temperature bosonic environment with a system coupling operator L in the system-environment interaction Hamiltonian, the resultant evolution equation is valid for both L = L † and L = L † cases, in contrast to those evolution equations valid only for L = L † case in the literature. The derived equation that generalizes the quantum regression theorem (QRT) to the non-Markovian case will have broad applications in many different branches of physics. We then give conditions on which the QRT holds in the weak system-environment coupling case, and apply the derived evolution equation to a problem of a two-level system (atom) coupled to a finite-temperature bosonic environment (electromagnetic fields) with L = L † .