Melanoma is the deadliest form of skin cancer that causes around 75% of deaths worldwide. However, most of the skin cancers can be cured, especially if detected and treated early. Existing approaches have employed various feature extraction methods, where different types of features are used individually for skin image classification which may not provide sufficient information to the classification algorithm necessary to discriminate between classes, leading to sub-optimal performance. This study develops a novel skin image classification method using multi-tree genetic programming (GP). To capture local information from gray and color skin images, Local Binary Pattern is used in this work. In addition, for capturing global information, variation in color within the lesion and the skin regions, and domain-specific lesion border shape features are extracted. GP with a multi-tree representation is employed to use multiple types of features. Genetic operators such as crossover and mutation are designed accordingly in order to select a single type of features at terminals in one tree of the GP individual. The performance of the proposed method is assessed using two skin image datasets having images captured from multiple modalities, and compared with six most commonly used classification algorithms as well as the standard (single-tree) wrapper and embedded GP methods. The results show that the proposed method has significantly outperformed all these classification methods. Being interpretable and fast in terms of the computation time, this method can help dermatologist identify prominent skin image features, specific to a type of skin cancer in real-time situations.