Hydraulic fracture networks (HFN) can be significantly affected by the anisotropy that is caused by bedding planes in shale. In this paper, a transverse isotropic model was built based on the parallel bond model (PBM) and smooth joint model (SJM) in the discrete element method (DEM) to simulate hydraulic fracturing. First, based on the micro characteristics of shale, a discrete fracture network (DFN) was embedded into the model to represent a non-penetrating bedding plane in shale. Second, the model that was constructed was compared with theoretical equations and physical experiments. The comparison results revealed that the model properly matched the breakdown pressure obtained in the theoretical equation as well as the injection pressure, lateral deformation and failure mode from physical experiments. Third, a systematic simulation was implemented to explore the effect of shale anisotropy on HFN under hydraulic fracturing conditions. Finally, sensitivity analysis was performed to examine the effects of the injection rate, initial aperture and permeability coefficient on HFN propagation under different anisotropy angles (θ). The results show that with the increase of θ from the vertical axis, the breakdown pressure increases first and then decreases, and the stimulated reservoir volume ratio (SRVR) shows the opposite trend. The fracturing efficiency (FE) at θ = 0° is far larger than at any other θ. In the crack propagation process, tensile cracks of rock matrix (TCRM) are far more extensive than other crack types, and shear cracks of bedding planes (SCBP) determine the crack orientation. As for the influence factors on HFN propagation, the injection rate has the greatest effect, followed by the permeability coefficient and initial aperture. The proposed model may help predict HFN propagation and optimize the fracturing treatment designs in transversely isotropic shale.