Functional groups (FGs) are molecular substructures that are served as a foundation for analyzing and predicting chemical properties of molecules. Automatic discovery of FGs will impact various fields of research, including medicinal chemistry and material sciences, by reducing the amount of lab experiments required for discovery or synthesis of new molecules. In this paper, we investigate methods based on graph convolutional neural networks (GCNNs) for localizing FGs that contribute to specific chemical properties of interest. In our framework, molecules are modeled as undirected relational graphs with atoms as nodes and bonds as edges. Using this relational graph structure, we trained GCNNs in a supervised way on experimentally-validated molecular training sets to predict specific chemical properties, e.g., toxicity. Upon learning a GCNN, we analyzed its activation patterns to automatically identify FGs using four different explainability methods that we have developed: gradient-based saliency maps, Class Activation Mapping (CAM), gradient-weighted CAM (Grad-CAM), and Excitation Back-Propagation. Although these methods are originally derived for convolutional neural networks (CNNs), we adapt them to develop the corresponding suitable versions for GCNNs. We evaluated the contrastive power of these methods with respect to the specificity of the identified molecular substructures and their relevance for chemical functions. Grad-CAM had the highest contrastive power and generated qualitatively the best FGs. This work paves the way for automatic analysis and design of new molecules.