We identify various contributors of systematic effects in the measurement of the neutron star (NS) tidal deformability and quantify their magnitude for several types of neutron star-black hole (NSBH) binaries. Gravitational waves from NSBH mergers contain information about the components' masses and spins as well as the NS equation of state. Extracting this information requires comparison of the signal in noisy detector data with theoretical templates derived from some combination of post-Newtonian (PN) approximants, effective one-body (EOB) models, and numerical relativity (NR) simulations. The accuracy of these templates is limited by errors in the NR simulations, by the approximate nature of the PN/EOB waveforms, and by the hybridization procedure used to combine them. In this paper, we estimate the impact of these errors by constructing and comparing a set of PN-NR hybrid waveforms, for the first time with NR waveforms from two different codes, namely, SPEC and SACRA, for such systems. We then attempt to recover the parameters of the binary using two non-precessing template approximants. As expected, these errors have negligible effect on detectability. Mass and spin estimates are moderately affected by systematic errors for near equal-mass binaries, while the recovered masses can be inaccurate at higher mass ratios. Large uncertainties are also found in the tidal deformability., due to differences in PN base models used in hybridization, numerical relativity NR errors, and inherent limitations of the hybridization method. We find that systematic errors are too large for tidal effects to be accurately characterized for any realistic NS equation of state model. We conclude that NSBH waveform models must be significantly improved if they are to be useful for the extraction of NS equation of state information or even for distinguishing NSBH systems from binary black holes.