We propose in this paper two kinds of continuity preserving discrete shape gradients of boundary type for PDE-constrained shape optimizations. First, a modified boundary shape gradient formula for shape optimization problems governed by elliptic Dirichlet problems was proposed recently based on the discrete variational outward normal derivatives. The advantages of this new formula over the previous one lie in the improved numerical accuracy and the continuity along the boundary. In the current paper we generalize this new formula to other shape optimization problems including the Laplace and Stokes eigenvalue optimization problems, the shape optimization of Stokes or Navier-Stokes flows, and the interface identification problems. We verify this new formula's numerical accuracy in different shape optimization problems and investigate its performance in several popular shape optimization algorithms. The second contribution of this paper is to propose a continuous discrete shape gradients of boundary type for Neumann problems, by using the ideas of gradient recovery techniques. The continuity property of the discrete boundary shape gradient is helpful in certain shape optimization algorithms and provides certain flexibility compared to the previous discontinuous ones, which are extensively discussed in the current paper.