In the field of structural optimization, a major objective is to maximize the service life of structures by avoiding unexpected brittle fracture that leads to premature failure. One effective approach is to incorporate tailored heterogeneous materials that can enhance the macroscopic fracture toughness through the interaction between cracks and heterogeneities. While optimization methods that consider stationary cracks with volume and/or stress constraints are well-established, methods that incorporate propagating cracks are still in their early stages of development. Our study proposes a Gaussian processes-based Bayesian optimization strategy to optimize the shapes and placement of stiff elliptical inclusions within a periodic microstructure in 2D, aiming to maximize its effective fracture toughness under mode-I loading. To ensure design feasibility, an inequality constraint is introduced, which requires a minimum clearance between the inclusions. The phase-field fracture method is utilized for bulk crack modeling, and the system potential energy functional is augmented by a softening cohesive potential to capture the effect of tough or weak interfaces. The extended finite element method is used for discretization of the evolving interfaces, which makes it suitable for use in a parametric framework without the need for finite element remeshing. To reduce the computational burden of the fracture simulation, an interior-point monolithic solver is used instead of the standard alternating minimization algorithm for phase-field fracture, along with adaptive mesh refinement near the crack tip and subsequent coarsening along the tail of the crack. In addition, the optimization algorithm considers multiple initial crack locations to account for the possibility of different crack patterns realized for a given design, and optimizes for the worst-case of the corresponding objective values. The numerical experiments demonstrate significant improvements in fracture resistance with justifiable computational cost.