Polycrystalline materials exhibit complex macroscopic mechanical behavior that originates from collective local interactions of individual grains. Crystal plasticity (CP) frameworks can capture these dependencies by incorporating physically motivated constitutive models that describe various deformation mechanisms. Although CP simulations offer physics-based insights into how microstructural heterogeneity influences global material responses, their high computational cost limit their application in extensive parametric studies to infer structure-property linkages across a diverse spectrum of different microstructures. This challenge motivates the development of efficient, data-driven surrogate models that can rapidly perform predictions without compromising on accuracy. Graph neural networks (GNNs) have emerged as a compelling alternative for modeling such complex systems by representing polycrystalline structures as graphs, where each grain is depicted as a node connected by edges that signify grain boundaries. This graph-based formulation naturally encodes the spatial heterogeneity and connectivity patterns inherent in the material. By embedding grain attributes, such as orientations, grain sizes, and aspect ratios, into the node features, the GNN processes the feature matrix through several message passing layers to aggregate local neighborhood information. Furthermore, GNNs also benefit from their scalability that only depends on the number of grains and, thus, overcomes resolution-dependent limitations of grid-based machine learning methods such as convolutional neural networks. In this work, a large synthetic dataset of microstructure volume elements (MVEs) is generated exhibiting diverse grain size distributions, aspect ratios, and textures. CP simulations with different loading conditions are performed on the dataset. The GNN model is trained to predict the homogenized stress-stain response. Additionally, Bayesian optimization is employed to design microstructural features toward configurations that maximize the stress response. The proposed approach provides a CP surrogate model that significantly reduces computational cost while maintaining the predictive performance, which accelerates materials design and optimization.