Abstract: 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.