A number of observations of stellar systems show a mild preference for anomalously fast cooling compared with what predicted in the standard theory, which leads to a speculation that there exists an additional energy loss mechanism originated from the emission of axions in stars. We revisit the global analysis of the stellar cooling anomalies by adopting conservative assessments on several systematic uncertainties and find that the significance of the cooling hints becomes weaker but still indicates a non-vanishing axion-electron coupling at around 2.4$,sigma$. With the revised analysis results, we explore the possibility that such excessive energy losses are interpreted in the framework of variant axion models, which require two Higgs doublets and flavor-dependent Peccei-Quinn charge assignments. These models resolve two fundamental issues faced in the traditional KSVZ/DFSZ models by predicting a sizable axion coupling to electrons required to explain the cooling anomalies and at the same time providing a solution to the cosmological domain wall problem. We also find that a specific structure of the axion couplings to electrons and nucleons slightly relaxes the constraint from supernova 1987A and enlarges viable parameter regions compared with the DFSZ models. It is shown that good global fits to the observational data are obtained for axion mass ranges of $0.45,mathrm{meV} lesssim m_a lesssim 30,mathrm{meV}$, and that the predicted parameter regions can be probed in the forthcoming helioscope searches.