  The full one-loop corrections, both the weak and QED corrections, to the
process $e^+ e^- \to Z \phi $ ($\phi=h^0,H^0$) in the two Higgs doublet model
(2HDM) at the Higgs factories are presented. Up to the $O(\alpha_{em})$ level,
the virtual corrections are evaluated by using the FeynArts/FormCalc packages.
The real emission corrections are computed using the Feynman Diagram
Calculation (FDC) package and the collinear divergences are regularized by the
structure functions of an electron. Using the FeynArts/FormCalc and the FDC
packages, we study the corrections in the Standard Model (SM) and the 2HDM,
respectively. Gauge dependence arising in the normalization of mixing angles is
removed by using the pinch technique. After taking into account experimental
constraints from the current LHC data, we propose four interesting benchmark
scenarios for future colliders. By using these benchmark scenarios, we evaluate
the deviation of $\Delta\sigma(e^+ e^- \to Z \phi)$ from their SM values. We
also examine Higgs boson decays $\phi\to b\bar{b}$ and $\phi\to \tau^+\tau^-$,
which can have large electroweak (EW) contributions from triple Higgs couplings
which are absent in the SM. It is found that for these benchmark scenarios,
both EW and real emission corrections are sizeable and could be measured at a
future $e^+ e^-$ colliders such as the ILC, CLIC, and CEPC.
