Aims. We constrain the space density and properties of massive galaxy candidates at redshifts of z ≥ 3.5 in the Great Observatories Origin Deep Survey North (GOODS-N) field. By selecting sources in the Spitzer + IRAC bands, a sample highly complete in stellarmass is assembled, including massive galaxies that are very faint in the optical/near-IR bands and would be missed by samples selected at shorter wavelengths. Methods. The z ≥ 3.5 sample was selected to m AB = 23 mag at 4.5 μm using photometric redshifts obtained by fitting the galaxies spectral energy distribution at optical, near-IR bands, and IRAC bands. We also require that the brightest band (in AB scale) in which candidates are detected is the IRAC 8.0 μm band to ensure that the near-IR 1.6 μm (rest-frame) peak is falling in or beyond this band. Results. We found 53 z ≥ 3.5 candidates, of masses in the range M ∼ 10 10 −10 11 M . At least ∼81% of these galaxies are missed by traditional Lyman Break selection methods based on ultraviolet light. Spitzer + MIPS emission is detected for 60% of the sample of z ≥ 3.5 galaxy candidates. Although in some cases this might suggest a residual contamination from lower redshift starforming galaxies or Active Galactic Nuclei, 37% of these objects are also detected in the sub-mm/mm bands in SCUBA, AzTEC, and MAMBO surveys, and have properties fully consistent with vigorous starburst galaxies at z ≥ 3.5. The comoving number density of galaxies with stellar masses of above 5 × 10 10 M (a reasonable stellar-mass completeness limit for our sample) is 2.6 × 10 −5 Mpc −3 (using the volume within 3.5 < z < 5), and the corresponding stellar mass density is ∼(2.9 ± 1.5) × 10 6 M Mpc −3 , or about 3% of the local density above the same stellar mass limit. For the subsample of MIPS-undetected galaxies, we measure a number density of ∼0.97 × 10 −5 Mpc −3 and a stellar mass density of ∼(1.15 ± 0.7) × 10 6 M Mpc −3 . Even in the unlikely case that these are all truly quiescent galaxies, this would imply an increase in the space density of passive galaxies by a factor of ∼15 between z ∼ 4 and z = 2, and by ∼100 to z = 0.