In QCD with ${N}_{c}$ colors the antisymmetric valence quark color space singlet state $\ensuremath{\sim}{\ensuremath{\epsilon}}_{{i}_{1}\ensuremath{\cdots}{i}_{{N}_{c}}}|{i}_{1},\dots{},{i}_{{N}_{c}}⟩$ of the proton corresponds to the reduced density matrix ${\ensuremath{\rho}}_{ij}=(1/{N}_{c}){\ensuremath{\delta}}_{ij}$ for a single color degree of freedom. Its degenerate spectrum of eigenvalues, ${\ensuremath{\lambda}}_{i}=1/{N}_{c}$, the purity $\mathrm{tr}{\ensuremath{\rho}}^{2}=1/{N}_{c}$, and the von Neumann entropy ${S}_{\mathrm{vN}}=\mathrm{log}({N}_{c})$ all indicate maximal entanglement of color. On the other hand, for ${N}_{c}\ensuremath{\rightarrow}\ensuremath{\infty}$ the spatial wave function of the proton factorizes into valence quark wave functions determined by a mean field [E. Witten, Nucl. Phys. B160, 57 (1979)] where there is no entanglement of spatial degrees of freedom. A model calculation at ${N}_{c}=3$ using a simple three quark model light front wave function by Brodsky and Schlumpf predicts percent level entanglement of spatial degrees of freedom. Using light cone perturbation theory, we also derive the density matrix associated with the four parton $|qqqg⟩$ Fock state. Tracing out the quarks, we construct the reduced density matrix for the degrees of freedom of the gluon, which encodes its entanglement with the sources. Our expressions provide the dependence of the density matrix on the soft cutoff $x$ for the gluon light cone momentum, and on the collinear and ultraviolet regulators. Numerical results obtained in a simple approximation indicate stronger entanglement for the gluon (with ${x}_{g}<⟨{x}_{q}⟩$) than for quarks in the three quark Fock state.