A regional hybrid gain data assimilation (HGDA) system is newly developed using Weather Research and Forecasting model (WRF). The WRF-HGDA augments an ensemble-based Kalman filter (WRF-LETKF) with information from the variational analysis system (WRF-3DVAR) by combining their gain matrices. The performance of WRF-HGDA is evaluated by assimilating the GNSS radio occultation (RO) observations from the FORMOSAT-3/ COSMIC (FS3/C) and the FORMOSAT-7/COSMIC2 (FS7/C2) under an Observing System Simulation Experiment (OSSE) framework. The results demonstrate that the variational correction improves the WRF-LETKF, with the equal-weighted WRF-HGDA outperforming its component DA systems in the moisture and wind fields when only conventional observations are assimilated. Assimilating additional RO data from FS7/C2 further improves the WRF-LETKF and WRF-HGDA systems. Although the variational correction for the mid-level temperature causes degradation in the WRF-HGDA analysis, this can be alleviated by adjusting the combination weight to include more flow-dependent information in WRF-HGDA at these levels. Further tuning of the static background error covariance used in WRF-3DVAR also brings some improvement in the WRF-HGDA wind analysis. Our results imply that a well-tuned variational system is critical for the accuracy of the regional HGDA analysis.