Symmetry protected topological (SPT) phases in free fermion and interacting bosonic systems have been classified, but the physical phenomena of interacting fermionic SPT phases have not been fully explored. Here, employing large-scale quantum Monte Carlo simulation, we investigate the edge physics of a bilayer Kane-Mele-Hubbard model with zigzag ribbon geometry. Our unbiased numerical results show that the fermion edge modes are gapped out by interaction, while the bosonic edge modes remain gapless at the $(1+1)d$ boundary, before the bulk quantum phase transition to a topologically trivial phase. Therefore, finite fermion gaps both in the bulk and on the edge, together with the robust gapless bosonic edge modes, prove that our system becomes an emergent bosonic SPT phase at low energy, which is, for the first time, directly observed in an interacting fermion lattice model.