We establish shape holomorphy results for general weakly- and hyper-singular boundary integral operators arising from second-order partial differential equations in unbounded two-dimensional domains with multiple finite-length open arcs. After recasting the corresponding boundary value problems as boundary integral equations, we prove that their solutions depend holomorphically upon perturbations of the arcs' parametrizations. These results are key to prove the shape (domain) holomorphy of domain-to-solution maps associated to boundary integral equations appearing in uncertainty quantification, inverse problems and deep learning, to name a few applications.