We formulate a new combined quantum mechanics/molecular mechanics (QM/MM) method based on a self-consistent polarizable embedding (PE) scheme. For the description of the QM region, we apply the popular coupled cluster (CC) method detailing the inclusion of electrostatic and polarization effects into the CC Lagrangian. Also, we consider the transformations required to obtain molecular properties from the linear and quadratic response functions and provide an implementation for the calculation of excitation energies, one- and two-photon absorption properties, polarizabilities and hyperpolarizabilities all coupled to a polarizable MM environment. In the process, we identify CC densitylike intermediates that allow for a very efficient implementation retaining a computational low cost of the QM/MM terms even when the number of MM sites increases. The strengths of the new implementation are illustrated by property calculations on different systems representing the frontier of the capabilities of the CC/MM method. We combine our method with a molecular dynamics sampling scheme such that statistical averages of different excited state solvated properties may be obtained. Especially, we systematically investigate the relative importance of multipoles and polarizabilities in the description of two-photon absorption activity for formamide in aqueous solution. Also, we demonstrate the strengths of the CC hierarchies by incorporating correlation effects both at the CC2, CCSD, and at the triples level in the so-called PE-CCSDR(3) model. Finally, we utilize the presented method in the description of a full protein by investigating the shift of the intense electronic excitation energy of the photoactive yellow protein due to the surrounding amino acids.