An extremely useful evolution equation that allows systematically calculating the two-time correlation functions (CFs) of system operators for non-Markovian open (dissipative) quantum systems is derived. The derivation is based on perturbative quantum master equation approach, so non-Markovian open quantum system models that are not exactly solvable can use our derived evolution equation to easily obtain their two-time CFs 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 eq 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 eq L^+.