  As hosts of living high-mass stars, Wolf-Rayet (WR) regions or WR galaxies
are ideal objects for constraining the high-mass end of the stellar initial
mass function (IMF). We construct a large sample of 910 WR galaxies/regions
that cover a wide range of stellar metallicity, by combining three catalogs of
WR galaxies/regions as previously selected from the SDSS and SDSS-IV/MaNGA
surveys. We measure the equivalent width of the WR blue bump at ~4650 \r{A}
from each spectrum, and make comparisons with predictions of both singular
population models in Starburst99 and binary population models in BPASS. We have
also applied a Bayesian inference code to perform full spectral fitting to the
WR spectra using the singular and binary stellar population models from BPASS
as spectral templates, and we make model selection for models of different IMF
slopes based on the Bayesian evidence ratios. These analyses have consistently
led to a positive correlation of IMF high-mass slope $\alpha$ with stellar
metallicity $Z$, i.e. with steeper IMF (more bottom-heavy) at higher
metallicities, and the conclusion holds even when binary population models are
adopted.
